2
$\begingroup$

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 = Developer`ToPackedArray[{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 = NDSolve`FiniteDifferenceDerivative[1, nx, opt][
   "DifferentiationMatrix"];
d2udx2 = NDSolve`FiniteDifferenceDerivative[2, nx, opt][
   "DifferentiationMatrix"];
(***************FOR SECOND SPATIAL VARIABLE***************)
dudy = NDSolve`FiniteDifferenceDerivative[1, ny, opt][
   "DifferentiationMatrix"]; 
d2udy2 = NDSolve`FiniteDifferenceDerivative[2, ny, opt][
   "DifferentiationMatrix"]; 
(***************FOR THIRD SPATIAL VARIABLE***************)
dudz = NDSolve`FiniteDifferenceDerivative[1, nz, opt][
   "DifferentiationMatrix"];
d2udz2 = NDSolve`FiniteDifferenceDerivative[2, nz, opt][
   "DifferentiationMatrix"];
(***************FOR FOURTH SPATIAL VARIABLE***************)
dudw = NDSolve`FiniteDifferenceDerivative[1, nw, opt][
   "DifferentiationMatrix"];
d2udw2 = NDSolve`FiniteDifferenceDerivative[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.

$\endgroup$
  • $\begingroup$ Those repeated calls to Interpolation seem costly. Are they truly necessary? $\endgroup$ – mmeent Oct 18 at 14:00
  • $\begingroup$ 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? $\endgroup$ – Fazlollah Oct 18 at 14:10
  • $\begingroup$ Compile will do nothing with high level functions like Interpolation. I all likelihood it will just slow things down. $\endgroup$ – mmeent Oct 18 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.