[Update notice: I had left in the previous code initial conditions for NDSolve
from working through the OP's problem and forgot to generalize them. They are now fixed.]
Introduction
The method of Frobenius obtains a series solution to an second-order, linear ODE $y''+P\,y'+Q\,y=0$ at a regular singular point $x_0$ in the form
$$y(x) = (x-x_0)^m u(x) = (x-x_0)^m \sum_{k=0}^\infty a_k (x-x_0)^k\,,$$
where $u(x)$ is analytic in a neighborhood $x_0$ and $m$ is a root of the indicial equation.
If we make the substitution $y(x) \mapsto x^m u(x)$ in the ODE, we should expect to get a differential equation for $u$ that is nonsingular in a neighborhood of $x_0$.
One should not expect to be able to solve every initial-value problem at a singularity. The factor $(x-x_0)^m$ may force $y(x)$ or $y'(x)$ to be zero or to go to infinity at $x = x_0$. If there are two independent solutions and one of them is infinite at $x_0$, then it is impossible to obtain every initial condition from a linear combination of the two. In the OP's problem, there are two independent Frobenius-series solutions; one of them it turns out satisfies $y(0)=y'(0)=0$.
OP's Example
The function (currently) frobeniusNDSolve
returns two independent solutions from the singular point r == 0
of the OP's ODE over the interval {-4, 4}
.
opODE = (-y''[r]/1880) + (470 ([email protected])^2 r^2 - 48 + 1/(1880 r^2)) y[r];
sol = frobeniusNDSolve[opODE, y, {r, 0, {-4, 4}}]
Two independent solutions are returned. They are highly oscillatory and one is much, much larger than the other.
Plot[(y[r] /. sol)/{1000000, 1} // Evaluate, {r, -0.05, 0.05}]
Any linear combination of them is also a solution. Here is a solution to the IVP y[1] == 1, y'[1] == -1
.
ivpsol = {a, b} /. First@NSolve[{1, -1} == ({a, b}.({y[1], y'[1]} /. sol))]
(* {-0.0236628, 14203.7} *)
Plot[ivpsol.(y[r] /. sol) // Evaluate, {r, 0.999995, 1.000005}, AspectRatio -> Automatic]
Plot[ivpsol.(y[r] /. sol) // Evaluate, {r, 0., 1.2}]
Code dump
This is a rather minimal implementation. It handles the easy case in which the indicial equation of a second-order, linear differential equation has distinct roots whose difference is not an integer. Otherwise, it will return unevaluated. Also missing is the full NDSolve
syntax for the return value.
frobeniusNDSolve[ode, y, {x, x0, x2}, ndsolveopts]
frobeniusNDSolve[ode, y, {x, x0, {x1, x2}}, ndsolveopts]
returns two independent solutions to the linear second-order ode
, valid over x0 <= x <= x2
or x1 <= x <= x2
if possible. The point x0
should be a regular singular point with distinct roots of the indicial equation not differing by an integer.
ClearAll[frobeniusODE];
frobeniusODE[ode_, y_ -> u_, {x_, x0_}, m_: \[FormalM]] :=
y''[x] - (y''[x] /. First@Solve[ode == 0, y''[x]]) /.
y -> ((# - x0)^m u[#] &) // Collect[#, (x - x0)^m] &
ClearAll[linearODEQ];
linearODEQ[ode_, y_, x_] := Length@ coefficientsODE[ode, y, x] === 2;
ClearAll[orderODE];
orderODE[ode_] := Module[{d},
d = Cases[ode, Derivative[n : __][_][_] :> {n}, Infinity];
If[Max[Length /@ d] === 1,
Max@ d, (* ODE *)
$Failed] (* not an ODE *)
]
ClearAll[coefficientsODE];
mem : coefficientsODE[ode_, y_, x_] :=
mem = With[{order = orderODE[ode]},
With[{yp = Derivative[order][y][x]},
CoefficientArrays[yp /. Solve[ode == 0, yp],
Table[Derivative[n][y][x], {n, 0, order - 1}]
]
]];
ClearAll[indicialCoefficients];
indicialCoefficients[ode_, y_, {x_, x0_}] :=
With[{res = -Limit[#, x -> x0] & /@
({(x - x0)^2, x - x0} First@ Last@ coefficientsODE[ode, y, x])},
res /; VectorQ[res, NumericQ]
];
ClearAll[indicialRoots];
indicialRoots[ode_, y_, {x_, x0_}, m_: \[FormalM]] /;
linearODEQ[ode, y, x] && orderODE[ode] == 2 :=
Module[{c},
c = indicialCoefficients[ode, y, {x, x0}];
Solve[m (m - 1) + c.{1, m} == 0, m] /; FreeQ[c, indicialCoefficients]
]
(* Produces two independent solutions *)
Clear[frobeniusNDSolve];
Options[frobeniusNDSolve] = Options[NDSolve];
abs[x_] := Piecewise[{{x, x >= 0}, {-x, x < 0}}];
frobeniusNDSolve[ode_, y_, {x_, x0_, x2_?NumericQ}, opts : OptionsPattern[]] :=
frobeniusNDSolve[ode, y, {x, x0, {x0, x2}}, opts];
frobeniusNDSolve[ode_, y_, {x_, x0_, {x1_, x2_}},
opts : OptionsPattern[]] :=
Module[{roots, u, m, ode2, ode3, icp, sol},
roots = indicialRoots[ode /. Equal -> Subtract, y, {x, x0}];
(ode2 = Collect[u''[x] - (u''[x] /.
First@ Solve[
frobeniusODE[ode /. Equal -> Subtract, y -> u, {x, x0}, #] == 0, u''[x]]),
{u[x], u'[x], u''[x]}, Simplify];
ode3 = ode2 /.
With[{up = Coefficient[ode2, u'[x]],
uq = Coefficient[ode2, u[x]]},
{up u'[x] -> Piecewise[{{up, x != x0}}] u'[x],
uq u[x] -> Piecewise[{{uq, x != x0}}] u[x]}
];
icp = Solve[SeriesCoefficient[ode2, {x, x0, -1}] == 0 /. u[x0] -> 1,
u'[x0]];
If[icp === {{}},
icp = 0,
If[MatchQ[icp, {{_Rule}}],
icp = u'[x0] /. First@icp,
icp = $Failed]];
sol = NDSolve[
{ode3 == 0, {u[x0], u'[x0]} == {1, icp}},
u, {x, x1, x2}, opts];
{y -> Function @@ {x, abs[x - x0]^# u[x] /. First[sol]}}) & /@
Flatten@ Values[roots] /;
FreeQ[roots, indicialRoots] && ! IntegerQ@ First@ Differences[Flatten@ Values[roots]]
]
The code also works if the indicial roots are conjugate complex numbers. To get independent solutions, one should take the real and imaginary parts. (This could be made automatic, but as yet, I haven't programmed it.)
==
instead of=
in your initial conditions. Clear everything first before doing this. $\endgroup$}
andy
. $\endgroup$