Finding a solution to this equation using Frobenius method

  1. Hi, I have this equation to solve.

    y'' + (1/x)y' + [(x^2) + k + (m^2 / x^2)]y = 0

    now, I've tried to solve this using frobenius method but cannot formulate a solution.

    I have that

    a_(n+4) = [-ka_(n+2) - a_(n)] / [n^2 +/- 2inm]

    is my recurrence relation, but now I'm stuck and don't know how to proceed, any help greatly appreciated.

    thanks.
     
  2. jcsd
  3. Using the standard convention of letting [itex]a_0=1[/itex], I get:

    [tex]a_1=0[/tex]
    [tex]a_2=-\frac{h a_0}{(2+c)(1+c)+1+m^2}[/tex]

    with the remaining odd a_n=0 and even a_n equal to:

    [tex]a_n=-\frac{a_{n-2}+a_{n-4}}{(n+c)(n+c+1)+1+m^2}[/tex]

    with a similar relation for the other root expressed in b_n and so I can write the solution as:

    [tex]y(x)=K_1 \sum_{n=0}^{\infty} a_n x^{n+mi}+K_2 \sum_{n=0}^{\infty} b_n x^{n-mi},\quad x>0[/tex]

    and I believe because of the conjugates in the solution, for real initial conditions [itex]y(x_0)=y_0[/itex], and [itex]y'(x_0)=y_1[/itex], the imaginary component is annihilated leaving the desired real solution.

    If you into it, here's some Mathematica code to check the solution. I believe it's correct but not as close as I would have expected. May be some convergence issues though.

    Code (Text):

    nmax = 75;
    x0 = 0.1;
    y0 = 0;
    y1 = 1.;
    h = 5;
    m = 4;

    mysol = NDSolve[{x^2*Derivative[2][y][x] + x*Derivative[1][y][x] + (x^4 + h*x^2 + m^2)*y[x] == 0, y[x0] == y0,
         Derivative[1][y][x0] == y1}, y, {x, x0, 2}];

    sol[x_] := Evaluate[y[x] /. mysol];
    p1 = Plot[y[x] /. mysol, {x, x0, 2}];

    c1 = I*m;
    Subscript[a, 0] = 1;
    Subscript[a, 1] = 0;
    Subscript[a, 2] = ((-h)*Subscript[a, 0])/((2 + c1)*(1 + c1) + 1 + m^2);
    Subscript[a, 3] = 0;

    Table[Subscript[a, n] = -(Subscript[a, n - 4] + h*Subscript[a, n - 2])/((n + c1)*(n + c1 - 1) + 1 + m^2),
       {n, 4, nmax}];

    f1[x_] := Sum[Subscript[a, n]*x^(n + c1), {n, 0, nmax}]
    f1d[x_] = D[f1[x], x];

    c2 = (-I)*m;
    Subscript[b, 0] = 1;
    Subscript[b, 1] = 0;
    Subscript[b, 2] = ((-h)*Subscript[b, 0])/((2 + c2)*(1 + c2) + 1 + m^2);
    Subscript[b, 3] = 0;

    Table[Subscript[b, n] = -((Subscript[b, n - 4] + h*Subscript[b, n - 2])/((n + c2)*(n + c2 - 1) + 1 + m^2)),
       {n, 4, nmax}];

    f2[x_] := Sum[Subscript[b, n]*x^(n + c2), {n, 0, nmax}];
    f2d[x_] = D[f2[x], x];

    myks = First[NSolve[{y0 == k1*f1[x0] + k2*f2[x0], y1 == k1*f1d[x0] + k2*f2d[x0]}, {k1, k2}]];

    myy[x_] := k1*f1[x] + k2*f2[x] /. myks;

    myd1[x_] = D[myy[x], x];
    myd2[x_] = D[myy[x], {x, 2}];

    N[x^2*myd2[x] + x*myd1[x] + (x^4 + x^2*h + m^2)*myy[x]] /. x -> 0.334

    p2 = Plot[myy[x], {x, x0, 2}, PlotStyle -> Red]
    Show[{p1, p2}]
     
     
    Last edited: Sep 19, 2011
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook