# Attempting to use NDSolve to plot harmonic oscillator solutions

NDSolve[{ (-y''[r]/1880) + (470 (0.04077)^2 r^2 - 48 + 1/(1880 r^2)) y[r] == 0,
y[0] = 0, y'[0] = 0}, y, {r, -4, 4}]


I use this but get errors and am unable to get a plot.

Update

Even after I fix the syntax error,

NDSolve[{(-y''[r]/1880) + (470 (0.04077)^2 r^2 - 48 + 1/(1880 r^2)) y[r] == 0,
y[0] == 0, y'[0] == 0}, y, {r, -4, 4}]


I still get Power::infy (1/0.^2), Infinity::indet, and NDSolve::ndnum (non-numerical derivative) errors. Is there a way to integrate and plot this ODE?

• Use == instead of = in your initial conditions. Clear everything first before doing this. Oct 12, 2015 at 11:48
• I've done as you said but been told 'NDSolve called with 2 arguments; 3 or more arguments are expected.' now Oct 12, 2015 at 11:51
• The code in your question has three arguments (list of equations, dependent variable, iterator), doesn't it? If you're getting that error, then that wasn't what you exactly entered. Oct 12, 2015 at 11:54
• @J.M.isback. Edditing this I've added a comma between } and y. Oct 12, 2015 at 11:56
• …and that, folks, is why we try not to touch the code in questions asking for debugging assistance. Oct 12, 2015 at 12:01

[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 (Rationalize@0.04077)^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.)

• The substitution must throw off NDSolve's error estimate. The precision is often only about half of what is normal. Suggestions, anyone? Oct 20, 2015 at 18:12

There is a singularity at r==0 you have to deal with.

{(-y''[r]/1880) + (470 (0.04077)^2 r^2 - 48 + 1/(1880 r^2)) y[r] == 0,
y[0] == 0, y'[0] == 0} /. r -> 0

{False, y[0] == 0, Derivative[1][y][0] == 0}
Power::infy : "\"Infinite expression \[NoBreak]1/0^2\[NoBreak] \
encountered.\""


One way around this maybe to integrate from a r > 0 to r==4.

    ode = NDSolveValue[{Derivative[2][y][
r] - (1880*470*0.04077^2*r^2 - 48 + 1/1880 r^2)*y[r] == 0,
y[10^(-10)] == 0, Derivative[1][y][10^(-10)] == 0},
y[r], {r, -4, 4}];
Plot[ode, {r, -4, 4}]


The solution seems to be 0. everywhere.

ode /. r -> 0
0.
`