4
$\begingroup$

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]];
$\endgroup$
  • 2
    $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Jul 2 '15 at 17:41
  • 3
    $\begingroup$ NDSolve solves a DE for the highest-order derivative. At the initial condition, it will divide by zero (i.e. by r^2 for r == 0). $\endgroup$ – Michael E2 Jul 2 '15 at 18:41
  • 1
    $\begingroup$ In these cases, I sometimes try to replace 0 with $MachineEpsilon. $\endgroup$ – J. M. will be back soon Jul 2 '15 at 21:38
5
$\begingroup$

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ Thank you so much! I greatly appreciate the input and assistance. $\endgroup$ – Nathan Melton Jul 2 '15 at 20:53
5
$\begingroup$

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

Mathematica graphics

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}}  *)
$\endgroup$
3
$\begingroup$

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] 

enter image description here

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.

$\endgroup$
  • $\begingroup$ Just for curiosity, how is this kind of problem solved in matlab? $\endgroup$ – xzczd Jul 7 '15 at 8:06
  • $\begingroup$ MATLAB has a number of numerical solvers specifically adapted to boundary value problems. These solvers include options that allow you to solve equations with singular points at $r = 0$ over the interval $[0,b]$. See the documentation here. $\endgroup$ – Michael Seifert Jul 7 '15 at 13:33
0
$\begingroup$

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

Mathematica graphics

$\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.