# NDSolve missing dynamics

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]]


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.

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]]


• Thanks for your answer Alex, that works a treat! +1 Commented Jun 21, 2018 at 10:27