# Solving a simple BVP (with an error)

I am trying to redo the method given in

Solving a simple BVP

which is very nice, to my equation

2 + 2 u'[x2]^2 + u[x2] u''[x2]


where x2=[-1.,1.] and u[-1.]=u[1.]=1/10. I copyed the steps and change the pamareters, obviously. My code in this case is

Manipulate[eq = 2 + 2 u'[x2]^2 + u[x2] u''[x2] == 0;
ic = {u[-1] == ic0, u[1] == ic1};
sol = First@NDSolve[Flatten[{eq, ic}], u[x2], {x2, -1, to}];
Plot[u[x2] /. sol, {x2, -1, to}, Frame -> True, PlotRange -> All,
FrameLabel -> {{u[x2], None}, {x2,
Style[Row[{"solution to ",
2 + 2 Derivative[1][u][x2]^2 +
u[x2] (u^\[Prime]\[Prime])[x2] == 0}], 12]}},
GridLines -> Automatic,
GridLinesStyle -> Directive[LightGray, Thickness[.001]]], {{to, 1,
"to?"}, 0, 1, .01, ImageSize -> Tiny,
Appearance -> "Labeled"}, {{ic0, 1/10, "u(x20)"}, 0, 1, .01,
ImageSize -> Tiny, Appearance -> "Labeled"}, {{ic1, 1/10, "u(x21)"},
0, 1, .01, ImageSize -> Tiny, Appearance -> "Labeled"}]


As a result, I get

Power::infy: Infinite expression 1/0. encountered.

Power::infy: Infinite expression 1/0.^2 encountered.

Power::infy: Infinite expression 1/0. encountered.

General::stop: Further output of Power::infy will be suppressed during this calculation.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

NDSolve::ndnum: Encountered non-numerical value for a derivative at x2 == -1..

and

ReplaceAll::reps: {2+2 (u^[Prime])[-0.999959]^2+u[-0.999959] (u^[Prime][Prime])[-0.999959]==0,u[-1]==1/10,u1==1/10} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

ReplaceAll::reps: {2. +2. (u^[Prime])[-0.999959]^2+u[-0.999959] (u^[Prime][Prime])[-0.999959]==0.,u[-1.]==0.1,u[1.]==0.1} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

ReplaceAll::reps: {2+2 (u^[Prime])[-0.959143]^2+u[-0.959143] (u^[Prime][Prime])[-0.959143]==0,u[-1]==1/10,u1==1/10} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation.

I ask the same question giving my whole code and parameters in

My code

where you can see that u[-1.]=u[1.]=1/100, but thats not the main point.

• you can use "Quiet" to get rid of 1/0. – Gopal Verma Mar 4 '18 at 20:44
• @GopalVerma It didnt work. – Patrick El Pollo Mar 5 '18 at 12:58
• But, using sol = First@Quiet[NDSolve[Flatten[{eq, ic}], u[x2], {x2, -1, to}]] we found that there is no Infinite expression 1/0. encountered. – Gopal Verma Mar 5 '18 at 14:56
• @GopalVerma I saw whats is acctually the function of Quiet in the code, and it look it only hides the message not eliminate the 1/0 singularity. I think the problem is that the range of x2 passes through 0, but I dont know how to modify that because that range is required. – Patrick El Pollo Mar 5 '18 at 16:13
• One can also avid the singulartiy by adding a small number u[x2]+10^-9. Did you find the the solution for x2>0?. – Gopal Verma Mar 5 '18 at 16:59

The automatic Shooting Method does not suffer errors well. It seems to give up when poorly chosen initial conditions lead to an error. In this case, there is a lower limit on the initial condition for u'[-1], below which the solution develops a singularity. It is very close to the actual solution, so the built-in shooting method inevitably runs into a singularity and fails. Thus a manual approach seems to be necessary. We add an extrapolation handler that will cause FindRoot to increase the initial condition when this happens.

Also, one cannot have boundary conditions u[-1] == 0 nor u[1] == 0, since in solving for u''[x2], the equation is divided by u[x2]. So I limited the input range on the sliders for the BCs.

