Message Boards Message Boards

Why NDSolve is producing unexpected results, independent of some input parameters?

Posted 6 months ago

I am trying to solve a set of coupled differential equation.

  P1 = 10^7;
  P2 = 10^15;
  y3 = 10^-6*(E^-x)*x^(-3/2);
 sol12 = NDSolve[Rationalize[SetPrecision[{y1'[x] + 4 y1[x]/x == 0,
      y2'[x] + (3/
          x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
                P1 ((y2[x] - y3)/y1[x])) - 
              P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
          x] == (-4 P2 (y2[x]^2 - y3^2) - 
          4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
            P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x),
      y1[10] == 10^-7, y2[10] == 0}, 1000]], {y1, y2}, {x, 10, 10^12}, 
   Method -> "Automatic", WorkingPrecision -> 30, 
   MaxSteps -> \[Infinity]]

 p1 = LogLogPlot[{Evaluate[{y1[x]/(y1[x] + y2[x]), y2[x]/(
       y1[x] + y2[x])} /. sol12]}, {x, 10, 10^12}, PlotRange -> Full]

but getting wrong results. The results are also not sensitive for any change in P1 and P2

POSTED BY: John Wick
10 Replies

Here's a way which shows off the "Splitting" and "LocallyExact" methods of NDSolve[]. I rarely use "Splitting" because it requires that one can split the vector field of the ODE $f=f_1+f_2$ in a nice way as well as, in practice, that the ordinary integration methods fail. Consequently I have to look things up to remind myself how it is implemented. It seems easier to me to use DSolve[] as I did in my first reply. But I wanted to remind myself how splitting works. Since "LocallyExact" uses DSolve[] under the hood to advance time integration, I was curious to see if I could get NDSolve[] to do it all for me.

"Splitting" is explained here. I followed an example there. One requirement is that the ODE be autonomous. I used Daniel's substitution $x=e^t$ in the form x == z[t], z'[t] == z[t], z[Log[10]] == 10, which will convert John Wick's nonautonomous system of two DEs into an autonomous system of three DEs. One can then split the vector field into a sum of the form $$\cases{y_1' = -4 y_1 &\cr y_2' = 0 & \cr z' = z &\cr} + \quad\cases{y_1' = 0 &\cr y_2' = *** & \cr z' = 0 &\cr}$$ The first can be integrated with "LocallyExact". For the second, we'll use one-step numerical method. "Splitting" seems to use a fixed step size, which you can set with StartingStepSize, or you can let NDSolve[] estimate it for you (presumably based on the order and error estimates at the beginning of the integration). If you want NDSolve[] to check and reset the step size periodically, I think adding an event like WhenEvent[Mod[t, 1] == 0, {t; "RestartIntegration"}] works. The no-op t seems necessary to trick NDSolve[] into recalibrating.

P1 = 10^7;
P2 = 10^15;
y3 = 10^-6*(E^-x)*x^(-3/2);
ode = {y1'[x] + 4 y1[x]/x == 0, 
   y2'[x] + (3/x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
             P1 ((y2[x] - y3)/y1[x])) - 
           P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) - 
       4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
         P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x)};
ode = Join[
   ode /. {f_'[x] :> f'[t]/z[t], f_[x] :> f[t], x -> z[t]} // 
    Simplify[# // Map[MultiplySides[#, z[t]] &], z[t] > 0] &,
   {z'[t] == z[t]}];
ode = ode /. {f_'[t] :> f'[t], f_[t] :> f[t], t -> z[t]};
(*Print[ode/.v_[t]:>v];*)
With[{sol = First@Solve[ode, {y1'[t], y2'[t], z'[t]}]},
  f1 = {y1'[t] == (y1'[t] /. sol), y2'[t] == 0, 
    z'[t] == (z'[t] /. sol)};
  f2 = {y1'[t] == 0, y2'[t] == (y2'[t] /. sol), z'[t] == 0};
  ];
(*Print[{f1,f2}];*)

sol12 = NDSolve[
  {ode, y1[Log[10]] == 10^-7, y2[Log[10]] == 0, z[Log[10]] == 10}
  , {y1, y2, z}, {t, Log[10], Log[10^12]}
   Method -> {"Splitting", "DifferenceOrder" -> 6,
    "Equations" -> {f1, f2, f1}
    , "Method" -> {
      "LocallyExact"                     (* f1 *)
      , {"ImplicitRungeKutta"            (* f2 *)
       , "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients"
       , "DifferenceOrder" -> 6
       , "ImplicitSolver" -> {FixedPoint, 
         AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, 
         "IterationSafetyFactor" -> 1}}
      , "LocallyExact"}}                 (* f1 *)
  , WorkingPrecision -> MachinePrecision
  (*,InterpolationOrder->All*)
  ]

p1 = LogLogPlot @@ {{y1[t]/(y1[t] + y2[t]), y2[t]/(y1[t] + y2[t])} /. 
     t -> Log[x] /. sol12, Flatten@{x, Exp[y1["Domain"] /. sol12]}, 
   PlotRange -> Full}
POSTED BY: Michael Rogers

Roughly speaking, an AccuracyGoal of $g$ means absolute errors less than $10^{-g}$ are accepted/acceptable. When the magnitude of the solution gets less than this threshold, then errors greater in magnitude may cause the solution to change sign. I think of AccuracyGoal as telling NDSolve the threshold below which a solution may safely be treated as zero plus or minus an error. In this case, you don't want that. Theoretically raising AccuracyGoal so that $10^{-g}$ is much less than the magnitude of the solution is a fix. However, the corrector step in LSODA sometimes/often will exhibit a convergence problem, which it did for me. I could get that to go away by keeping PrecisionGoal modestly low.

So here's the first solution:

sol12a = NDSolve[{y1'[x] + 4 y1[x]/x == 0, 
     y2'[x] + (3/x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - 
             P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) - 
         4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
           P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x), 
     y1[10] == 10^-7, y2[10] == 0}, {y1, y2}, {x, 10, 10^12}, 
    Method -> "Automatic", PrecisionGoal -> 6, 
    WorkingPrecision -> 100, MaxSteps -> \[Infinity]];

Now, the first equation is self-contained subsystem with an exact solution. So here's a second method that avoid the numerical issue entirely:

y1sol = DSolve[{y1'[x] + 4 y1[x]/x == 0, y1[10] == 10^-7}, y1, x];
y2sol = 
  NDSolve[{y2'[x] + (3/x) (((\[Sqrt](y1[x] + 
               y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - 
             P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) - 
         4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
           P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x), 
     y2[10] == 0} /. First@y1sol,
   {y2}, {x, 10, 10^12}, Method -> "Automatic", 
   WorkingPrecision -> 30, MaxSteps -> \[Infinity]];
 sol12b = Join[y1sol, y2sol, 2];

I usually do something like Daniel's approach, but in this case I'd scale the variable y1, which is the misbehaving one. Transforming $x=e^t$ might also be convenient — or not. It often is when a variable ranges from $10$ to $10^{12}$. Here's a third solution:

eqns0 = {y1'[x] + 4 y1[x]/x == 0, 
   y2'[x] + (3/
        x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
             P1 ((y2[x] - y3)/y1[x])) - 
           P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
       x] == (-4 P2 (y2[x]^2 - y3^2) - 
       4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
         P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x)};
eqns1 = eqns0 /. y1 -> Exp@*u1 // Simplify[#, u1[x] > 0] &;

sol12c = NDSolve[Join[eqns1, {u1[10] == Log[10^-7], y2[10] == 0}],  
    {u1, y2}, {x, 10, 1*^12}] /. (u1 -> usol_) :>  y1 -> Exp@*usol;
POSTED BY: Michael Rogers
Posted 6 months ago

Unfortunately, at the second step you have replaced x by t instead of Exp[t] and y'[x] -> y2'[t]*Exp[-t], after correcting these minor typos the solution remains the same as previous, i.e. negative valued at some places and independent of P1 and P2 enter image description here

POSTED BY: John Wick

Thanks for those corrections. So now it seems one issue is that NDSolve has trouble even with just the y1 ODE.

y1val = NDSolveValue[{(4 y1[t])/Exp[t] + E^-t Derivative[1][y1][t] == 
     0, y1[Log[10]] == 10^(-7)}, y1[t], {t, Log[10], 12*Log[10]}];

This will show it dips below zero in several places.

Plot[y1val, {t, 5 Log[10], 12*Log[10]}, PlotRange -> Full]

I don't know if this is a bug or some consequence of solving ODEs numerically. I get a decent result if I set Method->"Adams" but that does not fully carry over to the full system. It shows improvement but eventually the y1 value breaks down into what might be noise.

Possibly some different NDSolve Method settings and subsettings would do better.

POSTED BY: Daniel Lichtblau
Posted 6 months ago

Theoretically the coupled equations are devised such that the y1, y2 should always remain positive and, with increasing x at first y2 should increase and dominate over y1 but at higher x, y1 should become higher than y2. I also have tried to keep the P1, P2 and initial condition practical.

   P1 = 10^7;
   P2 = 10^12;
   y3 = 10^-6*BesselK[2, x]*x^(-3/2);

   sol12 = N[DSolve[Rationalize[SetPrecision[{y1'[x] + 4 y1[x]/x == 0,
         y2'[
            x] + (3/x) (((\[Sqrt](y1[x] + 
                   y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
                   P1 ((y2[x] - y3)/y1[x])) - 
                 P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
             x] == (-4 P2 (y2[x]^2 - y3^2) - 
             4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x),
         y1[0.1] == 10^7, y2[0.1] == 0.10}, 1000]], {y1, y2}, {x, 0.1,
        10^10}, 
      Method -> {"StiffnessSwitching", "ExtraPrecision" -> 1000, 
        InterpolationOrder -> \[Infinity]}], {25, 
      500000}](*,WorkingPrecision\[Rule]50,MaxSteps\[Rule]\[Infinity],\
   InterpolationOrder\[Rule]All]*)

by evaluating this I am getting at higher x the issue

POSTED BY: John Wick

I don't understand the SetPrecision. Your equations have infinite precision, except for the 0.1 in the initial data, which are remedied by writing 1/10. The SetPrecision will decrease the precision from infinnite to finite. But then Rationalize will probably revert to the infinite precision where you started.

I think that N[..., {25, 500000}] has no effect at all.

POSTED BY: Gianluca Gorni

I get more stable results by rescaling the horizontal axis using the substitution t=Log[x]. I believe I did it correctly below.

eqns0 = {y1'[x] + 4 y1[x]/x == 0, 
   y2'[x] + (3/
        x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
             P1 ((y2[x] - y3)/y1[x])) - 
           P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
       x] == (-4 P2 (y2[x]^2 - y3^2) - 
       4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
         P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x)};
eqns1 = eqns0 /. {y1'[x] -> y1'[t]*Exp[t], 
     y2'[x] -> y2'[t]*Exp[t]} /. x -> t;

sol12b = 
 NDSolve[Join[eqns1, {y1[Log[10]] == 10^-7, y2[Log[10]] == 0}], {y1, 
   y2}, {t, Log[10], 12*Log[10]}]

Possibly the LogPlot (not log-log since we already transformed one axis) will be more along the lines of what you expect.

POSTED BY: Daniel Lichtblau

Can you explain why you find the numerical solution "wrong"?

POSTED BY: Gianluca Gorni

Your numerical solution y1 has small fluctuations around 0. In the intervals where y1 is negative, the LogLogPlot will leave a gap. P1 and P2 do not affect y1. Luckily, the differential equation for y1 can be solved symbolically:

DSolve[{y1'[x] + 4 y1[x]/x == 0, y1[10] == 10^-7}, y1, x]
POSTED BY: Gianluca Gorni
Posted 6 months ago

I don't think the problem is due to my equations, but somehow, they are coming from NDSolve. Very similar to Error control for NDSolve. But I don't know how to resolve this.

POSTED BY: John Wick
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract