# Overflow with NDSolve

Backslide introduced in 10.0, persisting through 12.1.

I am trying to solve the following differential equation and plot the result:

ClearAll["Global*"]
rin = 10^-30;
sol = NDSolve[{Sqrt[x] D[χ[x], {x, 2}] == χ[x]^(3/2), χ[rin] == 1, χ == 0}, χ, {x, rin, 10}, PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 20,Method -> "StiffnessSwitching"]
Chi[x_] := χ[x] /. sol[]
Plot[Chi[x], {x, 0, 10}]


but I get the errors

General::ovfl: Overflow occurred in computation.

NDSolve::ndsz: At x== 5.4454867280576339492985227093479235677620., step size is effectively zero; singularity or stiff system suspected.

EDIT I use version 10.0

• Just tested in v9.0.1, and NDSolve works without any adjustion. So the underlying issue seems to be the backslide mentioned here: mathematica.stackexchange.com/a/130373/1871 Apr 16, 2020 at 11:23
• @xzczd ok thanks Apr 16, 2020 at 11:38

I'd like to extend my comments to an answer. For those in v12 or higher, FiniteElement is a possible choice for this problem, as shown in user21's answer. But, if you're in a version lower than v12 but higher than v9, it becomes a bit more troublesome, because

1. nonlinear FiniteElement isn't implemented yet.

2. Shooting method can't handle the problem well, which is an arguable backslide. As we can see, though ndsz warning is generated, NDSolve manages to find the desired result in v9.

OK, so what to do? Well, to be honest I don't know if the following solution will cause other problem in v10.0, because v10.0 is a quite unstable version, but it does work in v9 and v12.1:

rin = 10^-30;

psol = ParametricNDSolveValue[{D[χ[x], {x, 2}] == (χ[x]^(3/2))/Sqrt[x], χ'[
10] == d, χ == 0}, χ, {x, rin, 10}, d]

drule = FindRoot[psol[d][rin] == 1, {d, 0 (* choose -1/10 if in v9 *)}] // Quiet
(* {d -> -0.0116574} *)

Plot[psol[d /. drule][r] // Evaluate, {r, rin, 10}, PlotRange -> All]


Alternatively, we can turn to finite difference method (FDM). I'll use pdetoae for the generation of finite difference equations:

rin = 0;

eq = D[χ[x], {x, 2}] Sqrt[x] == (χ[x]^(3/2));

bc = {χ[rin] == 1, χ == 0};

points = 25; domain = {rin, 10}; grid = Array[# &, points, domain]; difforder = 2;

del = #[[2 ;; -2]] &;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[χ[x], grid, difforder];

ae = ptoafunc@eq // del;

initialguess[x_] = 0;

solrule = FindRoot[{ae, bc}, Table[{χ[x], initialguess[x]}, {x, grid}]]

sol = ListInterpolation[solrule[[All, -1]], grid]

• I actually tried what you show at the very top as first thing. But it fails on V 12.1. Not just warning, no result is given. Screen shot. !Mathematica graphics Strange it works in V 9 Apr 17, 2020 at 0:34
• @Nasser Yeah, that's why I consider this as a backslide, sadly WRI doesn't think so. (See the discussion in comments in 2nd link for more info. ) Apr 17, 2020 at 2:08
• I wonder how Matlab's bvp numerical solver handles this. May be I should try and find out sometime. Apr 17, 2020 at 4:46
• @xzczd There is a typo in your code it should be D[χ[x], {x, 2}] == χ[x]^(3/2)/Sqrt[x] Jun 6, 2020 at 23:57
• @Alex Oops… Fixed, thx for pointing out. Jun 7, 2020 at 0:54

Update:

Alex noted a typo in my answer. Sorry about this. You can still solve this equation with the nonlinear FEM solver, thought it's not a as straight forward as for the miss typed equation.

Using

NDSolveValue[{D[\[Chi][x], {x, 2}] == \[Chi][x]^(3/2)/Sqrt[x], \[Chi][
rin] == 1, \[Chi] == 0}, \[Chi], {x, 0, 10},
Method -> "FiniteElement"]


gives an error message: This is a fairly general failure message. The reason it is so general is that for the code it's impossible to say why it failed. One cause can the that there is a transition from the real to the complex plain. An easy way to try this is to either give an complex valued initial seed or, alternatively, to add a 0. I complex component to the equation:

sol = NDSolveValue[{D[\[Chi][x], {x, 2}] == \[Chi][x]^(3/2)/
Sqrt[x], \[Chi][rin] == 1, \[Chi] == 0}, \[Chi], {x, 0, 10},
Method -> "FiniteElement", InitialSeeding -> \[Chi][x] == 0. I];
Plot[Re[sol[x]], {x, 0, 10}, PlotRange -> All] Note the Re in the plot. Now, there is also a small complex component in the solution:

Plot[Im[sol[x]], {x, 0, 10}, PlotRange -> All] But the solution compares favorably to other solutions presented here:

Plot[Evaluate[psol[d /. drule][r] - Re[sol[r]]], {r, rin, 10},
PlotRange -> All] rin = 0;
sol = NDSolve[{D[\[Chi][x], {x, 2}] == (\[Chi][x]^3/2)/
Sqrt[x] , \[Chi][rin] == 1, \[Chi] == 0}, \[Chi], {x, rin,
10}, Method -> "FiniteElement"];
Chi[x_] := Evaluate[\[Chi][x] /. sol[]];
Plot[Chi[x], {x, 0, 10}] If you are interested in why I reformulated the equations then this section from the documentation is a good starting point.

• Thanks I will read it! Apr 15, 2020 at 18:38
• I have just tried to run your code, but I get the error NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve I use version 10.0 Apr 16, 2020 at 7:06
• @mattiav27 Nonlinear FEM is supported by NDSolve since v12. For earlier versions, playing with Shooting method should help, FDM is also a choice. For both methods we already have many related posts, just search in this site. Apr 16, 2020 at 9:08
• Although the plot looks reasonable (compare 219895), the first derivative is far from accurate near the origin without Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> .01}} or finer. Using higher-order finite elements also might be helpful. Apr 19, 2020 at 4:57
• @user21 There is a typo in your code it should be D[χ[x], {x, 2}] == χ[x]^(3/2)/Sqrt[x] Jun 6, 2020 at 23:58

