# Numerical solution for ODE using NDSolve

I want to solve the following ODE:

y'[z]==-(y[z]^2-x[z]^2) chi/z^2


with the initial condition

y[z0] == x[z0]


where

x[z_] := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) + 0.43999 E^(-0.965985 z);
chi = 5.5 10^12;
z0 = 20;


I know that the solution, i.e., y(z) should look like:

But instead of this, I got the next solution with Wolfram Mathematica (WM) using NDSolve:

My WM code is:

x[z_] := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) +
0.43999 E^(-0.965985 z);
chi = 5.5 10^12;
z0 = 20;
solution =
NDSolve[{ y'[z] == -(y[z]^2 - x[z]^2) chi/z^2, y[z0] == x[z0]},
y, {z, 10, 100},
Method -> {"EquationSimplification" -> "Residual"}];
LogLogPlot[{Legended[Evaluate[y[z] /. solution],
Placed[StyleForm["y", FontSize -> 12], {.7, .8}]],
Legended[x[z],
Placed[StyleForm["x", FontSize -> 12], {.7, .8}]]}, {z, 100/7,
100}, PlotStyle -> {{Thickness[.004], Red}, {Thickness[.004],
Dashed, Purple}}, FrameLabel -> {"z", "y, x"},
AspectRatio -> 0.95,
PlotRange -> {{100/7, 100}, {0.2 10^-11, 4 10^-7}}, Frame -> True,
PlotStyle -> Thick,
FrameTicksStyle ->
Directive[FontSize -> 12, FontWeight -> Plain, FontColor -> Black],
LabelStyle ->
Directive[FontSize -> 12, FontWeight -> Plain, FontColor -> Black],
AspectRatio -> 0.95]


I do not understand why my solution is getting unstable. Could you help me with this issue?

• You have either not defined yxeq[z] as a function or you've meant to do something different than you do do. Commented Nov 20, 2020 at 20:04
• Sorry, that was a typo, yxeq[z] - > x[z]. I have already corrected this. The basic problem is the correct integration of the ODE. The correct solution is shown in the first picture I included to this question, but I cannot reproduce it in WM. Commented Nov 20, 2020 at 20:45
• Adding AccuracyGoal -> 100 to NDSolve improves y for z>20 Commented Nov 20, 2020 at 20:54
• Thank You! It works! I have to set also z0 = 14 and I got a solution which reasonably well agrees with the correct one. Commented Nov 20, 2020 at 21:41
• I think @MelaGo has a good suggestion, since one probably wants the AccuracyGoal at (64-bit) machine precision to be at least eight digits more than, say, the "average" value of the solution. Commented Nov 20, 2020 at 22:21

Incorporating a number of comments, just so the question has an answer...

x[z_] := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) +
0.43999 E^(-0.965985 z);
chi = 5.5 10^12;
z0 = 14;

solution =
NDSolve[{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2, y[z0] == x[z0]},
y, {z, 10, 100},
Method -> {"EquationSimplification" -> "Residual"},
AccuracyGoal -> 100];

LogLogPlot[{Legended[Evaluate[y[z] /. solution],
Placed[StyleForm["y", FontSize -> 12], {.7, .8}]],
Legended[x[z],
Placed[StyleForm["x", FontSize -> 12], {.7, .8}]]}, {z, 100/7,
100}, PlotStyle -> {{Thickness[.004], Red}, {Thickness[.004],
Dashed, Purple}}, FrameLabel -> {"z", "y, x"},
AspectRatio -> 0.95,
PlotRange -> {{100/7, 100}, {0.2 10^-11, 4 10^-7}}, Frame -> True,
PlotStyle -> Thick,
FrameTicksStyle ->
Directive[FontSize -> 12, FontWeight -> Plain, FontColor -> Black],
LabelStyle ->
Directive[FontSize -> 12, FontWeight -> Plain, FontColor -> Black],
AspectRatio -> 0.95]


• y(z) and x(z) must agree in the region 14<z<20, and then y(z) should smoothly separate from x(z) as shown in the first picture I included to this question. Commented Nov 20, 2020 at 20:49

