10
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – xzczd Apr 16 at 11:23
  • $\begingroup$ @xzczd ok thanks $\endgroup$ – mattiav27 Apr 16 at 11:38
8
$\begingroup$

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.

enter image description here

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]
| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – Nasser Apr 17 at 0:34
  • $\begingroup$ @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. ) $\endgroup$ – xzczd Apr 17 at 2:08
  • 1
    $\begingroup$ I wonder how Matlab's bvp numerical solver handles this. May be I should try and find out sometime. $\endgroup$ – Nasser Apr 17 at 4:46
  • 1
    $\begingroup$ @xzczd There is a typo in your code it should be D[χ[x], {x, 2}] == χ[x]^(3/2)/Sqrt[x] $\endgroup$ – Alex Trounev Jun 6 at 23:57
  • $\begingroup$ @Alex Oops… Fixed, thx for pointing out. $\endgroup$ – xzczd Jun 7 at 0:54
10
$\begingroup$

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:

enter image description here

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]

enter image description here

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]

enter image description here

But the solution compares favorably to other solutions presented here:

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

enter image description here

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}]

enter image description here

If you are interested in why I reformulated the equations then this section from the documentation is a good starting point.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks I will read it! $\endgroup$ – mattiav27 Apr 15 at 18:38
  • $\begingroup$ 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 $\endgroup$ – mattiav27 Apr 16 at 7:06
  • 3
    $\begingroup$ @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. $\endgroup$ – xzczd Apr 16 at 9:08
  • 2
    $\begingroup$ 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. $\endgroup$ – bbgodfrey Apr 19 at 4:57
  • 1
    $\begingroup$ @user21 There is a typo in your code it should be D[χ[x], {x, 2}] == χ[x]^(3/2)/Sqrt[x] $\endgroup$ – Alex Trounev Jun 6 at 23:58
6
$\begingroup$

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]

Figure 1

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 Figure2

| improve this answer | |
$\endgroup$
5
$\begingroup$

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]

enter image description here

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ It's a bit surprising that WorkingPrecision -> 20 seems to be necessary to make Shooting method work. $\endgroup$ – xzczd Jun 7 at 4:16
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Jun 7 at 12:36
3
$\begingroup$

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.

Mathematica graphics

Note that Maple does not have FEM solver, only FDM.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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