Thanks for paying attention to my question and modifying my terrible formatting[facepalm], now I think I've found the answer. The point is the selection StopIntegration should be modified, for example, WhenEvent[z[ρ]=10` Sqrt[zt^2 - ρc^2], the problem can be solved.

The original question is as below:

I'm solving a differential equation to get a function z[$\rho$], the expected solution needs to start from $z[\rho_c]=z[\frac{13}{5}]=\frac{24195890215702774518563237183870329}{8112963841460668169578900514406400}\approx 2.98$, until when the function $z[\rho]$ hits on the $\rho$-axis. The code is as follows:

precTmp = 4 MachinePrecision;
R = 4;
sol = 
  (ρc = 13/5;
   vt = Sqrt[zt^2 - ρc^2] - zt; 
     {(1 - z[ρ]^3)^2/z[ρ] + 3/2 z[ρ]^2 Derivative[1][z][ρ]^2 + 
        (1 - z[ρ]^3)(Derivative[1][z][ρ]^2/z[ρ] + (z^′′)[ρ]) + 
        (169 (-(13/5) + zt)^2 (13/5 + zt)^2 z[ρ]^2 (-3 z[ρ]^2 + 2 (z^′′)[ρ])) / 
        (200 zt^2) == 0, 
      Derivative[1][z][13/5] == 
        -((13 (2 - (-(169/25) + zt^2)^(3/2)))/(10 Sqrt[-(169/25) + zt^2])), 
      z[13/5] == Sqrt[-(169/25) + zt^2],
      WhenEvent[z[ρ] == 10^-3, "StopIntegration"], 
      WhenEvent[z[ρ] == 1.1` Sqrt[zt^2 - ρc^2], "StopIntegration"]}, 
     z, {ρ, ρc, R + 1}, {zt}, 
     WorkingPrecision -> precTmp, 
     AccuracyGoal -> precTmp/2, 
     MaxSteps -> 10^6, 
     Method -> "StiffnessSwitching"]);

Note that we can adjust the parameter precTmp to ramp up the precision. To check the correctness of the solution, it's natural to plot the solution out, as follows:

ztValue = 
  {ρ, sol[ztValue]["Domain"][[1, 1]], sol[ztValue]["Domain"][[1, 2]]}, 
  PlotRange -> All, 
  WorkingPrecision -> precTmp]

The first thing I can't understand is that the shape of the solution $z[\rho]$ is different when I change the precision, for example, taking precTmp to be 3 MachinePrecision, 5 MachinePrecision, or 6 MachinePrecision, the solution is decreasing and hits the $\rho$-axis, but when taking precTmp to be 4 MachinePrecision, 8 MachinePrecision or 12 MachinePrecision, the solution turns out to be increasing.

My first question is which one of them is the correct answer?

Another thing I can't understand is when I provide a PrecisionGoal to the ParametricNDSolveValue function — for example, when I set precTmp = 4 MachinePrecision and give PrecisionGoal -> precTmp — the solution becomes a decreasing function, although it was originally an increasing one. I've already learned the difference between precision and accuracy, which can somehow be regarded to be equivalent as long as the function doesn't go to zero.

My second question is how to understand the difference brought by varying PrecisionGoal.

  • $\begingroup$ For R = \[Rho]c, I get the same solutions, but with different domains, which make the plots look different, unless you superimpose them. (Maybe should includ R in the question.) $\endgroup$
    – Michael E2
    Commented Jan 18, 2020 at 1:31
  • $\begingroup$ Ahh, I'm sorry, the parameter $R$ is 4. $\endgroup$
    – Xuao Zhang
    Commented Jan 18, 2020 at 9:34
  • $\begingroup$ You should put your answer in an answer, not in the question. Answering one's own question is explicitly allowed and encouraged. $\endgroup$
    – Michael E2
    Commented Jan 18, 2020 at 11:23
  • $\begingroup$ I see, thanks for your suggestion, let me write the answer in a better form in case someone makes similar mistake as I did. $\endgroup$
    – Xuao Zhang
    Commented Jan 18, 2020 at 11:25
  • $\begingroup$ 你的这个微分方程是解决什么问题的,感觉好复杂... $\endgroup$ Commented Jan 18, 2020 at 11:32

1 Answer 1


It turns out that the problem isn't about precision, but about the setup of StopIntegration, where I set to be WhenEvent[z[ρ] == 10^-3, "StopIntegration"], WhenEvent[z[ρ] == 1.1` Sqrt[zt^2 - ρc^2], "StopIntegration"], however, if we have a closer look on the solution $z[ρ]$, we would find it's increasing at first, then decreasing until hitting the ρ-axis. However, the upper bound to stop integration 10` Sqrt[zt^2 - ρc^2] is so small that for some specific set of precision, the function just exceeds the bound, which shouldn't be eliminated by it was. And the solution is simply to set the bound to be higher, for example, 1.1` Sqrt[zt^2 - ρc^2]. The bound was originally set to eliminate the solutions that don't hit on the ρ-axis, thus it doesn't matter to set it higher.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.