# Speeding up the process of NDSolve[] when a user-defined function is involved?

I am trying to tackle a (1+4 dimensional PDE) model at which the solution of the first PDE (with some interpolations and changing the domain) would be used in the second PDE.

In fact, I must choose around $$m_1=m_2=m_3=m_4=20$$ nodes along each dimension while solving the second PDE is challenging and time intensive due to calling a user-defined function named as And2[].

My try for a very few nodes is still time consuming when I solve the second PDE with NDSolve[]. See below:

ClearAll["Global*"];
(***************PARAMETER SETTINGS***************)
TT = 5.;
m1 = 6; m2 = 6; m3 = 6; m4 = 6;
k1 = .05;
size = m1*m2*m3*m4;
roRr = 0; roRz = 0; roRy = 0; royz = 0; sigmaR = 0; gammar = 1.; \
gammaz = -0.2; a = 0;
sigmay = 0.4; sigmar = 0.08; sigmaz = 0.1; rorz = 0; rory = 0;
\[Kappa] = 0.0001; r = 0.02; \[Theta] = -210.; b = 0.1; a1 = 0.08; b1 \
= 0.1; e = 1.15; e1 = 0.45;
(***************DOMAIN DISCRETIZATION***************)
Rmin = 0.; Rmax = 1.; rmin = 0.; rmax = 1.;
ymin = -6.; ymax = 0.; zmin = 0.; zmax = 4.;
pt = DeveloperToPackedArray[{0.45, 0.03, -4.089, 1.15}];
xs = Table[(i - 1.)/(m1 - 1), {i, 1, m1}];
nx = xgrid1 = Range[Rmin, Rmax, (Rmax - Rmin)/(m1 - 1)];
ny = ygrid1 = Range[rmin, rmax, (rmax - rmin)/(m2 - 1)];
nz = zgrid1 = Range[ymin, ymax, (ymax - ymin)/(m3 - 1)];
nw = wgrid1 = Range[zmin, zmax, (zmax - zmin)/(m4 - 1)];
dims = {m1, m2, m3, m4}; grids = {nx, ny, nz, nw}; leng =
Length[origrid];
origrid = Flatten[Outer[List, nx, ny, nz, nw], 3];
origrid1 =
Flatten[Outer[List, nx, ny (1 + gammar), nz, nw (1 + gammaz)], 3];

(***************FILLING SOME MATRICES***************)
Idx = SparseArray[{{i_, i_} -> 1.}, {m1, m1}, 0]; Idy =
SparseArray[{{i_, i_} -> 1.}, {m2, m2}, 0];
Idz = SparseArray[{{i_, i_} -> 1.}, {m3, m3}, 0]; Idw =
SparseArray[{{i_, i_} -> 1.}, {m4, m4}, 0];
DR = KroneckerProduct[(SparseArray@DiagonalMatrix@nx), Idy, Idz, Idw];
Dr = KroneckerProduct[Idx, (SparseArray@DiagonalMatrix@ny), Idz, Idw];
Dy = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@nz), Idw];
Dz = KroneckerProduct[Idx, Idy, Idz, (SparseArray@DiagonalMatrix@nw)];
Id = KroneckerProduct[Idx, Idy, Idz, Idw];

Dy2 = KroneckerProduct[Idx, Idy, (SparseArray@DiagonalMatrix@Exp[nz]),
Idw];
opt = "DifferenceOrder" -> 2;
(***************FOR FIRST SPATIAL VARIABLE***************)
dudx = NDSolveFiniteDifferenceDerivative[1, nx, opt][
"DifferentiationMatrix"];
d2udx2 = NDSolveFiniteDifferenceDerivative[2, nx, opt][
"DifferentiationMatrix"];
(***************FOR SECOND SPATIAL VARIABLE***************)
dudy = NDSolveFiniteDifferenceDerivative[1, ny, opt][
"DifferentiationMatrix"];
d2udy2 = NDSolveFiniteDifferenceDerivative[2, ny, opt][
"DifferentiationMatrix"];
(***************FOR THIRD SPATIAL VARIABLE***************)
dudz = NDSolveFiniteDifferenceDerivative[1, nz, opt][
"DifferentiationMatrix"];
d2udz2 = NDSolveFiniteDifferenceDerivative[2, nz, opt][
"DifferentiationMatrix"];
(***************FOR FOURTH SPATIAL VARIABLE***************)
dudw = NDSolveFiniteDifferenceDerivative[1, nw, opt][
"DifferentiationMatrix"];
d2udw2 = NDSolveFiniteDifferenceDerivative[2, nw, opt][
"DifferentiationMatrix"];

(****************BUILDING THE FIRST SYSTEM MATRIX********************)

