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[(#[]) &, 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*)
0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2), 0]], {i, 1,
n1}, {j, 1, n2}]; // AbsoluteTiming
If[radial <= 1*epsilon, -((
56 Max[cl - radial, 0]^5 (cl + 5 radial) rx[i, j])/cl^8),
0]], {i, 1, n1}, {j, 1, n2}];
If[radial <= 1*epsilon, -((
56 Max[cl - radial,
0]^4 (cl^2 + 4 cl*radial - 5 (radial^2 + 6 rx[i, j]^2)))/
cl^8), 0]], {i, 1, n1}, {j, 1, n2}];
If[radial <= 1*epsilon, -((
56 Max[cl - radial, 0]^5 (cl + 5 radial) ry[i, j])/cl^8),
0]], {i, 1, n1}, {j, 1, n2}];
If[radial <= 1*epsilon, -((
56 Max[cl - radial,
0]^4 (cl^2 + 4 cl*radial - 5 (radial^2 + 6 ry[i, j]^2)))/
cl^8), 0]], {i, 1, n1}, {j, 1, n2}];
If[radial <= 1*epsilon, -((
56 Max[cl - radial, 0]^5 (cl + 5 radial) rz[i, j])/cl^8),
0]], {i, 1, n1}, {j, 1, n2}];
If[radial <= 1*epsilon, -((
56 Max[cl - radial,
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[]) 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]

initc = Thread[u == payoff];
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 == initc[[All, 2]]},
v[t], {t, 0, T}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]; // AbsoluteTiming
{0.197291, Null}

compared to your method:

Monitor[lines =
NDSolve[{eqns, initc}, u[t], {t, 0, T},
Method -> {"EquationSimplification" -> "Residual"},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]; // AbsoluteTiming
{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. – M.J.2 Sep 7 '16 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. – user21 Sep 7 '16 at 15:03
• Your revised comment and answer is very useful. You are right. Now, can you please let me know if there is a way in order to speed up the process of filling those matrices in the code? Because after enforcing your comment, the other tremendous portion of the elapsed time belongs to that! Can we compile the process of filling the matrices in order to speed them up? – M.J.2 Sep 7 '16 at 16:25
• Yes, that should be possible. But I'd suggest to make a new question for that. Simplify the problem, make an example where one needs to change one variable (size) that shows the problem. Get rid of the ParallelTable - fist optimize for your machine then think about parallelization. Make use of compile. – user21 Sep 7 '16 at 16:35