Two different solutions for the same ContourPlot while using AccuracyGoal and WorkingPrecission

I was making a contour plot, of some differential equation however I get two different solutions If I use "AccuracyGoal" and "Working Precision". The question will appear long but most of the post is code and images.

First I will put the solution that I think is the right one:

Initial conditions

Clear["Global*"]
fot = 6.580813053912583*^-19; zp = 1000; lu = 8.418054414588785*^-33;

Differential equation

pr = ParametricNDSolve[{(1 + x)^5 D[ (r[x])/(1 + x)^4, x] ==
l024 (r[x])^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l024},
AccuracyGoal -> 75];

Plot

ab1 = ContourPlot[((r[l024][x] /. pr)/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2}, PlotLegends -> Automatic]

ab2 = ContourPlot[((r[l024][x] /. pr)/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2}, PlotLegends -> BarLegend[Automatic, LegendMarkerSize -> 180,  LegendFunction -> "Frame", LegendMargins -> 5,  LegendLabel -> "\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)"], Frame -> True,  FrameLabel -> {{"\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)",  ""}, {"\!$$\*SubscriptBox[\(\[Lambda]$$, $$0$$]\)", ""}},  BaseStyle -> {FontWeight -> "Bold", FontSize -> 14}, Contours -> {5}, ContourStyle -> Directive[Thick, Black], ContourShading -> None, PlotRange -> All];

Show[ab1, ab2]
ab2

Then I get: However if I use "WorkingPrecission" (This will take around 2 minutes)

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

ab11 = ContourPlot[((r[l024][x] /. pr2)/(lu))^(1/4), {l024, 0,1.2*10^-22}, {x, 0, 2}, PlotLegends -> Automatic]

ab22 = ContourPlot[((r[l024][x] /. pr2)/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2},PlotLegends ->  BarLegend[Automatic, LegendMarkerSize -> 180, LegendFunction -> "Frame", LegendMargins -> 5,  LegendLabel -> "\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)"], Frame -> True,  FrameLabel -> {{"\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)",  ""}, {"\!$$\*SubscriptBox[\(\[Lambda]$$, $$0$$]\)", ""}},  BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},  Contours -> {5}, ContourStyle -> Directive[Thick, Black],  ContourShading -> None, PlotRange -> All];

Show[ab11, ab22]

ab22 Possible solution

I think that the problem is something related with "PrecisionGoal", but I'm not sure.

• Why does ab22 refer to pr versus pr2? Is it a potential typo? Jun 17 '20 at 4:51
• Yes is a Typo, I will correct it. Jun 17 '20 at 5:22
• Your equation is quite non-linear. I changed the method to StiffnessSwitching and it seemed to help. Jun 17 '20 at 5:24

Your solution in the second case returned complex numbers when evaluated. When I switched the method to StiffnessSwitching, the complex number went away. Note that I changed to ParametricNDSolveValue, because that is what I usually work with.

pr2 = ParametricNDSolveValue[{(1 + x)^5 D[(r[x])/(1 + x)^4, x] ==
l024 (r[x])^(1/2), r[zp] == fot}, r, {x, 0, 10^8}, {l024},
WorkingPrecision -> 75, Method -> "StiffnessSwitching"];
ab11 = ContourPlot[((pr2[l024][x])/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2}, PlotLegends -> Automatic, PlotRange -> All]

ab22 = ContourPlot[((pr2[l024][x])/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2},
PlotLegends ->
BarLegend[Automatic, LegendMarkerSize -> 180,
LegendFunction -> "Frame", LegendMargins -> 5,
LegendLabel -> "\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)"],
Frame -> True,
FrameLabel -> {{"\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)",
""}, {"\!$$\*SubscriptBox[\(\[Lambda]$$, $$0$$]\)", ""}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},
Contours -> {5}, ContourStyle -> Directive[Thick, Black],
ContourShading -> None, PlotRange -> Full];
Show[ab11, ab22]
ab22 In the following, I will do some basic analysis and scaling of the differential equation. I will the subscript $$d$$ to denote the variable/parameter has dimensions. Here is OP's initial equation:

$${\left( {{x_d} + 1} \right)^5}\frac{\partial }{{\partial {x_d}}}\frac{{{r_d}\left( {{x_d}} \right)}}{{{{\left( {{x_d} + 1} \right)}^4}}} = {\lambda _d}\sqrt {{r_d}\left( {{x_d}} \right)} ;{x_d} \geq 0$$

We can use Mathematic to evaluate and simplify the equation to obtain:

$$\frac{{\partial {r_d}\left( {{x_d}} \right)}}{{\partial {x_d}}} = \frac{{{\lambda _d}\sqrt {{r_d}\left( {{x_d}} \right)} + 4{r_d}\left( {{x_d}} \right)}}{{\left( {{x_d} + 1} \right)}}$$

We can define dimensionless variables and parameters like so:

$$x = \frac{{{x_d}}}{{{z_p}}};r = \frac{{{r_d}}}{{{f_{ot}}}};\lambda = \frac{{{\lambda _d}}}{{\sqrt {{f_{ot}}} }}$$

Now, we can create a non-dimensionalized equation like so:

$$\frac{{dr}}{{dx}} = \frac{{4r + \lambda \sqrt r }}{{\frac{1}{{{z_p}}} + x}}$$

We know that at $$r(x=1)=1$$, which implies that the right hand side of the equation is real and positive. Beyond $$x=1$$, $$r$$ is a monotonically increasing function. If we look backwards from $$x=1$$, then $$r$$ should be monotonically decreasing. A singularity occurs at $$x=-\frac{-1}{z_d}$$, but we are always above that point since $$x \geq 0$$. Examining the equation in simplified non-dimensional form, it is difficult to see how $$r$$ could turn complex since the right hand side should be positive.

Here is an example workflow using the non-dimensionalized form. I increased the MaxRecursions in the plot to eliminate the small spikes. Also, I imported the NDSolveUtilities package to look at the timesteps taken by the solver.

Needs["DifferentialEquationsNDSolveUtilities"];
eq = r'[x] == (4 r[x] + λ Sqrt[r[x]])/(1/zp + x);
pr3 = ParametricNDSolveValue[{eq, r == 1},
r, {x, 0, 2}, {λ}, WorkingPrecision -> 75,
Method -> "StiffnessSwitching"];
ab111 = ContourPlot[((pr3[λd/Sqrt[fot]][xd/zp])/(lu/fot))^(1/
4), {λd, 0, 1.2*10^-22}, {xd, 0, 2},
PlotLegends -> Automatic, PlotRange -> All]
ab222 = ContourPlot[((pr3[λd/Sqrt[fot]][xd/zp])/(lu/fot))^(1/
4), {λd, 0, 1.2*10^-22}, {xd, 0, 2}, MaxRecursion -> 4,
PlotLegends ->
BarLegend[Automatic, LegendMarkerSize -> 180,
LegendFunction -> "Frame", LegendMargins -> 5,
LegendLabel -> "\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)"],
Frame -> True,
FrameLabel -> {{"\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)",
""}, {"\!$$\*SubscriptBox[\(λ$$, $$0$$]\)", ""}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},
Contours -> {5}, ContourStyle -> Directive[Thick, Green],
ContourShading -> None, PlotRange -> All];
Show[ab111, ab222]
ab222
StepDataPlot[pr3[(1.2*10^-22)/(2 Sqrt[fot])]] With the StiffnessSwitching method activated, we see a nice smooth transition to the timestep. The following plots show the timestep control for 4 cases that I ran. Setting the AccuracyGoal only looks like a coarse description of when the StiffnessSwitching is turned on. The WorkingPrecision only setting appear to give up on adjusting the timestep when the solution moves away from the initial boundary condition.

