2
$\begingroup$

I am trying to solve a non linear differential equation with variable parameter.

β[b_, c_] = (4 c - 14)/(12 b - 6 c - 9);
α[b_, c_] = (10 - 4 b)/(12 b - 6 c - 9);
eqn3[b_, c_, A_] = β[b, c] x f[x] - α[b, c] x^2 D[f[x], x] + 2 A^(b - 1) D[Sqrt[x] D[f[x]^b x^(2 + c), x], x];

b1 = 5/3;
c1 = -(1/2);

NDSolve[{eqn3[b1, c1, A] == 0, f[xmin] == 1, f[1] == p}, 
 f[x], {x, xmin, 1}, 
 Method -> {"Shooting", 
   "StartingInitialConditions" -> {f[xmin] == 1, 
     Derivative[1][f][xmin] == -((2 + c1)/(b1 xmin))}}, 
 MaxSteps -> 1000]

with variable xmin, A and p->0. I am getting the error for p<0.01

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

Is there any way to sort the problem?.

$\endgroup$
  • $\begingroup$ Please format your code nicely in code blocks by using either back-ticks or indenting four spaces: for help, click the grey question mark on the right side of the editing toolbar. $\endgroup$ – march Dec 20 '15 at 18:40
  • $\begingroup$ Please give example values for xmin, A, and p. Setting them equal to zero causes problems. $\endgroup$ – march Dec 20 '15 at 18:44
  • 1
    $\begingroup$ Although better formatted, this questions is essentially identical to 102418. Please delete one or the other. $\endgroup$ – bbgodfrey Dec 20 '15 at 19:15
4
$\begingroup$

The error occurs, because f[x] vanishes near x == 0.4976549855466093. This problem can be overcome by experimenting to obtain a better estimate for Derivative[1][f][xmin]. To illustrate, supply as the missing constants,

xmin = .01; p = 10^-3; A = 10;

Then,

s = NDSolve[{eqn3[b1, c1, A] == 0, f[xmin] == 1, f[1] == p}, f[x], {x, xmin, 1}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {f[xmin] == 1, 
    Derivative[1][f][xmin] == -94.81}}];
LogPlot[f[x] /. s, {x, xmin, 1}]

enter image description here

Obtaining smaller p

As noted at the beginning of this answer, the ODE becomes singular at f[x] = 0. Hence, seeking a solution with p = 0 or very near to zero causes the Shooting Method problems. An automated approach for obtaining very small p is as follows. First, obtain a numerical solution parameterized in terms of Derivative[1][f][xmin] == bc.

s = ParametricNDSolveValue[{eqn3[b1, c1, A] == 0, f[xmin] == 1, 
    Derivative[1][f][xmin] == bc}, f, {x, xmin, 1}, {bc}]

and then Plot the value of bc at which f[x] = 0 or the endpoint x = 1 is reached, whichever occurs first.

Plot[(Quiet@s[bc]["Domain"])[[1, 2]], {bc, -100, -90}]

enter image description here

A precise value of the break in the slope of the curve is the desired value of bc, which can be obtain by

bl = -100; bu = -90; inf = 1;
Do[bt = (bl + bu)/2; If[Quiet@s[bt]["Domain"][[1, 2]] < inf, bl = bt, bu = bt], {i, 30}]
N[bu, 12]
(* -94.8105141707 *)
s[bu][1]
(* 4.96698*10^-8 *)

which is a reasonable approximation of zero. For this approach to work, the initial estimates of bl and bu must bracket the desired value of bc.

$\endgroup$
  • $\begingroup$ This is good but can we have a solution for even smaller value of p, for example p<10^(-6) with the range 10^(-4)<xmin<10^(-1). I have always obtained a point of singularity in this which does not seem to be an intrinsic issue. Probably this is a boundary value problem as i can find a solution if p>=0.05. $\endgroup$ – magesh Dec 21 '15 at 4:34
  • 1
    $\begingroup$ @magesh I have added a procedure that yields p < 10^-7. It works equally well for xmin = 10^-3, but a much smaller value of bc is needed, -921.268789973. $\endgroup$ – bbgodfrey Dec 21 '15 at 9:35

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.