Message Boards Message Boards

Error with NDSolveValue: step size is effectively zero

Amin bkh
Posted 2 years ago

Hi all, I am trying to solve a second order equation, in Mathematica as follows:

ysol = NDSolveValue[{y''[x] == -4 Pi (0.0021) Exp[-y[x]] + (0.0021) 4 Pi Exp[y[x]] , y[6.94444] == 0, 
   y'[1.38889] ==-4.125}, y, {x, 1.38889, 6.94444}]

and I get the following error message:

At x == 1.3888888888888888`, step size is effectively zero; \ singularity or stiff system suspected.

Could you please help me to realize where the problem is and how I can solve it?

Thank you so much.
POSTED BY: Amin bkh
18 Replies
Sort By:
Amin bkh
Amin bkh
Posted 2 years ago

Thank you so much for all your help.
POSTED BY: Amin bkh
Amin bkh
Amin bkh
Posted 2 years ago

Sorry Mariusz, this way does not work

ysol = NDSolve[{y''[x] == -4 Pi (a) Exp[-y[x]] + (a) 4 Pi Exp[y[x]], 
   y'[c] == 10^-6, y'[b] == d}, y, {x, b, c}, 
  Method -> {"BoundaryValues" -> {"Shooting", 
      "StartingInitialConditions" -> {y[c] == 0, y'[c] == -1/2}}}, 
  WorkingPrecision -> 50]
POSTED BY: Amin bkh

Try:

$Version
(*"13.0.0 for Microsoft Windows (64-bit) (December 3, 2021)"*)
ClearAll["`*"]; Remove["`*"];

a = Rationalize[0.0021, 0];
b = Rationalize[1.38889, 0];
c = Rationalize[6.94444, 0];
d = Rationalize[-4.125, 0];

ysol = NDSolve[{y''[x] == -4 *Pi*a*Exp[-y[x]] + 4*Pi*a*Exp[y[x]], 
   y'[c] == 0, y'[b] == d}, y, {x, b, c}, 
  Method -> {"BoundaryValues" -> {"Shooting", 
      "StartingInitialConditions" -> {y[c] == Rationalize[1.653, 0], 
        y'[c] == 0}}}, WorkingPrecision -> 50]; Plot[{y[x] /. ysol, 
  y'[x] /. ysol}, {x, b, c}, PlotLabels -> {"y[x]", "y'[x]"}]
 {y'[x] /. ysol /. x -> c, y'[x] /. ysol /. x -> b}

Or: 

ysol1 = NDSolve[{y''[x] == -4 *Pi*a*Exp[-y[x]] + 4*Pi*a*Exp[y[x]], 
   y'[c] == 0, y'[b] == d}, y, {x, b, c}, 
  Method -> {"BoundaryValues" -> {"Shooting", 
      "StartingInitialConditions" -> {y[b] == Rationalize[5.792, 0], 
        y'[b] == d}}}, WorkingPrecision -> 50]; Plot[{y[x] /. ysol1, 
  y'[x] /. ysol1}, {x, b, c}, PlotLabels -> {"y[x]", "y'[x]"}]
{y'[x] /. ysol1 /. x -> c, y'[x] /. ysol1 /. x -> b}

Or: 

ysol3 = NDSolve[{y''[x] == (-4*Pi*a*Exp[-y[x]] + a*4*Pi*Exp[y[x]]) + 
     NeumannValue[0, x == c] + NeumannValue[d, x == b]}, y, {x, b, c},
   Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2}, 
    "IntegrationOrder" -> 5, 
    "MeshOptions" -> {MaxCellMeasure -> 0.001}}, 
  InitialSeeding -> {y[x] == x}]
Plot[{y[x] /. ysol3, y'[x] /. ysol3}, {x, b, c}, 
 PlotLabels -> {"y[x]", "y'[x]"}]
{y'[x] /. ysol3 /. x -> c, y'[x] /. ysol3 /. x -> b}

Regards.
POSTED BY: Mariusz Iwaniuk
Amin bkh
Amin bkh
Posted 2 years ago

Is it possible I come up with a method that Mathematica by itself put an approximated point for starting? for y'[c] == -1/2? Once again thank you so much.
POSTED BY: Amin bkh

Probably not, for more information read this.

Regards.
POSTED BY: Mariusz Iwaniuk
Amin bkh
Amin bkh
Posted 2 years ago

Thanks a lot for all your help.
POSTED BY: Amin bkh
Amin bkh
Amin bkh
Posted 2 years ago

And one more thing, is it possible to make y'[c]==0?
POSTED BY: Amin bkh

With y'[c]=0 no way, but with: y'[c] == 10^-6 works.
POSTED BY: Mariusz Iwaniuk
Amin bkh
Amin bkh
Posted 2 years ago

@Mariusz Thank you so much for your help and answer, I deeply appreciate it.
POSTED BY: Amin bkh

With FEM:

$Version
(*"13.0.0 for Microsoft Windows (64-bit) (December 3, 2021)"*)

ysol2 = NDSolve[{y''[x] == (-4 *Pi*a*Exp[-y[x]] + a *4 *Pi *Exp[y[x]]) + 
 NeumannValue[d, x == b], DirichletCondition[y[x] == 0, x == c]}, y, {x, b, c}, 
Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2}, 
"IntegrationOrder" -> 5, 
"MeshOptions" -> {MaxCellMeasure -> 0.001}}, 
InitialSeeding -> {y[x] == x}]
Plot[{y[x] /. ysol2, y'[x] /. ysol2}, {x, b, c}, 
 PlotLabels -> {"y[x]", "y'[x]"}]

 y[x] /. ysol2 /. x -> c
 (*{0}*)
 y'[x] /. ysol2 /. x -> b
 (*-4.125*)
POSTED BY: Mariusz Iwaniuk
Amin bkh
Amin bkh
Posted 2 years ago

Thank you so much for all your kind help and attention.
POSTED BY: Amin bkh

One way:

a = Rationalize[0.0021, 0]; 
b = Rationalize[1.38889, 0]; 
c = Rationalize[6.94444, 0]; 
d = Rationalize[-4.125, 0];

ysol = NDSolve[{y''[x] == -4 Pi (a) Exp[-y[x]] + (a) 4 Pi Exp[y[x]], 
   y[c] == 0, y'[b] == d}, y, {x, b, c}, 
  Method -> {"BoundaryValues" -> {"Shooting", 
      "StartingInitialConditions" -> {y[c] == 0, y'[c] == -1/2}}}, 
  WorkingPrecision -> 50]

Plot[{y[x] /. ysol, y'[x] /. ysol}, {x, b, c}, PlotLabels -> {"y[x]", "y'[x]"}]

We can check boundary conditions:

y[x] /. ysol /. x -> c
(*{0.*10^-50}*)
y'[x] /. ysol /. x -> b
(*{-4.12500000000000000000071556836358857851310492}*)
POSTED BY: Mariusz Iwaniuk
Amin bkh
Amin bkh
Posted 2 years ago

Unfortunately, this is not the case for the equation I have written. Numerical result does not show singularity whatsoever. Also for Sinh function between these two point there are not singularities.
POSTED BY: Amin bkh

You should expect explosions. Take for example the equation y''[t] == y[t]^3: an explicit solution is Sqrt[2]/( t - Sqrt[2] ), which is has a singularity.
POSTED BY: Gianluca Gorni

You get an analytical solution this way:

DSolveValue[{y''[x] == Rationalize[(0.0021)]*8 \[Pi] Sinh[y[x]]}, 
 y[x], x]

but it is not very easy to deal with, because of the branch cuts.
POSTED BY: Gianluca Gorni
Amin bkh
Amin bkh
Posted 2 years ago

Yes, also I can obtain it on paper, so the problem is I see no reason I encounter singularity :(
POSTED BY: Amin bkh

I am no expert in boundary value problems, but those two exponentials in the differential equation makes it very prone to explosions, as soon as the value of y gets larger than a few units.
POSTED BY: Gianluca Gorni
Amin bkh
Amin bkh
Posted 2 years ago

It has a closed from analytical equation, but I wonder why Mathematica is not able to solve it.
POSTED BY: Amin bkh
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract