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?.

Improve this question
$\endgroup$
3
  • 1
    $\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
    Commented Dec 20, 2015 at 18:40
  • $\begingroup$ Please give example values for xmin, A, and p. Setting them equal to zero causes problems. $\endgroup$
    – march
    Commented Dec 20, 2015 at 18:44
  • 1
    $\begingroup$ Although better formatted, this questions is essentially identical to 102418. Please delete one or the other. $\endgroup$
    – bbgodfrey
    Commented Dec 20, 2015 at 19:15

1 Answer 1

Reset to default
6
$\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.

Improve this answer
$\endgroup$
2
  • $\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
    Commented Dec 21, 2015 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
    Commented Dec 21, 2015 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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