3
$\begingroup$

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: enter image description here

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

enter image description here

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?

$\endgroup$
5
  • $\begingroup$ You have either not defined yxeq[z] as a function or you've meant to do something different than you do do. $\endgroup$
    – Jagra
    Nov 20, 2020 at 20:04
  • $\begingroup$ 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. $\endgroup$
    – student
    Nov 20, 2020 at 20:45
  • $\begingroup$ Adding AccuracyGoal -> 100 to NDSolve improves y for z>20 $\endgroup$
    – MelaGo
    Nov 20, 2020 at 20:54
  • $\begingroup$ Thank You! It works! I have to set also z0 = 14 and I got a solution which reasonably well agrees with the correct one. $\endgroup$
    – student
    Nov 20, 2020 at 21:41
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Nov 20, 2020 at 22:21

3 Answers 3

2
$\begingroup$

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]

enter image description here

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – student
    Nov 20, 2020 at 20:49
5
$\begingroup$

Update: Incorporating comments

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

enter image description here

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.

Original answer

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 -> {{Charting`ScaledTicks[{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}}]
 ]

enter image description here

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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Nov 20, 2020 at 21:34
  • $\begingroup$ 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]. $\endgroup$
    – Michael E2
    Nov 20, 2020 at 21:41
  • $\begingroup$ Thank you for your detailed explanation! Setting z0 = 14 and adding AccuracyGoal -> 100 finally solved my problem. $\endgroup$
    – student
    Nov 20, 2020 at 21:43
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Nov 20, 2020 at 22:13
  • $\begingroup$ @MichaelE2 Very impressive idea! Using this transformation y->u , NDSolve is able to solve this problem without assumptions AccuracyGoal ->…, Method->… ! $\endgroup$ Nov 21, 2020 at 12:38
2
$\begingroup$

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

->ready to solve

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]

enter image description here

$\endgroup$
3
  • $\begingroup$ 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) $\endgroup$
    – Michael E2
    Nov 21, 2020 at 13:59
  • $\begingroup$ Thanks. That's true, but you have to know the solution. $\endgroup$ Nov 21, 2020 at 14:01
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    Nov 21, 2020 at 17:04

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.