2
$\begingroup$

I am trying to use NDSolve to solve some a set ODEs coupled in a slightly slightly unusual, and am running into some problems.

The following code generates the equations of interest, where one the solution of one ODE is used as the initial condition for a second, generating a collection of 2D interpolating functions (apologies for the large matrix, this is the minimal example I could think of!):

(*Prepare vector of dependent variables*)
Rho[t_]= Flatten@Array[Subscript[r,#1,#2][t]&,{3,3}];
Lamb[t_, tau_]=Flatten@Array[Subscript[La,#1,#2][t,tau]&,{3,3}];

(*Include propagator matrix*)
lio ={{-0.0000506308, 0. + 1. I, 0, 0. - 1. I, 0, 0, 0, 0, 0}, {0. + 1.I, -9.72405, 0, 0, 0. - 1. I, 0, 0, 0, 0}, {0,  0, -4.72405, 0, 0, 0. - 1. I, 0, 0, 0}, {0. - 1. I, 0, 0, -9.72405,   0. + 1. I, 0, 0, 0, 0}, {0, 0. - 1. I, 0, 0. + 1. I, -10., 0, 0, 0,   0}, {0, 0, 0. - 1. I, 0, 0, -5., 0, 0, 0}, {0, 0, 0, 0, 0,   0, -4.72405, 0. + 1. I, 0}, {0, 0, 0, 0, 0, 0, 0. + 1. I, -5.,   0}, {0.0000506308, 0, 0, 0, 10., 0, 0, 0, 0}};

(*Solve the first ODE*) 
rho0 = {1, 0, 0, 0, 0, 0, 0, 0, 0};
Rsol = First@NDSolve[{D[Rho[t],t]==lio.Rho[t], Rho[0]==rho0}, Rho[t],{t,0,100}];

(*Solve the second ODE using the first as initial condition*) 
Lam0 = {0., 0., 0., 0., 0., 0., 0. + 1. Subscript[r, 2, 1][t], 0. + 1. Subscript[r, 2, 2][t], 0. + 1. Subscript[r, 2, 3][t]}/.Rsol;
LambSol = NDSolveValue[{D[Lamb[t,tau],tau]==lio.Lamb[t,tau], Lamb[t,0] == Lam0},Lamb[t,tau],{t,0,100},{tau,0,3}];

The issue arises when we compare the solutions in the t variable. Plotting RSol[[5]] and the matching element LambSol[[-2]], with tau = 0, we find:

Show[LogLinearPlot[{Evaluate[Rho[t][[5]]/.Rsol],Evaluate[LambSol[[-2]]/.tau->0]},{t,10^-3,20}],
ListLogLinearPlot[Transpose[{LambSol[[-2,0]]["Coordinates"][[1]],
LambSol[[-2,0]]["ValuesOnGrid"][[All,1]]}],PlotStyle->Red]]

enter image description here

As you can see the dynamics found from the second ODE solution (yellow) doesn't match the converged solution from RSol (blue). I think the reason is that NDSolve is not taking enough points the along the t-axis in second NDSolve (given by the red points). I can fix this to an extent by decreasing MaxStepFraction, however taking this too small or for systems of equations that are too large makes the code take excessive amounts of time, or the Kernel crash.

Does anybody have insight into using NDSolve in this way? I'm at a loss of how to tackle this problem.

$\endgroup$
2
$\begingroup$

You are right that with decreasing the integration step, the results get closer.Here is an example with a step of 0.1 - the left figure.On the other hand, if the step is reduced to 0.05, the integration time increases noticeably and a memory of 13940 MB is required - the second figure.I see the way out is to go to the logarithm of time - see the code and the bottom pictures with an automatic step and with a step of 0.1.

    (*Prepare vector of dependent variables*)
Rho[t_] = Flatten@Array[Subscript[r, #1, #2][t] &, {3, 3}];
Lamb[t_, tau_] = 
  Flatten@Array[Subscript[La, #1, #2][t, tau] &, {3, 3}];

(*Include propagator matrix*)
lio = {{-0.0000506308, 0. + 1. I, 0, 0. - 1. I, 0, 0, 0, 0, 
    0}, {0. + 1. I, -9.72405, 0, 0, 0. - 1. I, 0, 0, 0, 0}, {0, 
    0, -4.72405, 0, 0, 0. - 1. I, 0, 0, 0}, {0. - 1. I, 0, 
    0, -9.72405, 0. + 1. I, 0, 0, 0, 0}, {0, 0. - 1. I, 0, 
    0. + 1. I, -10., 0, 0, 0, 0}, {0, 0, 0. - 1. I, 0, 0, -5., 0, 0, 
    0}, {0, 0, 0, 0, 0, 0, -4.72405, 0. + 1. I, 0}, {0, 0, 0, 0, 0, 0,
     0. + 1. I, -5., 0}, {0.0000506308, 0, 0, 0, 10., 0, 0, 0, 0}};

(*Solve the first ODE*)
rho0 = {1, 0, 0, 0, 0, 0, 0, 0, 0};
Rsol = First@
   NDSolve[{D[Rho[t], t] == lio.Rho[t]*Exp[t], Rho[-5] == rho0}, 
    Rho[t], {t, -5, 5}, MaxStepSize -> 0.1];

(*Solve the second ODE using the first as initial condition*)
Lam0 = {0., 0., 0., 0., 0., 0., 0. + 1. Subscript[r, 2, 1][t], 
    0. + 1. Subscript[r, 2, 2][t], 0. + 1. Subscript[r, 2, 3][t]} /. 
   Rsol;
LambSol = 
  NDSolveValue[{D[Lamb[t, tau], tau] == lio.Lamb[t, tau]*Exp[tau], 
    Lamb[t, -5] == Lam0}, Lamb[t, tau], {t, -5, 5}, {tau, -5, 5}, 
   MaxStepSize -> .1];
Show[Plot[{Evaluate[Rho[t][[5]] /. Rsol], 
   Evaluate[LambSol[[-2]] /. tau -> -5]}, {t, -5, 5}], 
 ListPlot[Transpose[{LambSol[[-2, 0]]["Coordinates"][[1]], 
    LambSol[[-2, 0]]["ValuesOnGrid"][[All, 1]]}], PlotStyle -> Red]]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for your answer Alex, that works a treat! +1 $\endgroup$ – Jiles Jun 21 '18 at 10:27

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.