There is also wavelet method for BVP. It is an example with Haar wavelets. It takes 0.36 s to solve this problem with 64 colocation points:

ClearAll["Global*"]
L = 10; A = 0; B = 1; J = 5; M =
2^J; dx = (B - A)/(2 M);
h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}];
p1[x_, n_] := (1/n!)*(x - A)^n;
h[x_, k_, m_] :=
Piecewise[{{1,
Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {-1,
Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0];
p[x_, k_, m_, n_] :=
Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!,
Inequality[k/m, LessEqual, x,
Less, (1 + 2*k)/(2*m)]}, {((-(k/m) + x)^n -
2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, (1 + 2*k)/(2*m) <=
x <= (1 + k)/
m}, {((-(k/m) + x)^n + (-((1 + k)/m) + x)^n -
2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0];
xl = Table[A + l dx, {l, 0, 2 M}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
f2[x_] :=
Sum[af[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0 h1[x];
f1[x_] :=
Sum[af[i, j] p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0 p1[x, 1] + f10;
f0[x_] :=
Sum[af[i, j] p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0 p1[x, 2] + f10 x + f00;
bc1 = {f0 == 1};
bc2 = {f0 == 0};
var = Flatten[Table[{af[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
varM = Join[{a0, f10, f00}, var];
eqq[x_] := Sqrt[x] f2[x]/L^(3/2) - f0[x]^(3/2);
eq = Flatten[Table[{eqq[x] == 0}, {x, xcol}]];
eqM = Join[eq, bc1, bc2];
sol = FindRoot[eqM, Table[{varM[[i]], 0.1}, {i, Length[varM]}],
MaxIterations -> 1000]; lst =
Table[{L x, Evaluate[f0[x] /. sol]}, {x, 0, 1, .01}];

ListLinePlot[lst, PlotRange -> All]


Now we can compare with solution by xzczd psol[d /. drule][r], here we show difference f0[x]-psol[d /. drule][x] And as we can see, combination of ParametricNDSolveValue[] and FindRoot[] is still good Thinking for a moment suggests χ == 0 is a numerically problematic boundary condition for two reasons: First, the χ[x]^(3/2) term, which means that χ[x] should not go negative or the integration will run into branch cut difficulties. Second, the solution is concave up so that if it has a turning point (a positive minimum), it will increase; in fact, since the power 3/2 is greater than 1, most solutions, if not all, will go off to infinity in finite time. These in turn suggests starting the shooting method at the x == 10 boundary condition, instead of the default lesser value x == rin chosen automatically by NDSolve.

rin = 10^-30;
sol = NDSolve[
{Sqrt[x] D[χ[x], {x, 2}] == χ[x]^(3/2), χ[rin] == 1, χ == 0},
χ, {x, rin, 10}, PrecisionGoal -> 8, AccuracyGoal -> 8,
WorkingPrecision -> 20,
Method -> {"Shooting",
"StartingInitialConditions" -> {χ == 0, χ' == -1/100}}];

Chi[x_] := χ[x] /. sol[]
Plot[Chi[x], {x, 0, 10}, PlotRange -> All] • It's a bit surprising that WorkingPrecision -> 20 seems to be necessary to make Shooting method work. Jun 7, 2020 at 4:16
• @xzczd For reasons I can't fathom, the shooting method is varying x near x == 0 instead of integrating to x == rin. I'd think any programmer would understand that stepping outside the specified domain is inviting trouble. You can fix it by using Sqrt[x + $MinMachineNumber] instead of Sqrt[x], which won't make a difference if x > 10^-292, since x +$MinMachineNumber rounds to exact x in machine precision. It's a pretty ugly fix though. Jun 7, 2020 at 12:36

This is just to confirm @user21 result using Maple. I used bvp[middefer] (For BVP problem like this one) in Maple to get solution. I am not sure what this maps to in NDSolve now if any.

restart;
ode := diff(X(x),x\$2)=X(x)^(3/2)/sqrt(x);
bc  := X(0)=1,X(10)=0;
sol := dsolve([ode,bc],X(x),numeric,method=bvp[middefer],abserr=0.001,maxmesh=8192*4);
plots:-odeplot(sol, x=0..10);
`

Similar (but not exact) looking solution is generated. Note that Maple does not have FEM solver, only FDM.