# Why is this non-vectorized NDSolve faster than vectorized?

I'm trying to get a feeling for stochastic population models. A simple starting place is a birth-death process whose deterministic limit is the logistic equation. The idea is to model the probability p[n][t] of having n individuals at time t, which results in an infinite set of linear differential equations (truncated at nmax to reduce it to a large but finite set).

I have implemented this two ways: non-vectorized and vectorized. I thought the vectorized implementation would improve the speed, but it is actually much slower. Does anyone know why, or have other tips for speeding this up? The code below is fast enough in both cases, but it gets worse for larger nmax.

(* set up parameters *)
nmax = 120; (* max population size*)
k = 60; (* carrying capacity *)
b0 = 2; (* density-independent birth rate *)
d0 = 1; (* density-independent death rate *)
d1 = (b0 - d0) / k; (* density-dependent death rate *)

Do[
b[n] = b0*n; (* birth rate *)
d[n] = (d0 + d1*n)*n; (* death rate *)
, {n, 0, nmax}];

tmax = 10^9;  (* max time *)


Here's the non-vectorized version:

(* set up equations, ICs, and unknowns *)
eqns = Join[Table[
p[n]'[t] == b[n - 1] p[n - 1][t] - b[n] p[n][t] + d[n + 1] p[n + 1][t] - d[n] p[n][t]
, {n, 1, nmax - 1}],
{p'[t] == d p[t],
p[nmax]'[t] == b[nmax - 1] p[nmax - 1][t] - d[nmax] p[nmax][t]}
];
ics = Table[p[n] == If[n == k, 1, 0], {n, 0, nmax}];
unks = Table[p[n], {n, 0, nmax}];

(* solve it *)
AbsoluteTiming[
sol = NDSolve[Flatten[{eqns, ics}], unks, {t, 0, tmax}][];
]

(* {0.028615, Null} *)


Here's the vectorized version that I thought would be faster:

i[n_] := n + 1; (* add 1 to indices in array *)

(* set up transition matrix *)
a = Normal@SparseArray[Join[
Table[{i[n], i[n]} -> -b[n] - d[n], {n, 0, nmax}],
Table[{i[n + 1], i[n]} -> b[n], {n, 0, nmax - 1}],
Table[{i[n - 1], i[n]} -> d[n], {n, 1, nmax}]
]];

(* set up ICs *)
p0 = Normal@SparseArray[i[k] -> 1, {nmax + 1}];

(* solve it *)
AbsoluteTiming[
sol2 = NDSolve[{p'[t] == a.p[t], p == p0}, p, {t, 0, tmax}][];
]

(* {0.586452, Null} *)


Both give the same results. For fun, here they are:

LogLogPlot[Evaluate[Table[p[n][t], {n, nmax, 0, -1}] /. sol], {t, 0.001, tmax}] ListPlot[Table[{n, p[n][10^5]}, {n, 0, nmax}] /. sol, PlotRange -> All]
ListPlot[Table[{n, p[n][10^9]}, {n, 0, nmax}] /. sol, PlotRange -> All]  You can see that the system quickly reaches a quasi-stationary distribution around the carrying capacity of n=60, but over a long timespan reaches the true stationary distribution of everyone dead (p==1).

Update:

I tried out some NDSolve options and found two major improvements for the vectorized code that put it on par with the non-vectorized. First, remove the Normal from the definition of a as per @user21's comment. Then try some different options:

a = SparseArray[
Join[Table[{i[n], i[n]} -> -b[n] - d[n], {n, 0, nmax}],
Table[{i[n + 1], i[n]} -> b[n], {n, 0, nmax - 1}],
Table[{i[n - 1], i[n]} -> d[n], {n, 1, nmax}]]];

AbsoluteTiming[
sol2 = NDSolve[{p'[t] == a.p[t], p == p0}, p, {t, 0, tmax},
Jacobian -> {Automatic, Sparse -> a["PatternArray"]}
][];]

(* {0.043259, Null} *)

AbsoluteTiming[
sol2 = NDSolve[{p'[t] == a.p[t], p == p0}, p, {t, 0, tmax},
Jacobian -> {a}
][];]

(* {0.16169, Null} *)

AbsoluteTiming[
sol2 = NDSolve[{p'[t] == a.p[t], p == p0}, p, {t, 0, tmax},
Method -> {"IDA", "ImplicitSolver" -> {"Newton",
"LinearSolveMethod" -> {"Band", "BandWidth" -> {1, 1}}}}
][];]

(* {0.034921, Null} *)


Looks like the first and third versions are winners. Too bad they don't seem to stack:

AbsoluteTiming[
sol2 = NDSolve[{p'[t] == a.p[t], p == p0}, p, {t, 0, tmax},
Method -> {"IDA", "ImplicitSolver" -> {"Newton",
"LinearSolveMethod" -> {"Band", "BandWidth" -> {1, 1}}}},
Jacobian -> {Automatic, Sparse -> a["PatternArray"]}
][];]

(* {0.033892, Null} *)

• Why do do call Normal on the SparseArray? Feb 10, 2017 at 14:47
• @user21 Without the Normal on the definition of the ICs p, I get the error NDSolve::ndinnt: Initial condition SparseArray[Automatic,{121},0.,{1,{{0,1},{{61}}},{1.}}] is not a number or a rectangular array of numbers. Leaving the Normal on the ICs but removing it from a runs, but still takes 0.317334 AbsoluteTiming. Feb 10, 2017 at 14:53
• Because Length@Flatten@a/ Length@Flatten[{eqns, ics}]=121/2? Feb 10, 2017 at 15:55
• @Feyre Not sure I follow your comment. What you say is true, but a is the transition matrix (mostly zeros). The size of the problem seems identical in both cases, because Length@Flatten@p0==Length@Flatten@ics==Length@Flatten@eqns==121. Feb 10, 2017 at 16:05
• Similar phenomenon here, not sure if the nature is the same. Feb 11, 2017 at 4:19