B = SparseArray[
+(1/2 sigmaR^2 DR.(Id - DR)).KroneckerProduct[d2udx2, Idy, Idz, Idw]
+ (1/2 sigmar^2 Dr).KroneckerProduct[Idx, d2udy2, Idz, Idw]
+ (1/2 sigmay^2)*KroneckerProduct[Idx, Idy, d2udz2, Idw]
+ (1/2 sigmaz^2 Dz^2)*KroneckerProduct[Idx, Idy, Idz, d2udw2]
+ ((roRr*sigmaR*sigmar)*Sqrt[(DR.(Id - DR)).Dr]).KroneckerProduct[
dudx, dudy, Idz, Idw]
+ ((roRz*sigmaR*sigmaz)*(Dz.Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, Idz, dudw]
+ ((rorz*sigmar*sigmaz)*(Dz.Sqrt[Dr])).KroneckerProduct[Idx, dudy,
Idz, dudw]
+ ((roRy*sigmaR*sigmay)*(Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, dudz, Idw]
+ ((rory*sigmar*sigmay)*(Sqrt[Dr])).KroneckerProduct[Idx, dudy,
dudz, Idw]
+ ((royz*sigmay*sigmaz)*Dz).KroneckerProduct[Idx, Idy, dudz, dudw]
+ (a*(b*Id - DR)).KroneckerProduct[dudx, Idy, Idz, Idw]
+ (a1 (b1*Id - Dr)).KroneckerProduct[Idx, dudy, Idz, Idw]
+ ((r*Id - Dr).Dz).KroneckerProduct[Idx, Idy, Idz, dudw]
+ (\[Kappa] (\[Theta]*Id - Dy)).KroneckerProduct[Idx, Idy, dudz,
Idw]
- (r*Id)
];
B0 = B;

payoff = Flatten@
Table[(1 - nx[[i]]) nw[[l]] (1 + gammaz), {i, 1, m1}, {j, 1,
m2}, {k, 1, m3}, {l, 1, m4}];
initc = Thread[v[0] == payoff];

s = NDSolveValue[{D[v[t], t] == B.v[t], v[0] == initc[[All, 2]]},
v, {t, 0, TT}, Method -> {"FixedStep", "StepSize" -> k1,
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"StiffnessTest" -> False}},
PrecisionGoal -> 5, AccuracyGoal -> 5];

And1[time_Real] :=
Block[{sol2 = s[time]},
Interpolation@Join[origrid1, Partition[sol2, 1], 2]
]
And2[t2_Real] := Quiet@(And1[t2] @@@ origrid)

(****************BUILDING THE SECOND SYSTEM MATRIX********************)

B = SparseArray[
+(1/2 sigmaR^2 DR.(Id - DR)).KroneckerProduct[d2udx2, Idy, Idz, Idw]
+ (1/2 sigmar^2 Dr).KroneckerProduct[Idx, d2udy2, Idz, Idw]
+ (1/2 sigmay^2)*KroneckerProduct[Idx, Idy, d2udz2, Idw]
+ (1/2 sigmaz^2 Dz^2)*KroneckerProduct[Idx, Idy, Idz, d2udw2]
+ ((roRr*sigmaR*sigmar)*Sqrt[(DR.(Id - DR)).Dr]).KroneckerProduct[
dudx, dudy, Idz, Idw]
+ ((roRz*sigmaR*sigmaz)*(Dz.Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, Idz, dudw]
+ ((rorz*sigmar*sigmaz)*(Dz.Sqrt[Dr])).KroneckerProduct[Idx, dudy,
Idz, dudw]
+ ((roRy*sigmaR*sigmay)*(Sqrt[DR.(Id - DR)])).KroneckerProduct[
dudx, Idy, dudz, Idw]
+ ((rory*sigmar*sigmay)*(Sqrt[Dr])).KroneckerProduct[Idx, dudy,
dudz, Idw]
+ ((royz*sigmay*sigmaz)*Dz).KroneckerProduct[Idx, Idy, dudz, dudw]
+ (a*(b*Id - DR)).KroneckerProduct[dudx, Idy, Idz, Idw]
+ (a1 (b1*Id - Dr)).KroneckerProduct[Idx, dudy, Idz, Idw]
+ ((r*Id - Dr).Dz).KroneckerProduct[Idx, Idy, Idz, dudw]
+ (\[Kappa] (\[Theta]*Id - Dy)).KroneckerProduct[Idx, Idy, dudz,
Idw]
- (r*Id)
];

B0 = B;
B = SparseArray[
B0 - Dy2 - gammaz*(Dy2.(Dz.KroneckerProduct[Idx, Idy, Idz, dudw]))];
payoff1 =
Flatten@Table[0., {i, 1, m1}, {j, 1, m2}, {k, 1, m3}, {l, 1, m4}];
initc = Thread[v2[0] == payoff1];

k2 = 4 k1;
fu = ParallelTable[Dy2.And2[t2], {t2, 0, TT, k2}]; // AbsoluteTiming

ss = NDSolveValue[{D[v2[t2], t2] == B.v2[t2] + Dy2.And2[t2],
v2[0] == initc[[All, 2]]}, v2, {t2, 0, TT},
Method -> {"FixedStep", "StepSize" -> k2,
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"StiffnessTest" -> False}},
PrecisionGoal -> 5, AccuracyGoal -> 5]; // AbsoluteTiming


This code takes around 12 second in my system to solve the last NDSolveValue[].

Can anyone suggest how we can speed up the process of NDSolve[] as a time-stepping method, when we call a user-defined function inside it?

I checked for any CUDA applicability or something like that, but I failed. Maybe a compilation of the function And2[] into C or another idea could help.

Thank you.

• Those repeated calls to Interpolation seem costly. Are they truly necessary? – mmeent Oct 18 at 14:00
• Yes, I thought so. But we need to interpolate the previous solution in a new domain and then use it in the process of the second PDE solving. Can Compiling the function help? – Fazlollah Oct 18 at 14:10
• Compile will do nothing with high level functions like Interpolation. I all likelihood it will just slow things down. – mmeent Oct 18 at 14:12