# singularity in boundary value problem

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 == p},
f[x], {x, xmin, 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {f[xmin] == 1,
Derivative[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?.

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

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[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 == p}, f[x], {x, xmin, 1},
Method -> {"Shooting", "StartingInitialConditions" -> {f[xmin] == 1,
Derivative[f][xmin] == -94.81}}];
LogPlot[f[x] /. s, {x, xmin, 1}] 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[f][xmin] == bc.

s = ParametricNDSolveValue[{eqn3[b1, c1, A] == 0, f[xmin] == 1,
Derivative[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}] 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]
(* 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.

• 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. – magesh Dec 21 '15 at 4:34
• @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`. – bbgodfrey Dec 21 '15 at 9:35