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
  • 1
    $\begingroup$ At $x=b$, the coefficient of $y''$ vanishes; hence, you have a singularity. $\endgroup$
    – Michael E2
    May 15, 2019 at 0:09

1 Answer 1

3
$\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$
2
  • $\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, 2019 at 14:51
  • $\begingroup$ @Ninpou Read this site first mathematica.stackexchange.com/… $\endgroup$ May 20, 2019 at 14:55

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.