1
$\begingroup$

I have a 2nd Order ODE problem I am trying to solve numerically in Wolfram:
$$ (k^2-1+\frac{2x}{x^2+c^2})\frac{d^2y}{dx^2}+\frac{c^2-x^2}{c^2+x^2}\frac{dy}{dx}+x\frac{x^2-3c^2}{(x^2+c^2)^3}y=0 $$ where $c=0.7$,$k=\sqrt{1-\frac{2b}{b^2+c^2}}$ for $b=50$. My initial conditions are $y(b)=1$ $\frac{dy}{dx}(b)=0$.

I get the following error:

NDSolve::ndsz: At x == 50.`, step size is effectively zero; singularity or stiff system suspected.

My code is here:

sola = 
  NDSolve[
    {(k^2 - 1 + 2*x/(x^2 + c^2))*
       D[y[x], x, x] + ((c^2 - x^2)/(x^2 + c^2)^2)*
        D[y[x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*y[x] == 0, 
     y[b] == 1, y'[b] == 0}, 
    y, {x, 0, 100}]

Any idea of what to do here?

$\endgroup$
  • 1
    $\begingroup$ At $x=b$, the coefficient of $y''$ vanishes; hence, you have a singularity. $\endgroup$ – Michael E2 May 15 at 0:09
2
$\begingroup$

We can use the method of the false transient to solve this problem. We transform the equation to the form of the heat equation and then integrate it on both sides of the singular point using FEM:

Needs["NDSolve`FEM`"]

c = .7; b = 50; k = Sqrt[1 - 2*b/(b^2 + c^2)]; x0 = 100; x1 = 50;
sola = NDSolveValue[{-D[y[t, x], t] + 
     D[y[t, x], x, 
      x] + (((c^2 - x^2)/(x^2 + c^2)^2)*
         D[y[t, x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*
         y[t, x])/(k^2 - 1 + 2*x/(x^2 + c^2)) == 
    NeumannValue[0, x == x1], y[0, x] == 0, 
   DirichletCondition[y[t, x] == 1, x == x1]}, 
  y, {t, 0, 10}, {x, x1, x0}, 
  Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2}, 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];
solb = NDSolveValue[{-D[y[t, x], t] + 
     D[y[t, x], x, 
      x] + (((c^2 - x^2)/(x^2 + c^2)^2)*
         D[y[t, x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*
         y[t, x])/(k^2 - 1 + 2*x/(x^2 + c^2)) == 
    NeumannValue[0, x == x1], y[0, x] == 0, 
   DirichletCondition[y[t, x] == 1, x == x1]}, 
  y, {t, 0, 10}, {x, 0, x1}, 
  Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2}, 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];

The solution quickly converges on t and it looks like this

{Plot[Evaluate[Table[sola[t, x], {t, 1, 10, 1}]], {x, x1, x0}, 
  PlotRange -> All], 
 Plot[Evaluate[Table[sola[t, x], {t, 1, 10, 1}]], {x, x1, x1 + 1}, 
  PlotRange -> All]}

{Plot[Evaluate[Table[solb[t, x], {t, 1, 10, 1}]], {x, 0, x1}, 
  PlotRange -> All], 
 Plot[Evaluate[Table[solb[t, x], {t, 1, 10, 1}]], {x, x1 - 1, x1}, 
  PlotRange -> All]}

fig1

Combining solutions we find

Plot[If[x < 50, solb[10, x], sola[10, x]], {x, 0, 100}, 
 PlotRange -> All]

fig2

$\endgroup$
  • $\begingroup$ Could please recommend any material (Books, pdf...) where I may read about this method applied to Ordinary differential equations like this example? $\endgroup$ – Ninpou May 20 at 14:51
  • $\begingroup$ @Ninpou Read this site first mathematica.stackexchange.com/… $\endgroup$ – Alex Trounev May 20 at 14:55

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.