# How to accelerate the NDSolve for system of dense ODEs?

I am solving a three-dimensional time-dependent linear PDE with the radial basis function approach. In order not to waste the time of the reader I just include my whole code which run in Mathematica 10.4 quite well when the number of nodes in each direction is low.

ClearAll["Global*"];
t1 = AbsoluteTime[];

(*Setting up the problem values*)
r = 0.04; q = 0; T = 1.; n = 10; size =
n*n*n; e = 100.; σ1 = 0.3; σ2 = 0.35; σ3 = 0.4; \
ρ1 = 0.5; ρ2 = 0.5;
q11 = 0; q22 = 0; q33 = 0; tol = 10^-7; ρ3 = 0.5; beta1 =
1/3; beta2 = 1/3; beta3 = 1/3;
bet = (ρ1*ρ2 - ρ3)/(
1 - ρ1^2); Cmat = {{1, ρ1, ρ2}, {ρ1,
1, ρ3}, {ρ2, ρ3, 1}};
(*Computing the LDLT factorization that we need later*)
LDLT[mat_?SymmetricMatrixQ] :=
Module[{n = Length[mat], mt = mat, v, w},
Do[If[j > 1, w = mt[[j, ;; j - 1]];
v = w Take[Diagonal[mt], j - 1]; mt[[j, j]] -= w.v;
If[j < n, mt[[j + 1 ;;, j]] -= mt[[j + 1 ;;, ;; j - 1]].v]];
mt[[j + 1 ;;, j]] /= mt[[j, j]], {j, n}];
{LowerTriangularize[mt, -1] + IdentityMatrix[n], Diagonal[mt]}];
{lpart, dpart} = LDLT[Cmat];

xsmin = 1.; xsmax = 300.; ysmin = 1.; ysmax = 300.; zsmin = 1.; zsmax \
= 300.;
xgrid1 = Range[xsmin, xsmax, (xsmax - xsmin)/(n - 1)]; ygrid1 =
Range[ysmin, ysmax, (ysmax - ysmin)/(n - 1)]; zgrid1 =
Range[zsmin, zsmax, (zsmax - zsmin)/(n - 1)];
origrid = Flatten[Outer[List, xgrid1, ygrid1, zgrid1], 2];

(*Transformed domain*)
grid = Apply[({(1/σ1) Log[#1/
e], (1/σ2) Log[#2/e] - (ρ1/σ1) Log[#1/e],
(1/σ3) Log[#3/e] + (bet/σ2) Log[#2/
e] - (bet*ρ1 + ρ2) (1/σ1) Log[#1/e]}) &,
origrid, {1}];
xgrid2 = Map[First, grid]; ygrid2 = Map[(#[[2]]) &, grid]; zgrid2 =
Map[Last, grid];
n1 = Length[xgrid2]; n2 = Length[ygrid2]; n3 = Length[zgrid2];

(*The shape parameter inside the RBF approach*)
epsilon =
0.815 (1/size) Total@
Table[Norm[{grid[[i]], grid[[i + 1]]}], {i, 1, size - 1}]
cl = epsilon*T;
rx[k_, l_] := xgrid2[[k]] - xgrid2[[l]];
ry[k_, l_] := ygrid2[[k]] - ygrid2[[l]];
rz[k_, l_] := zgrid2[[k]] - zgrid2[[l]];
rad[k_, l_] := Sqrt[(xgrid2[[k]] - xgrid2[[l]])^2 + (ygrid2[[k]] -
ygrid2[[l]])^2 + (zgrid2[[k]] - zgrid2[[l]])^2]

(*Filling the matrices we need later*)
n1}, {j, 1, n2}]; // AbsoluteTiming
0]], {i, 1, n1}, {j, 1, n2}];
0]^4 (cl^2 + 4 cl*radial - 5 (radial^2 + 6 rx[i, j]^2)))/
cl^8), 0]], {i, 1, n1}, {j, 1, n2}];
0]], {i, 1, n1}, {j, 1, n2}];
0]^4 (cl^2 + 4 cl*radial - 5 (radial^2 + 6 ry[i, j]^2)))/
cl^8), 0]], {i, 1, n1}, {j, 1, n2}];
0]], {i, 1, n1}, {j, 1, n2}];
0]^4 (cl^2 + 4 cl*radial - 5 (radial^2 + 6 rz[i, j]^2)))/
cl^8), 0]], {i, 1, n1}, {j, 1, n2}]; // AbsoluteTiming