Manipulate[
eq = 2 + 2 u'[x2]^2 + u[x2] u''[x2] == 0;
ic = {u[-1] == ic0, u[1] == ic1};
With[{pen = 3/ic0^2},   (* slightly informed guess *)
psol = ParametricNDSolveValue[
Flatten[{eq, {u[-1] == ic0, u'[-1] == p0}}], u, {x2, -1, 1}, {p0},
"ExtrapolationHandler" -> {p0 - pen &, "WarningMessage" -> False}]
];
Quiet[
usol = psol[p0] /. FindRoot[psol[p0][1] == ic1,
{p0, 2/ic0^2, E/ic0^2}],  (* starting values from inspection of psol *)
ParametricNDSolveValue::ndsz];
Dynamic@Plot[usol[x2], {x2, -1, to}, Frame -> True,
PlotRange -> {{-1, 1}, All}, PlotRangePadding -> Scaled[.02],
FrameLabel -> {{u[x2], None}, {x2,
Style[Row[{"solution to ",
2 + 2 u'[x2]^2 + u[x2] u''[x2] == 0}], 12]}},
GridLines -> Automatic,
GridLinesStyle -> Directive[LightGray, Thickness[.001]]],
{{to, 1, "to?"}, 0, 1, .01, ImageSize -> Tiny, Appearance -> "Labeled"},
{{ic0, 1/10, "u(x20)"}, 0.001, 1, .01, ImageSize -> Tiny, Appearance -> "Labeled"},
{{ic1, 1/10, "u(x21)"}, 0.001, 1, .01, ImageSize -> Tiny, Appearance -> "Labeled"}]


usol'[-1]  (*  p0  found by FindRoot for  ic0 = 0.001  is approx  E/ic0^2 *)
(*  2.78641*10^6  *)


Note how large the derivative is compared to the size of the solution. It might be hard to lower the limit on the slider for ic0, without increasing precision.

• Thanks for your answer. I'd like to ask how did you get those values for p0, [p0, 2/ic0^2, E/ic0^2]? – Patrick El Pollo Mar 12 '18 at 11:38
• For p0, when I got the answer for ic0 = 0.1, I thought p0 is about 100 or so, or about 1/ic0^2 (it was actually 200+, so 2/ic0^2). I solved it for ic0 = 0.01 and again it was about 2/ic0^2. The second point was originally 1/ic0^2, which worked, but I realized that all the solutions I saw were about 2.7 or 2.8 times 1/ic0^2. For the penatly pen, I picked a number that was bigger than that. The value returned will be negative when p0 is too small and integration did not make it to x2 == 1; since ic1 > 0, the secant method will then increase p0 on the next step. – Michael E2 Mar 12 '18 at 13:11
• @resanrom For an initial value for u'[-1] in the IVP, pick successively things like 1, 10, 100, etc. and examine the solution, what it looks like, and how close to the BC at 1 you get. One could do that inside Manipulate, too and adjust the IC with a slider. As for "in principle," I just guessed, from two data points, maybe more. I forget. I didn't have time to do the analysis; maybe someone else can. (Actually, all I could see from the calculations, is that it's close enough a guess for FindRoot to succeed, which I verified by moving the slider.) – Michael E2 Mar 13 '18 at 22:10
• @resanrom I learned of "ExtrapolationHandler" here. Extrapolation is what happens when the input to an interpolating function lies outside the domain. The handler is a function that is applied to such an input. Normally, a warning message is issued, but that can be turned off. A penalty is a common method to enforce a constraint when solving an equation. In this case, extrapolation happens when the integration did not reach the end at x2 == 1; the penalty indicates the – Michael E2 Apr 6 '18 at 0:45
• parameter caused the solution to evaluate to a value <= p0. If the integration completed (reached x2 == 1) the boundary value will be >= p0. The desired value of the parameter will be between those two values, unless in the first case, the solution improbably evaluates exactly to p0. – Michael E2 Apr 6 '18 at 0:47