3
$\begingroup$

I am trying to solve a non-linear second order boundary value problem in a finite interval. The differential equation is $$y''-\frac{a}{b}y-\frac{u_n}{b}y^3-\frac{ge_0}{b}x=0,$$ with $x\in[0,L]$, along with the boundary conditions, $$y'(0)-\frac{o}{b}y(0)=0,\qquad y'(L)+\frac{l}{b}y(L)=\frac{ge_0L}{b}.$$ In the code attached below, with the chosen parameters, the solution is regular and shows up until $a=1.2$. On using $a=1.3$ and beyond, I get the error "At x == 0.99505895226764`, step size is effectively zero; singularity or stiff system suspected."

I am changing the input right after $z$ in Plot[Evaluate...

Working:

ClearAll;
Clear[Derivative]
fneeonnd[z_?NumericQ, a_?NumericQ, b_?NumericQ, o_?NumericQ, 
  l_?NumericQ, g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ, 
  un_?NumericQ, eo_?NumericQ] := (
  ss = NDSolve[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 
      0, y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, 
    y, {x, 0, L}
    , Method -> "StiffnessSwitching", WorkingPrecision -> 100]; 
  y[z] /. ss)
Plot[Evaluate[
  fneeonnd[z, 1.2, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1]], {z, 0, 
  1}, PlotStyle -> {Thickness[0.007]}, BaseStyle -> {FontSize -> 20}, 
 Frame -> {True}, FrameLabel -> {"z", \[Eta]}, 
 FrameStyle -> {{Black, Directive[Thick]}, {Black, 
    Directive[Thick]}, {Black, Directive[Thick]}, {Black, 
    Directive[Thick]}}]

NOT Working:

ClearAll;
Clear[Derivative]
fneeonnd[z_?NumericQ, a_?NumericQ, b_?NumericQ, o_?NumericQ, 
  l_?NumericQ, g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ, 
  un_?NumericQ, eo_?NumericQ] := (
  ss = NDSolve[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 
      0, y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, 
    y, {x, 0, L}
    , Method -> "StiffnessSwitching", WorkingPrecision -> 100]; 
  y[z] /. ss)
Plot[Evaluate[
  fneeonnd[z, 1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1]], {z, 0, 
  1}, PlotStyle -> {Thickness[0.007]}, BaseStyle -> {FontSize -> 20}, 
 Frame -> {True}, FrameLabel -> {"z", \[Eta]}, 
 FrameStyle -> {{Black, Directive[Thick]}, {Black, 
    Directive[Thick]}, {Black, Directive[Thick]}, {Black, 
    Directive[Thick]}}]

a=1.25 a=1 a=0.5

$\endgroup$
3
  • 1
    $\begingroup$ Look at your ODE. The main terms are: y'' = y^3. What does this mean? y'' has to to with the curvature, that is, if y>0, the larger y, the more it is upwards bent and the faster it grows. The same for y<0. Therefor you can solve only for values where y stays small. $\endgroup$ Nov 6, 2020 at 21:23
  • 1
    $\begingroup$ As far as I could see, y was actually getting smaller as I was increasing the input a. One can obtain the plots for a=0.5,0.75...1.2 to see this. $\endgroup$ Nov 6, 2020 at 21:26
  • $\begingroup$ As Daniel Huber says, curvature explodes at certain x. It shows to depend strongly on parameter b. Chose b == .006 and you can go up with a to a ==2 . $\endgroup$
    – Akku14
    Nov 7, 2020 at 4:39

2 Answers 2

3
$\begingroup$

Often, Method -> "Shooting" with a good guess for the "StartingInitialConditions" is necessary to solve ODE boundary value problems. Here, we use a similar do-it-yourself shooting approach. As an example, use the parameters given for the "not working" case in the question:

{a, b, o, l, g, L, c, u, un, eo} =
    {1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1};

sp = ParametricNDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0, 
    y'[0] - (o/b)*y[0] == 0, y[0] == y0}, {y[x], y'[L] + (l/b)*y[L] - (g*eo*L/b)}, 
    {x, 0, L}, {y0}];

FindRoot[Last[sp[y00]], {y00, -.02}, Evaluated -> False]
First[sp[y00 /. %]];
Plot[%, {x, 0, L}, ImageSize -> Large, AxesLabel -> {x, y}, 
    LabelStyle -> {15, Bold, Black}, PlotRange -> All]
(* {y00 -> -0.041088} *)

enter image description here

Of course, a reasonable guess is needed for y00 in FindRoot, here somewhere in the range {-.01, -.08}. A more ambitious case, a = 3, also works, yielding

(* {y00 -> -0.0163729} *)

enter image description here

Addendum: Simplified Approach

It turns out that solving the linearized ODE yields a sufficiently good seeding for the standard Method -> "Shooting" to produce solutions of the nonlinear ODE for all values of a. The example below provides results for interger values of a between 0 and 10.

Clear[a]; 
Table[
    tem = DSolveValue[{y''[x] - (a/b)*y[x] - (g*eo/b)*x == 0, 
    y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, y[0], {x, 0, L}];
    NDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0, 
    y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, y[x], {x, 0, L}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {y[0] == tem, 
    y'[0] == (o/b) y[0]}}], {a, 0, 10}];
Plot[Evaluate@%, {x, 0, L}, ImageSize -> Large, AxesLabel -> {x, y}, 
    LabelStyle -> {15, Bold, Black}, PlotRange -> All, 
    PlotLegends -> Placed[Range[0, 10], {.3, .6}]]

enter image description here

$\endgroup$
2
$\begingroup$

One can use the FEM method or use the FEM method to find good starting initial conditions. I do the latter in fneeonnd but save the FEM solution in solFEM. The check at the end shows that the shooting method did a better job of satisfying the boundary conditions.

fneeonnd[z_, a_?NumericQ, b_?NumericQ, o_?NumericQ, l_?NumericQ, 
  g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ, un_?NumericQ, 
  eo_?NumericQ] := (
  Clear[x, y];
  solFEM = 
   NDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x ==
       NeumannValue[(o/b)*y[x], x == 0] + 
       NeumannValue[-(g*eo*L/b) + (l/b)*y[x], x == L]}, 
    y, {x} \[Element] Line[{{0}, {L}}]];
  NDSolveValue[
    {y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0,
     bcs = {y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}}, 
   y, {x, 0, L}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {y[0] == solFEM[0], 
       y'[0] == solFEM'[0]}}(*,WorkingPrecision\[Rule]100*)])

solFN = fneeonnd[z, 1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1];
Plot[solFN[z], {z, 0, 1}, PlotStyle -> {Thickness[0.007]}, 
 BaseStyle -> {FontSize -> 20}, Frame -> {True}, 
 FrameLabel -> {"z", \[Eta]}, 
 FrameStyle -> {{Black, Directive[Thick]}, {Black, 
    Directive[Thick]}, {Black, Directive[Thick]}, {Black, 
    Directive[Thick]}}]

enter image description here

Plot[{solFN[z], solFEM[z]}, {z, 0, 1}]

enter image description here

bcs /. Equal -> Subtract /. {{y -> solFN}, {y -> solFEM}}
(*  {{0., 2.55954*10^-6},        <-- Shooting
     {-0.00153343, -0.0726405}}  <-- FEM      *)
$\endgroup$

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.