2

I saw some strange behavior as I was solving a differential equation, so I decided to plot the solution in three different conditions.

First, I defined the initial conditions and some constants:

fot = 6.580813053912583`*^-19;
zp = 1000;
lu = 8.418054414588785`*^-33;

Then I defined the three differential equations:

pr1 = ParametricNDSolve[{(1 + x)^5 D[ (r[x])/(1 + x)^4, x] == l0 (r[x] + (1 + x)^3)^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l0}, AccuracyGoal -> 75];
pr2 = ParametricNDSolve[{(1 + x)^5 D[ (r[x])/(1 + x)^4, x] == l0 (r[x] + (1 + x)^3)^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l0}, WorkingPrecision -> 75];
pr3 = ParametricNDSolve[{(1 + x)^5 D[ (r[x])/(1 + x)^4, x] == l0 (r[x] + (1 + x)^3)^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l0}, PrecisionGoal -> 75];

Next I plot the solutions:

plotpr1 = Plot[Evaluate[r[1*10^-22][x] /. pr1], {x, 0, 1}]
plotpr2 = Plot[Evaluate[r[1*10^-22][x] /. pr2], {x, 0, 1}]
plotpr3 = Plot[Evaluate[r[1*10^-22][x] /. pr3], {x, 0, 1}]

the three resulting plots

They give two different plots: in this case AccuracyGoal and WorkingPrecision give the same answer. However in my previous post, I showed that they give different answers (although it was not exactly the same problem).

Question When should I use each one of these options?

Rohit Namjoshi
  • 10,212
  • 6
  • 16
  • 67
No name
  • 585
  • 1
  • 5
  • 17
  • 1
    You may want to read this: https://mathematica.stackexchange.com/q/118249/1871 – xzczd Jun 18 '20 at 02:11
  • I read that post before. I find the solution of my problem with the help of another question that I made [https://mathematica.stackexchange.com/questions/224091/two-different-solutions-for-the-same-contourplot-while-using-accuracygoal-and-wo] . In easy terms, the equation is "stiff". – No name Jun 18 '20 at 02:16

1 Answers1

4

From the documentation of PrecisionGoal: "In most cases, you must set WorkingPrecision to be at least as large as PrecisionGoal". You didn't, so you obtained iffy results, as simple as that.

Note also that, when you request a high working precision, your input should also match or exceed that precision; that's why you get those warnings when you come to evaluating your second ParametricNDSolve upon plotting it:

ParametricNDSolve::precw: The precision of the differential equation (...) is less than WorkingPrecision (75.`).

That's because fot and lu are defined at machine precision. You can fix that by defining them to have higher precision. Here I chose a precision of 80 digits to meet and exceed any extra precision requirements that may be encountered during evaluation at WorkingPrecision -> 75:

fot = 6.580813053912583`80*^-19;
lu = 8.418054414588785`80*^-33;

That fixes the second plot.

For the third one using PrecisionGoal, let's follow the instructions in the docs and set an appropriate WorkingPrecision as well, to a value much higher than the requested PrecisionGoal. I often use the rule of thumb of setting the working precision to at least twice the precision goal to be on the safe side. This comes from the fact that the default setting (PrecisionGoal -> Automatic) "normally yields a precision goal equal to half the setting for WorkingPrecision", again from the Details section of the docs of PrecisionGoal. So that's what I'll do here:

pr3new = ParametricNDSolve[
            {(1 + x)^5 D[(r[x])/(1 + x)^4, x] == l0 (r[x] + (1 + x)^3)^(1/2), 
             r[zp] == fot}, r, {x, 0, 10^8}, {l0},
            PrecisionGoal -> 30, WorkingPrecision -> 60
         ]

plotpr3 = Plot[Evaluate[r[1*10^-22][x] /. pr3new], {x, 0, 1}]

enter image description here

I do not feel sufficiently qualified to thoroughly address your final question, i.e. when you should use each of these options, so I will limit myself to telling you that, in my practice, I use WorkingPrecision much more often than AccuracyGoal and PrecisionGoal.

MarcoB
  • 67,153
  • 18
  • 91
  • 189
  • 3
    Let me help: AccuracyGoal -> a roughly means that error estimates less than 10^-a are acceptable. When the solution is on the order of 10^{-23}, then you need an AccuracyGoal much greater than 23. For instance for AccuracyGoal -> 23 would mean a relative error of around 100% would be acceptable, which it should not be. The first two solutions have accuracy goals of 75, 37.5, resp., each of which is sufficiently large to ensure an accurate solution. The third solution has an accuracy goal of around 8, which is too small. – Michael E2 Jun 16 '20 at 21:13
  • 1
    Also: precw is just a warning. In fact, the system raises the precision for you automatically. The warning is just pointing out an inconsistency in the input, in case it's a user mistake. – Michael E2 Jun 16 '20 at 21:17
  • Thank you, your solutions is what I was looking for. Just for curiosity you saw the other post? – No name Jun 16 '20 at 21:23
  • 1
    Thank you @Michael_E2, your second comment is really useful. – No name Jun 16 '20 at 21:24
  • I have a reason to suspect that this is not the right answer, maybe I'm doing a mistake again, sorry if I'm wrong. Make this plotpr1 = Plot[Evaluate[r[0][x] /. pr1], {x, 0, 1}] and this plotpr3 = Plot[Evaluate[r[0][x] /. pr3new], {x, 0, 1}], once again they don't give the same answer, actually the right answer is the first one, you can prove it solving the diferential equation analitically for the case $l0=0$ – No name Jun 17 '20 at 00:01
  • The solution of the differential equation with $l0=0$ is psolution[x_]:=6.55456*10^-31*(1+x)^4 – No name Jun 17 '20 at 00:13
  • @MarcoB can you look at my last comment please? – No name Jun 17 '20 at 01:21