Let's check the assumptions of the previous analysis that said r was monotonically increasing and positive by plotting r vs x and r(0) vs $$lambda_d$$ with the following code:

Plot[((pr3[0.6*10^-22/Sqrt[fot]][xd/zp])/(lu/fot))^(1/4), {xd, 0,
2 zp}]
Plot[((pr3[λd/Sqrt[fot]][0/zp])/(lu/fot))^(1/4), {λd,
0, 1.2*10^-22}] The results seem to be consistent with our previous statements.

Finally, let's compare the "ab2" plots of AccuracyGoal Only (red), WorkingPrecision++StiffnessSwitching (green), and WorkingPrecision+AccuracyGoal+StiffnessSwitching (blue).

Show[ab2, ab222, ab2222] The blue curve took the longest, but had the most control and probably the most accurate. One needs to determine if the extra cost it worth it.

• Sorry for my basic question, but How do you know that the complex number is not "real", I mean, you eliminate in an arbitrary way or there is some reason behind it? Jun 17 '20 at 5:24
• I will read the documentation of StiffnessSwitching. If your solution is a real solution you save me! Thank you. Jun 17 '20 at 5:26
• @Cruz If you expect a complex solution, then you won't be able to plot the solutions with ContourPlot unless you transform it into a real. I am not familiar with your problem, but I assumed that you were expecting a real solution. If you look under possible issues of NDSolveValue, it recommended to use Method->"Extrapolation" when using high working precision. I tried that and received a warning that BDF or StiffnessSwitching might be better and that is what I tried. Jun 17 '20 at 5:33
• @Cruz In any case, your numerical solution is dependent on the numerical methods used and there are many possible methods to use. If your equation is meant to describe a physical observable, generally that means it should be real. If it is a purely mathematical question, all bets are off. To make it easier on the numerics, you may want to consider dimensional scaling to keep the values closer to 1. Jun 17 '20 at 5:44
• @Cruz I added an Analysis and Scaling section. The solver time steps for AccuracyGoal only is a coarse approximation to when the StiffnessSwitching is turned on. WorkingPrecision only in this case appears to stop making adjustments to the time step away from the boundary condition. Jun 17 '20 at 20:28

Here is an exact solution:

SetPrecision[
Hold[
fot = 6.580813053912583*^-19;
zp = 1000;
lu = 8.418054414588785`*^-33;],
Infinity] // ReleaseHold

ode = (1 + x)^5 D[(r[x])/(1 + x)^4, x] == l024 (r[x])^(1/2);

rsol = Last@DSolve[{ode, r[zp] == fot}, r, x]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

Solve::ifun: Inverse functions are being used....

(*
{r -> Function[{x}, (Sqrt -
18050427306500947968000 l024 + 2 Sqrt x +
36028797018963968 l024 x + Sqrt x^2 +
18014398509481984 l024 x^2)^2/
5213097220524497879008234839265467517785604096]}
*)

Plots:

ab11 = ContourPlot[((r[x] /. rsol)/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2}, PlotLegends -> Automatic]

ab22 = ContourPlot[((r[x] /. rsol)/(lu))^(1/4), {l024, 0,
1.2*10^-22}, {x, 0, 2},
PlotLegends ->
BarLegend[Automatic, LegendMarkerSize -> 180,
LegendFunction -> "Frame", LegendMargins -> 5,
LegendLabel -> "\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)"],
Frame -> True,
FrameLabel -> {{"\!$$\*SubscriptBox[\(z$$, $$Lss$$]\)",
""}, {"\!$$\*SubscriptBox[\(\[Lambda]$$, $$0$$]\)", ""}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 14},
Contours -> {5}, ContourStyle -> Directive[Thick, RGBColor["#FF2000"]],
ContourShading -> None, PlotRange -> All];

Show[ab11, ab22]

ab22 