3
$\begingroup$

Here is the code I want to numerically integrate my nonlinear ODE"

sol = NDSolve[{y''[x] + y[x] == -1/(1 + y[x])^2, y[0] == 1, y'[0] == 1}, y, {x, 0, 5}]
Plot[y[x] /. sol, {x, 0, 3}]

NDSolve surely fails since it passes through the singular at the point x=2.532686648521246.

I tried another method from here to bypassing a singularity:

odey = R''[y] + R[y] + 1/(1 + R[y])^2;
BCm = {R[0] == 1, R'[0] == 1};
projODE = {odey /. R -> (p[#]/q[#] &) // Together // Numerator, 
 D[p[y]^2 + q[y]^2, {y, 2}]} == 0 // Thread;

projICS = BCm /. {R[y0_] == r0_ :> {p[y0], q[y0]} == {Numerator[r0], Denominator[r0]}, 
R'[y0_] == 
  rp0_ :> ({D[p[y]/q[y], y] == rp0, D[p[y]^2 + q[y]^2, y] == 0} /. 
   y -> y0)};

 rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, 
 WorkingPrecision -> 50, MaxSteps -> 20000};
 projSOL = NDSolve[{projODE, projICS}, {p, q}, {y, 0, 5}, rules, 
 Method -> {"EquationSimplification" -> {Automatic,
  "SimplifySystem" -> True}}]

 Plot[{p[y]/q[y] /. First@projSOL}, {y, 0, 4}, PlotLegends -> {p/q}]

but dosen't work.

My question is that, is there a convenient numerical continuation method to bypass the singularity so that NDSolve does not fail from range: {x,0,5} ?

$\endgroup$
1
  • $\begingroup$ I think a good place to start is to try to copy/paste the code Michael posted in the thread you refer to. I tried and ran into a mess: missing variables and lots of errors. So I think first get his code working using the problem in that thread, then apply it to yours. For example, variable "rat" is not defined anywhere in that thread but I think it can be 10^-10 for example. $\endgroup$
    – josh
    Dec 16 '21 at 23:11
3
$\begingroup$

You need a method that is adapted to stiff DE. E.g. "LinearlyImplicitEuler"

Here is your adapted example:

sol = NDSolve[{y''[x] + y[x] == -1/(1 + y[x])^2, y[0] == 1, 
   y'[0] == 1}, y, {x, 0, 5}, Method -> "LinearlyImplicitEuler"]
Plot[y[x] /. sol, {x, 0, 3}, PlotRange -> All]

enter image description here

$\endgroup$
3
  • $\begingroup$ Clever solution, fits well to the phaseplot in my answer! $\endgroup$ Dec 17 '21 at 11:24
  • $\begingroup$ The result of Plot[y[x] /. sol, {x, 0, 5}, PlotRange -> All] looks very strange. $\endgroup$
    – user64494
    Dec 17 '21 at 11:56
  • $\begingroup$ Could you also look at what happens if you try Method -> "StiffnessSwitching"? $\endgroup$ Dec 17 '21 at 15:19
2
$\begingroup$

To long for a comment:

Perhaps StreamPlot gives an idea of the possible solution?

StreamPlot[{yp, -y - 1/(1 + y)^2}, {y, -2, 1}, {yp, -3, 3},FrameLabel -> {y[x], y'[x]}]  

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.