Several errors with InterpolatingFunction [closed]

I'm getting some error messages that don't make sense. Apologies, as the code is somewhat complex. Basically what I'm trying to do is solve a system of differential equations, then fill a table with values for time and ρsol21[time], with a few other things going on. This is done via a Do loop.

h2 = {{0, -(Ω/2)}, {-(Ω/2), -k v + Δ}}

ρ2 = {{ρ11[t], ρ12[t]}, {ρ21[t], ρ22[t]}}

ρtderiv := -I (h2.ρ2 - ρ2.h2) + {{1/2 γ ρ22[t], -γ ρ12[t]}, {-γ ρ21[t], -(1/2) γ ρ22[t]}}

replace3 = {Δ -> 10, γ -> 1, Ω -> 10, m -> 10^-25, ℏ -> 1.0545718*10^-34, k -> (2 N[Pi, 10])/10^-9};

txvarray = Table[{0, 0, 0}, 200];

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 N[Pi, 10])/10^-10;

Do[
{ρsol11, ρsol12, ρsol21, ρsol22} =
NDSolveValue[{ρ11'[t] == ρtderiv[[1, 1]], ρ12'[t] == ρtderiv[[1, 2]], ρ21'[t] == ρtderiv[[2, 1]], ρ22'[t] == ρtderiv[[2, 2]], ρ11[t0] == ρ110, ρ12[t0] == ρ120, ρ21[t0] == ρ210, ρ22[t0] == ρ220} /. replace3, {ρ11, ρ12, ρ21, ρ22}, {t, t0, t0 + 0.01}];
t0 += 0.01;
δ0 -= 0.01*(2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3);
txvarray[[i, 1]] = t0;
txvarray[[i, 2]] = 2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3;
txvarray[[i, 3]] = δ0;
ρ120 = ρsol12[t0];
ρ210 = ρsol21[t0];
ρ220 = ρsol22[t0];
ρ110 = ρsol11[t0],
{i, 1, 200}]


I'm getting a few error messages that I don't quite understand. The first is

NDSolveValue::mxst: Maximum number of 10000 steps reached at the point t == 0.26002828028276154.


And I don't have any idea what that means. The second is

InterpolatingFunction::dmval: Input value {0.27} lies outside the range of data in the interpolating function. Extrapolation will be used.


This doesn't make sense because when t0=0.27, the domain for the NDSolveValue of that iteration is {0.27,0.28}, so 0.27 should be in the range data. The third message is

An unknown box name (ToBoxes) was sent as the BoxForm for the expression. Check the format rules for the expression.


and I have no idea what this means either. Thanks for your time.

• EDIT: I have identified the source of the third error. It is the line that goes δ0 -= 0.01*(2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3); – Caleb Horwitz Aug 8 '18 at 1:48
• The error message shows that (one of) the integrations stopped at t == 0.26.., which is less than 0.27. That means the interpolating function returned a domain also had a domain that stopped at 0.26+, so that in the update steps, ρ120 = ρsol12[t0], which occurs after t0 = 0.26 is incremented t0 += 0.01 to t0 = 0.27, the input is not in the domain. – Michael E2 Aug 8 '18 at 1:49
• @MichaelE2 why would this happen at t=0.26? In theory the range data should be updated with each iteration so that this never happens, and it doesn't for the first 25 iterations. – Caleb Horwitz Aug 8 '18 at 1:53
• The message says it's because "Maximum number of 10000 steps reached". That can happen because of stiffness, a singularity, or a high number of oscillations. You can plot the 26th solution to see what happened. – Michael E2 Aug 8 '18 at 2:07
• You didn't specify v, so the problem cannot be reproduced and investigated by others. – Michael E2 Aug 11 '18 at 17:51

In principle, this is the working code. There are no mistakes, it is necessary only to define one parameter and avoid numerical divergence. All this is achieved if we put v=10^-9. Then we get a beautiful picture on the output

h2 = {{0, -(Ω/2)}, {-(Ω/2), -k*v + Δ}};

ρ2 = {{ρ11[t], ρ12[t]}, {ρ21[t], ρ22[t]}};

ρtderiv =
-I*(h2.ρ2 - ρ2.h2) + {{1/2*γ*ρ22[t], -γ*ρ12[t]}, {-γ*ρ21[t], -(1/2)*γ*ρ22[t]}};

replace3 =
{Δ -> 10, γ -> 1, Ω -> 10, m -> 10^-25, ℏ -> 1.0545718*10^-34,
k -> (2 N[Pi, 10])/10^-9, v -> 10^-9};

txvarray = Table[{0, 0, 0}, 200];

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = 2*Pi/10^-10;

Do[
{ρsol11, ρsol12, ρsol21, ρsol22} =
NDSolveValue[
{ρ11'[t] == ρtderiv[[1, 1]], ρ12'[t] == ρtderiv[[1, 2]],
ρ21'[t] == ρtderiv[[2, 1]], ρ22'[t] == ρtderiv[[2, 2]],
ρ11[t0] == ρ110, ρ12[t0] == ρ120, ρ21[t0] == ρ210,
ρ22[t0] == ρ220} /.
replace3, {ρ11, ρ12, ρ21, ρ22},
{t, t0, t0 + 0.01},
MaxSteps -> Infinity];
t0 += 0.01;
δ0 -= 0.01*(2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3);
txvarray[[i, 1]] = t0;
txvarray[[i, 2]] = 2 Ω^2 k^2 ℏ Im[ρsol21[t0]/Ω]/m /. replace3;
txvarray[[i, 3]] = δ0;
ρ120 = ρsol12[t0];
ρ210 = ρsol21[t0];
ρ220 = ρsol22[t0];
ρ110 = ρsol11[t0], {i, 1, 200}]

txabs = Table[{txvarray[[i, 1]], Abs[txvarray[[i, 2]]]}, {i, 1, 200}];
ListPlot[txabs, PlotStyle -> Orange]


• I see this got an upvote, but I'm having trouble understanding how this solves the OP's problem. Instead of solving the problem the OP wants to solve (for which v has not been given), your advice seems to be to solve a different problem with a different v chosen to let the integration reach t == 2 (presumably your v is different since the OP's integration fails just past t = 0.26). Is that right? – Michael E2 Aug 11 '18 at 17:48
• Maybe we can just ask the author to give us the secret value of the parameter $v$ at which he received the message and why he did not use the option MaxSteps -> Infinity` and did not investigate the problem of stability? – Alex Trounev Aug 12 '18 at 2:44