2
$\begingroup$

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, 
  ImagePadding -> 50, 
  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.

Improve this question
$\endgroup$
5
  • $\begingroup$ you can use "Quiet" to get rid of 1/0. $\endgroup$
    – Gopal Verma
    Commented Mar 4, 2018 at 20:44
  • $\begingroup$ @GopalVerma It didnt work. $\endgroup$
    – Patrick El Pollo
    Commented Mar 5, 2018 at 12:58
  • $\begingroup$ 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. $\endgroup$
    – Gopal Verma
    Commented Mar 5, 2018 at 14:56
  • $\begingroup$ @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. $\endgroup$
    – Patrick El Pollo
    Commented Mar 5, 2018 at 16:13
  • $\begingroup$ One can also avid the singulartiy by adding a small number u[x2]+10^-9. Did you find the the solution for x2>0?. $\endgroup$
    – Gopal Verma
    Commented Mar 5, 2018 at 16:59

1 Answer 1

Reset to default
3
$\begingroup$

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], 
   ImagePadding -> 50, 
   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"}]

Mathematica graphics

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.

Improve this answer
$\endgroup$
15
  • $\begingroup$ Thanks for your answer. I'd like to ask how did you get those values for p0, [p0, 2/ic0^2, E/ic0^2]? $\endgroup$
    – Patrick El Pollo
    Commented Mar 12, 2018 at 11:38
  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Mar 12, 2018 at 13:11
  • 1
    $\begingroup$ @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.) $\endgroup$
    – Michael E2
    Commented Mar 13, 2018 at 22:10
  • 1
    $\begingroup$ @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 $\endgroup$
    – Michael E2
    Commented Apr 6, 2018 at 0:45
  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Apr 6, 2018 at 0:47

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.