Overflow with NDSolve

Question

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, χ[10] == 0}, χ, {x, rin, 10}, PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 20,Method -> "StiffnessSwitching"]
Chi[x_] := χ[x] /. sol[[1]]
Plot[Chi[x], {x, 0, 10}]

but I get the errors

General::ovfl: Overflow occurred in computation.

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

EDIT I use version 10.0

Answer 1

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, χ[10] == 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, χ[10] == 0};

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

del = #[[2 ;; -2]] &;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)    
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]

Answer 2

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][10] == 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][10] == 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]

Old answer:

How about:

rin = 0;
sol = NDSolve[{D[\[Chi][x], {x, 2}] == (\[Chi][x]^3/2)/
      Sqrt[x] , \[Chi][rin] == 1, \[Chi][10] == 0}, \[Chi], {x, rin, 
    10}, Method -> "FiniteElement"];
Chi[x_] := Evaluate[\[Chi][x] /. sol[[1]]];
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.

Answer 3

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[0] == 1};
bc2 = {f0[1] == 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

Answer 4

Thinking for a moment suggests χ[10] == 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, χ[10] == 0},
   χ, {x, rin, 10}, PrecisionGoal -> 8, AccuracyGoal -> 8, 
   WorkingPrecision -> 20, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {χ[10] == 0, χ'[10] == -1/100}}];

Chi[x_] := χ[x] /. sol[[1]]
Plot[Chi[x], {x, 0, 10}, PlotRange -> All]

Answer 5

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.