Solving system of differential equations in loop with recursive variable

Backslide introduced in 11.3, fixed in 12.0.

I've been posting various questions about similar pieces of code all year, so if this looks familiar you may have seen another one of my posts, but this is a unique issue that has not been asked yet, so I am posting it here. Basically I am using a loop so solve a system of coupled differential equations, using those results to change the value of a variable accordingly, then solve the system again, this time with the new value of the variable that was changed, and I want to do this in small steps 200 times, then plot the output. Here is the code:

h2 = {{0, -(Ω/2)}, {-(Ω/2), δ0 + Δ}};

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

ρdecay = {{1/2*γ*ρ22[t], -γ*ρ12[t]}, {-γ*ρ21[t], -(1/2)*γ*ρ22[t]}};

ρtderiv = -I*(h2.ρ2 - ρ2.h2) + ρdecay;

replace3 = {Δ -> -1*10^9, γ -> 1.6*10^9, Ω -> 1, m -> 10^-25, ℏ -> 1*10^-34, k -> (2 π)/(500*10^-9), v -> 10^3};

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

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 π*10^3)/(500*10^-9);

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];
Δt = 0.01;
t0 += Δt;
fscatt = ℏ k^2 γ Re[ρsol22[t0]]/m /. replace3;
δ0 -= Δt fscatt;
txvarray[[i, 1]] = t0;
txvarray[[i, 2]] = fscatt;
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, 3]]]}, {i, 1, 200}];

ListPlot[txabs, PlotStyle -> Green]

Note the δ0 in the first line, in h2. This is the variable that changes after each iteration (this happens in the line δ0 -= Δt fscatt;). The problem is that whenever I try to run this code, it gets stuck and I have to quit the kernal. Any ideas why this is happening/how to fix it? Thanks!

• Cannot reproduce the issue in v9.0.1 and v11.2. (in v9.0.1 the definition of txvarray should be txvarray = Table[{0, 0, 0}, {200}];) though. ) Have you Cleared the variables before executing the code? (Especially δ0.) Please double check the sample. – xzczd Feb 12 at 3:32
• So you're saying it executes fine in v9.0.1 and 11.2? I have tried clearing all variables with this script, but I get the same result. – Caleb Horwitz Feb 12 at 4:42
• Please check if you've posted the correct sample here. – xzczd Feb 12 at 5:07
• A general comment. I suggest breaking up your code into smaller pieces which makes checking the pieces easier. Also use Print[] statements to track the progress of evaluation. – Somos Feb 12 at 19:40
• @xzczd this is definitely the correct sample. I just tried running it with cleared variables again, and I even copy-pasted the code from this post into a new notebook just to see if I messed up. – Caleb Horwitz Feb 13 at 8:37

Update

Just tested on Wolfram cloud, the sample no longer gets stuck in v12.0.

OK, I manage to reproduce the issue in v11.3, but not in v11.2 and v9. I guess this is related to the truth that "CatchMachineUnderflow" option is removed in this version, because after adding WorkingPrecision -> 16 to NDSolveValue together with Rationalize[..., 0], the problem is resolved:

h2 = {{0, -(Ω/2)}, {-(Ω/2), δ0 + Δ}};
ρ2 = {{ρ11[t], ρ12[t]}, {ρ21[t], ρ22[t]}};

ρdecay = {{1/2*γ*ρ22[t], -γ*ρ12[t]}, {-γ*ρ21[t], -(1/2)*γ*ρ22[t]}};

ρtderiv = -I*(h2.ρ2 - ρ2.h2) + ρdecay;

replace3 = {Δ -> -1*10^9, γ -> 1.6*10^9, Ω -> 1,
m -> 10^-25, ℏ -> 1*10^-34, k -> (2 π)/(500*10^-9), v -> 10^3};

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

t0 = 0; ρ120 = 0; ρ210 = 0; ρ220 = 0; ρ110 = 1; δ0 = (2 π*10^3)/(500*10^-9);

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 // Rationalize[#, 0] &, {ρ11, ρ12, ρ21, ρ22}, {t, t0,
t0 + 0.01}, WorkingPrecision -> 16];
Δt = 0.01;
t0 += Δt;
fscatt = ℏ k^2 γ Re[ρsol22[t0]]/m /. replace3;
δ0 -= Δt fscatt;
txvarray[[i, 1]] = t0;
txvarray[[i, 2]] = fscatt;
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, 3]]]}, {i, 1, 200}];

ListPlot[txabs, PlotStyle -> Green]

The output is just the same as in v9 and v11.2:

v11.3 result after adding WorkingPrecision -> 16: v9 result (Definition of txvarray is modified to txvarray = Table[{0, 0, 0}, {200}]): • It works in machine precision for me if I only rationalize the equations. – Michael E2 Apr 10 at 12:54
• @MichaelE2 Which version are you in? In v11.3, Win 7 64bit, WorkingPrecision->16 is necessary, or mxst warning pops up and the kernel crashes. Rationalize[..., 0] isn't necessary, though precw warning will pop up, the result is the same. – xzczd Apr 10 at 13:01
• "11.3.0 for Mac OS X x86 (64-bit) (January 22, 2018)" -- Note the date: the most recent microversion is from March 2018, I think. For some reason, I can't download it from the user portal. – Michael E2 Apr 10 at 13:14
• @MichaelE2 Indeed, mine is "11.3.0 for Microsoft Windows (64-bit) (March 27, 2018)". – xzczd Apr 10 at 13:31
• @xzczd Thank you! This appears to have fixed the issue of the kernal crashing. However, how do I know for certain that it is doing what I want, i.e. how do I know it is in fact using the new value of δ0 and not the first one assigned earlier in the script? – Caleb Horwitz Apr 23 at 4:48