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*)
phiMat = SparseArray@ParallelTable[With[{radial = rad[i, j]},
If[radial <= 1*epsilon,
Max[(1 - radial/cl),
0]^6 (3 + 18 radial/cl + 35 (radial/cl)^2), 0]], {i, 1,
n1}, {j, 1, n2}]; // AbsoluteTiming
phiMat1x = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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}];
phiMat2x = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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}];
phiMat1y = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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}];
phiMat2y = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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}];
phiMat1z = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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}];
phiMat2z = SparseArray@ParallelTable[With[{radial = rad[i, j]},
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[[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]
initc = Thread[u[0] == 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.
