Infinite Expression Error NDSolve

Question

I'm trying to figure out why I am consistently getting an infinite expression error when I'm trying to solve this DE using NDSolve. Any help will be appreciated.

k = NDSolve[{r^2 k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) - 
   2 (k[r] + 1) (2 - k[r]) == 0, k[0] == -1, k[1] == 0}, 
   k, {r, 0, 1}, 
   Method -> {"Shooting", "StartingInitialConditions" -> {k[1] == 0}}, 
   MaxSteps -> \[Infinity]];

Answer 1

One approach is to begin the integration very slightly beyond the singularity at the origin. That, plus correction of the shooting initial conditions, provides an accurate solution.

rs = .001;
s = NDSolveValue[{r^2 k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) -
    2 (k[r] + 1) (2 - k[r]) == 0, k[rs] == -1, k[1] == 0}, k, {r, rs, 1}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {k[rs] == -1, k'[rs] == .001}}];
Plot[s[t], {t, rs, 1}, AxesLabel -> {t, k}, LabelStyle -> {Black, Bold, 12}]

It is not difficult to show analytically that k has the form -1 + a r^2 near the origin (a a constant), and the numerical solution has this character too, as desired.

Answer 2

There are two problems with the ODE at r == 0: (1) To solve for k''[r] you have to divide by r^2, that is, divide by zero, (2) and k'[0] == 0 for all solutions, which frustrates the shooting method:

ode = r^2 k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) - 2 (k[r] + 1) (2 - k[r]) == 0;

Solve[ode, k''[r]]
% /. k[r] -> -1 // Apart  (* initial condition k == -1 at r == 0 *)
(*
  {{k''[r] ->
    -2 (2 - k[r] - 2 k[r]^2 + k[r]^3 - r k'[r] + r k[r] k'[r] - r^2 k'[r]^2) / 
     (r^2 (-1 + k[r]))}}
  {{k''[r] -> -(2 k'[r])/r - k'[r]^2}}
*)

If k'[r] approaches a finite, non-zero value at r == 0, then k''[r] has a pole, which would imply k'[r] goes to infinity; therefore k'[r] must approach zero. As I understand how the shooting method is implemented, NDSolve would use FindRoot to find an initial value k'[0] such that the solution satisfies the boundary condition k[1] == 0. But that's going to be hard to do, if the ODE is simultaneously constraining k'[0] to be zero. The coefficient that needs to be varied is k''[0]. (An alternative approach, such as bbgodfrey's, is to move the starting point far enough away from r == 0 so that the numerics allow NDSolve to hit the target condition k[1] == 0.)

To allow NDSolve to shoot with k''[0], we need a third-order equation. We can differentiate the ODE to get one, and since we know we want k[0] == -1 and k'[0] == 0, the boundary condition k[1] == 0 gives us the third condition we need for a solution. That leaves the problem of dividing by zero.

Using the shooting method with the third-order equation solves the numerics problem, and that allows us to bring the starting value for r close to zero.

Block[{r0 = 1*^-6},
 sol1 = NDSolve[
   {D[ode, r],
    k[r0] == -1, k'[r0] == 3.512`16 r0, k[1] == 0},
   k, {r, r0, 1}]]

In fact we can help NDSolve out a little by using Taylor's theorem to approximate k[r0], k'[r0]. Inspecting the solution above, we can see that k''[r] is about 3.512 near r == 0.

Block[{r0 = 1*^-6, q0 = 3.512},
 sol2 = NDSolve[
  {D[ode, r],
   k[r0] == -1 + q0 r0^2/2, k'[r0] == q0 r0, k[1] == 0},
  k, {r, r0, 1}]]

Plot[k[r] /. First@sol2, {r, 1*^-6, 1}]

If we check the boundary conditions, they look pretty good:

Block[{r0 = 1*^-6},
 {{k[r0] + 1, k'[r0], k''[r0]}, {k[1]}} /. First@sol2
 ]
(*  {{1.75604*10^-12, 3.512*10^-6, 3.51366}, {-4.08313*10^-7}}  *)

Answer 3

The division-by-zero error you're getting occurs when $r \to 0$, and it's probably because your ODE has a regular singular point there. You can see this because if you remove the r^2 term from in front of the ODE, Mathematica can solve it nearly instantaneously:

NDSolve[{ k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) - 
     2 (k[r] + 1) (2 - k[r]) == 0, k[0] == -1, k[1] == 0}, k, {r, 0, 1}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {k[1] == 0}}, 
  MaxSteps -> \[Infinity]]]
Plot[k[r] /. sol, {r, 0, 1}, PlotRange -> All] 

My guess is that Mathematica is dividing through to put your ODE in some kind of canonical form where the highest-order derivative term has a coefficient of 1; and that's why you're getting these division-by-zero errors.

Unfortunately, I don't know of any particularly good techniques to get around this in Mathematica; whenever I've dealt with equations like this (arising from cylindrical or spherical coordinates), I've ended up using MATLAB instead. I'll be interested to see if anyone has good ideas on how to solve such equations.

Answer 4

At least for this specific problem, it seems not to be a bad idea to make a step backward i.e. turn to finite difference method. pdetoae is a general purpose function for discretizing differential equations to algebraic equations, its definition can be found in this post.

points = 25;
{rL, rR} = {0, 1};
{kL, kR} = {-1, 0};
grid = Array[# &, points, {rL, rR}];
xdifforder = 4;
ptoa = pdetoae[k[r], grid, xdifforder];
eqn = ptoa[r^2 k''[r] + 2 r k'[r] + 2 r^2 (k'[r]^2)/(1 - k[r]) - 
     2 (k[r] + 1) (2 - k[r]) == 0];
bc = {k[rL] == kL, k[rR] == kR};
var = k /@ grid;
init = Array[# &, points, {kL, kR}];
sol = ListInterpolation[
   var /. FindRoot[{#[[2 ;; -2]] &@eqn, bc}, {var, init}\[Transpose]], grid];

(* Alternatively *)
(*
lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] := 
 FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

sol = ListInterpolation[
   var /. lSSolve[Subtract @@@ Flatten@{eqn, bc}, {var, init}\[Transpose], 
      MaxIterations -> 500][[2]], grid];
 *)

Plot[sol[r], {r, rL, rR}]