Originally I took for granted that the expected result was correct. A quick look at the IVP shows it starts on the nullcline $$y(z)=x(z)$$: $$y'(z)=-\frac{\chi\cdot\left(y(z)^2-x(z)^2\right)}{z^2}, \quad y(z_0)=x(z_0) \,.$$ Since $$y'(x)$$ changes sign and $$\chi$$ is large, the quick reversal at $$z=z_0$$ should be expected.

The waviness to the right of $$z=z_0=20$$ in the OP is a numerical issue, but the issue is not stability. The default AccuracyGoal is around 8. Given an accuracy goal $$a$$, error estimates less than around $$10^{-a}$$ are accepted. When the magnitude of the solution is not much more or less than $$10^{-a}$$, then errors are accepted that are significantly large relative to the solution. I think of the accuracy goal as telling NDSolve the largest number that I would consider equivalent to zero. If it produces jittery noise that the user finds unacceptable, the user should increase AccuracyGoal. In this case, since the first run produced a solution that spends a lot of time around $$10^{-11}$$, one should expect to need an AccuracyGoal of around 19 or 20 to get single-precision results.

So two modifications are needed: Change the location of the initial condition to be to the left of the integration interval, and increase the AccuracyGoal.

However, being interested in the mathematical problem, I tried to see how close I could get to x[t]. We can do pretty well if we start at the maximum of x[t]. Since both y'[t] and x'[t] at this point, this is probably the best we can do; at all other initial conditions on the graph of x[t], the solution will cross the graph. The following does a good job until y[z] tries to change sign:

x[z_] = Rationalize[-0.226679 E^(-0.991987 z) -
0.226679 E^(-0.991987 z) + 0.43999 E^(-0.965985 z), 10^-6/2];
chi = 55/10*10^12;
z0 = z /. First@NSolve[x'[z] == 0 && 0 < z < 5, WorkingPrecision -> 32]
solution = NDSolve[
{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2, y[z0] == x[z0]},
y, {z, 0, 100},
AccuracyGoal -> 20,
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}]


NDSolve::ndcf: Repeated convergence test failure at z == 1.1604817342549352; unable to continue.

Changing variables $$t=e^z$$ and using the method "StiffnessSwitching", we can integrate further into y[z] < 0, but the integrator struggles when z < 0.7. The left edge of the plot below represents noise and not the true solution. A log plot is not possible because the solution changes sign, so I scaled the y coordinate by the varying factor Exp[z/2]. Exponential grid lines are shown with the scale indicated on top. The separation of y from x shows y[z] approaching $$\approx 5.15 \times 10^{-12}$$ as $$z \rightarrow +\infty$$.

z0 = z /. First@NSolve[x'[z] == 0 && 0 < z < 5, WorkingPrecision -> 32]
solution =
NDSolve[{Exp[-t] y'[t] == -(y[t]^2 - x[Exp[t]]^2) chi/Exp[t]^2,
y[Log@z0] == x[z0]}, y, {t, -Log@100, Log@100},
AccuracyGoal -> 20, PrecisionGoal -> 10,
Method -> "StiffnessSwitching",
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}, WorkingPrecision -> 30]

Show[
Plot[Evaluate[Exp[z/2] Join[
Table[10^y0/Exp[1/2], {y0, -12, -3}]]],
{z, 0, 43},
PlotStyle -> Join[
Table[Directive[AbsoluteThickness[0.5], Gray], {y0, -12, -3}]],
PlotRange -> {-0.015, 0.015}, Frame -> True,
FrameTicks -> {{Automatic, Automatic}, {Automatic,
Table[{1 - 2 (2 + y0) Log[10], 10.^y0}, {y0, -12, -3, 2}]}}],
Plot[Evaluate[Exp[z/2] Join[
{y[Log@z] /. solution, x[z]}]],
{z, 0, 42},
PlotStyle -> Join[
{AbsoluteThickness[3], AbsoluteThickness[1.5]}],
PlotRange -> All,
PlotLegends -> {y, x}]
]


It should be clear that the desired plot is somewhat misleading. On the one hand, the IC y[z0] == x[z0] always results in a solution that goes down to the left as in the OP. On the other, the plot makes it look like y[z] and x[z] are asymptotic to each other (or even that they agree) as z decreases toward zero.