phiMatInv = PseudoInverse[phiMat];

(*Initial condition*)
payoff0 = (beta1*Exp[σ1*xgrid2] +
beta2*Exp[σ2 (ygrid2 + ρ1*xgrid2)] +
beta3*Exp[σ3 (zgrid2 - bet*ygrid2 + ρ2*xgrid2)]) - 1;
payoff = Table[Max[payoff0[[i]], 0], {i, 1, size}];

coe = (r - q33 - σ3^2/2)/σ3 +
bet (r - q22 - σ2^2/2)/σ2 - (bet*ρ1 + ρ2) (
r - q11 - σ1^2/2)/σ1;
B = -SparseArray@
Chop[(r*phiMat - ((
r - q11 - σ1^2/2)/σ1) phiMat1x - ((
r - q22 - σ2^2/2)/σ2 - ρ1 (
r - q11 - σ1^2/2)/σ1) phiMat1y
- (coe) phiMat1z - 1/2 phiMat2x -
1/2 (1 - ρ1^2) phiMat2y -
1/2 (dpart[[3]]) phiMat2z).phiMatInv, tol];
u[t_] = Flatten@
Table[Subscript[u, i, j, k][t], {i, 1, n}, {j, 1, n}, {k, 1, n}];
righthand = B.u[t];
eqns = Chop@Thread[D[u[t], t] == righthand];
MatrixPlot[B]

lines = NDSolve[{eqns, initc}, u[t], {t, 0, T},
Method -> {"EquationSimplification" ->
"Residual"}]; // AbsoluteTiming
solution =
Flatten[Table[{xgrid2[[i]], ygrid2[[j]], zgrid2[[k]],
First[Subscript[u, i, j, k][t] /. lines]}, {i, 1, n}, {j, 1,
n}, {k, 1, n}], 2];
list1 = With[{t = T}, Evaluate@solution]; T12 = Map[Last, list1];
set1 = Flatten[
Table[{xgrid1[[i]], ygrid1[[j]], zgrid1[[k]],
e*T12[[n*n*(i - 1) + n (j - 1) + k]]}, {i, 1, n}, {j, 1, n}, {k,
1, n}], 2];

t2 = AbsoluteTime[] - t1;
Print["Whole computational time=", t2];

g = Interpolation@set1;
Print["The final sought-after solution=", g[e, e, e]];


In the above code, for example, when I choose n=10, it means that the size of the discretized system of ODEs is 10*10*10=1000. There are many methods that we can use inside the NDSolve[] for solving time integration system of ODEs. I am using the "MethodSimplification", which works in my code but it is quite slow when $n$ is bigger.

Actually I want to run my code for n=20, but it seems filling the involved matrices and particularly solving the corresponding linear system of ODEs is too time-consuming. So, I would be thankful if someone give me some help in order to handle such dense linear system of ODEs in a quick way.

Two things come to mind. First, you could try to set this up as a matrix equation of the type:

matrix1.u'[t] + stiffness.u[t] == 0


And the other thing you could try is to specify the Jacobian as a pattern SparseArray and see if that helps. Something like

Jacobian -> {Automatic, Sparse-> YourJacSparePattern}


Both are explained in the Finite Element documentation but should work for other discretization methods as well. Have a look here in the section "Transient PDE with Stationary Coefficients and Stationary Boundary Conditions".

Here is a version I'd start from:

Monitor[lines =
NDSolve[{D[v[t], t] == B.v[t], v[0] == initc[[All, 2]]},
v[t], {t, 0, T}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]; // AbsoluteTiming
{0.197291, Null}


Monitor[lines =
{20.153026, Null}

• I wrote my system of ODEs in the given code in the form $u'(t)=B u(t)$. So, it has already been written in the form you mentioned. Furthermore, here my PDE is linear and so there is no Jacobian in viewpoint of nonlinear PDEs as well as since epsilon in my code is changed when $n$ varies, then finding the structure of this matrix is almost impossible. Commented Sep 7, 2016 at 14:41
• @M.J.2, u[t] is a list of {u111[t],.... I think it were more efficient if it were just u[t] (not evaluating to anything) and the Jacobian is still computed, regardless if this is a non-linear problem or not - because the underlying BDF method needs it. Commented Sep 7, 2016 at 15:03