NDSolve divergence error

Question

I am trying to solve an ODE using NDSolve, but get a divergence error. Could someone please help me wi it? The following is my code.

b = 0.001/0.03;
a = 0.007/0.04;
sg = NDSolve[{r^3 (a (s'[r] + 2 s[r]/r) + (s''[r] + 2 s'[r]/r - 
          2 s[r]/r^2)/
        r + (s'''[r] + 2 s''[r]/r - 4 s'[r]/r^2 + 4 s[r]/r^3) - 
       b (1 + s[r] s[r]) (s'[r] + 2 s[r]/r) - 2 b s[r] s'[r]) == 0, 
   s[0] == 0, s'[0] == 0, s''[0] == 0}, s, {r, 0, 30}]

Just to make sure, your boundary conditions state that the function, its first and second derivatives are all zero for $r=0$. Is that correct and reasonable for your model? Your equation contains a few $1/r$ terms so, at a glance, it may not be well defined at $r=0$. — MarcoB, May 23 '22 at 19:20
FWIW, a solution is the trivial one $s(r)=0$ (identically). It's the one NDSolve would compute. Are you interested in other initial conditions? If so, maybe you should edit the post to change them. — Michael E2, May 23 '22 at 19:43
I only know that s[0]=0. I am not sure about the other conditions but do not know what to put there so that ndsolve can solve it. @MarcoB — sara nj, May 23 '22 at 20:00
I only know that s[0]=0. So do not know what other conditions i can use @MichaelE2 — sara nj, May 23 '22 at 20:01

score 2 · Accepted Answer · answered May 23 '22 at 21:25

The explanation would be really long solution...

{ode = r^3 (a (s'[r] + 2 s[r]/r) + (s''[r] + 2 s'[r]/r - 2 s[r]/r^2)/
        r + (s'''[r] + 2 s''[r]/r - 4 s'[r]/r^2 + 4 s[r]/r^3) - 
       b (1 + s[r] s[r]) (s'[r] + 2 s[r]/r) - 2 b s[r] s'[r]) == 0, 
  s[0] == 0, s'[0] == 0, s''[0] == 0} // ExpandAll

(* Substitute s = u/r^m and solve for a value of m
   that admits a nonzero value for one of u[0], u'[0], u''[0].
   Taking m=1 makes u''[0] free. There are two solutions that
   make u'[0] free (left as an exercise for the reader, if
   there are any).  There are no solutions that make u[0] free
   unless a=b, which is not the case here.                  *)

uode = ode /. s -> Function[r, u[r]/r^m];
defaultICS = {u[r] -> 0, u'[r] -> 0, u''[r] -> 0};
freeterm = u''[r];
uppp = u'''[r] /. First@Solve[uode, u'''[r]];
conds = Expand@CoefficientList[uppp, freeterm] /. defaultICS // 
  Collect[#, r] &
msol = Solve[
     Replace[Series[#, {r, 0, -1}], sd_SeriesData :> sd[[3]]] == 
      0] & /@ conds // Flatten
(*
{0, (-3 + 3 m)/r}
{m -> 1}
*)

(* Take tiny, tiny steps from r=0 until we get to a reasonable
   location for an initial condition.  Keep halving the
   step size until we get a stable solution.
   And hope it's a good solution                            *)

PrintTemporary@Dynamic@{Clock@Infinity, foo};
Block[{a, b, m, f0},
 m = m /. msol;
 b = 0.001/0.03 // Rationalize;
 a = 0.007/0.04 // Rationalize;
 f0 = 1/100;
 pwode = u'''[r] == Piecewise[{{uppp, r != 0}},
    uppp /. freeterm -> f0 /. defaultICS /. r -> 0];
 pwics = {u[0] == u[r], u'[0] == u'[r], u''[0] == u''[r]} /. 
    freeterm -> f0 /. defaultICS;
 ndsol1 = 
  NDSolve[{pwode, pwics}, u, {r, 0, 0.01}, 
   Method -> {"FixedStep", Method -> "ExplicitRungeKutta", 
     "DiscontinuityProcessing" -> None}, 
   StartingStepSize -> 0.00000001, WorkingPrecision -> 16, 
   MaxSteps -> 10^8, StepMonitor :> (foo = r)]
 ]

(* Final solution: *)

PrintTemporary@Dynamic@{Clock@Infinity, foo};
Block[{a, b, m, f0},
 m = m /. msol;
 b = 0.001/0.03 // Rationalize;
 a = 0.007/0.04 // Rationalize;
 f0 = 1/100; (* arbitrarily chose IC for u''[r] *)
 pwode = u'''[r] == Piecewise[{{uppp, r != 0}},
    uppp /. freeterm -> f0 /. defaultICS /. r -> 0];
 pwics = {u[0.01], u'[0.01], 
    u''[0.01]} ==
   ({u[0.01], u'[0.01], u''[0.01]} /. First[ndsol1]);
 ndsol = 
  NDSolve[{pwode, pwics}, u, {r, 0, 30}, StepMonitor :> (foo = r)]
 ]
(* Don't forget that the solution is s[r] = u[r]/r *)
Plot[u[r]/r /. ndsol // Evaluate, {r, 0, 30}]

N.B. The final solution is integrated down to 2.84*10^-16 only, not all the way to 0. No error message was given though, which seems odd.

NDSolve divergence error

1 Answers1

Linked