I think it's thought it was just numerical instability. [The analysis below the graph shows why the OP's plot should have a quick turn at the initial condition.] I had to change the y variable to $$y = 10^{u/10}$$, because StreamPlot is finicky and needs a plot range that is roughly square.

changeVar = y -> (10^(u[#]/10) &);
First@*Solve @@@ (
NestList[
D[#, z] /. First@Solve[First@#, u'[z]] &,
{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2 /. changeVar, u'[z]}, 1]
) /. u[z] -> u // Flatten;
MapThread[Set[#[u_, z_], #2] &, {{up, upp}, {u'[z], u''[z]} /. %}];
?up
?upp

Show[
ContourPlot[
{up[u, z] == 0, upp[u, z] == 0},
{z, 15, 40}, {u, -115, -90},
ContourStyle -> {Directive[Dashed, Darker@Green],
Directive[Dashed, Darker@Magenta]},
FrameTicks -> {{ChartingScaledTicks[{10 # &, #/10 &}],
Automatic}, {Automatic, Automatic}},
PlotLegends -> {u'[z] == 0, u''[z] == 0}
],
StreamPlot[
{1, up[u, z]},
{z, 15, 40}, {u, -115, -90},
StreamPoints -> {{{sepIC + 0.033, Red}, Automatic}}]
]


Moving to the left, if a solution steps across the green line u'[z] == 0, the solution heads down. Step across the purple line u''[z] == 0 and the solution becomes concave down and it will eventually step across the green line.

• I would add that there are clearly solutions that cross u''[z] == 0 and solutions that do not. The boundary between these must also be a solution that does not cross u''[z] == 0. Commented Nov 20, 2020 at 21:34
• Another remark: I didn't pay attention to the ODE. It turns out that u'[z] == 0 is equivalent to y'[z] == 0 and y[z] == x[z]. Commented Nov 20, 2020 at 21:41
• Thank you for your detailed explanation! Setting z0 = 14 and adding AccuracyGoal -> 100 finally solved my problem. Commented Nov 20, 2020 at 21:43
• @student You're welcome. One more remark: The IC y[z0] == x[z0] sets the starting value for y to be the turning point (maximum) of the solution. Since you plot only from 100/7 = 14.2857, you do no see it turn down at z == 14. So it's a good choice for approximating the solution you seek. Perhaps even better is z0 = 10 or 9, since you're integrating from 10 to 100. Commented Nov 20, 2020 at 22:13
• @MichaelE2 Very impressive idea! Using this transformation y->u , NDSolve is able to solve this problem without assumptions AccuracyGoal ->…, Method->… ! Commented Nov 21, 2020 at 12:38

The problem might be solved without additional assumptions AccuracyGoal->…and Method->… using the transformation y[z]->10^u[z] which @MichaelE2 (Thanks for the idea!) showed in his very interesting answer:

x[z_] := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) + 0.43999 E^(-0.965985 z);
z0=10;
chi = 5.5 10^12;

changeVar = y -> (  10^(u[#]/10) &);(**)


transformed ode:

odeu = (-y'[z] - (y[z]^2 - x[z]^2) chi/z^2)     /. changeVar


new boundary condition:

uz0 = u[z0] /. Solve[y[z0] == x[z0] /. changeVar, u[z0]][[1]]


U = NDSolveValue[{odeu == 0, u[z0] == uz0}, u, {z, 10, 100}]

plot of y[z]= 10^(U[z]/10)
LogLogPlot[   10^(U[z]/10)  , {z, z0, 100}, PlotRange -> All]

• Thank you for your other comment. For me, I think z0 = 10; solution = NDSolve[{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2, y[z0] == x[z0]}, y, {z, 10, 100}, AccuracyGoal -> 20] is the simplest approach. (As I mentioned [in another comment](), the goal of 20 is suggested by -Log10[Quartiles[y[Subdivide[10., 100., 100]] /. First@solution]] + 8., which can be easily seen from the plot, too. (+1) Commented Nov 21, 2020 at 13:59
• Yep, but from a test solution, you get a hint about AccuracyGoal. In this case, the OP knew the theoretical asymptotic limit around 5 * 10^-12, which suggests AccuracyGoal -> 20; even if it was unknown, the initial condition x[20.] = 7 * 10^-10 is known. Commented Nov 21, 2020 at 17:04