1
$\begingroup$

I am interested in solving differential equations in the form of power series. Let's say we have following equation:

$$f^{\prime \prime} (\rho) + \left( \frac{2 e^{-k \rho}}{\rho} - \varepsilon \right) f(\rho) = 0$$

f''(rho) + (2 Exp(-k*rho)/rho - epsilon) f(rho) = 0   (1)

We can write f(rho) in the form of power series:

f(rho) = Sum[c[n] rho^n, {n,0,max}]
Exp[-k*rho] = Sum[(-k*rho)^n/n!, {n,0,max}]

We can collect terms next to 1/rho, 1, rho, rho^2, ... to get equations for coefficients c[0], c[1], ...:

expression = Collect[Sum[c[n + 2] (n + 1) (n + 2) rho^n, {n, 0, max}] + 
  2/rho Sum[(-k rho)^n/n!, {n, 0, max}] Sum[
    c[n] rho^n, {n, 0, max}] - epsilon Sum[
    c[n] rho^n, {n, 0, max}], rho]

After specifying c[0] = 0 and c[1] = 1 (c[0] = 0 is required by equations, c[1] = 1 is arbitrary) we get set of equations (from fact that expression = 0, so every coefficient in power expansion have to be zero):

s3 = Solve[(k^2 c[0] - 2 k c[1] - \[Epsilon] c[1] + 2 c[2] + 6 c[3] == 
     0) /. {c[0] -> 0, c[1] -> 1, c[2] -> -1}, c[3]];

s4 = Solve[(-(1/3) k^3 c[0] + k^2 c[1] - 2 k c[2] - \[Epsilon] c[2] + 
       2 c[3] + 12 c[4] == 0) /. 
    Flatten@{c[0] -> 0, c[1] -> 1, c[2] -> -1, s3}, c[4]];

etc.

Is there a way to automatize this procedure, so that the input would be an equation, or some expression in the form of power series and desired order N and the output would be solution for coefficients c[0], ..., c[N]?

I thought about something like this:

rho*expression/.rho->0

which gives coefficient next to 1/rho (and kills every other order) but this gives only the lowest order and again I would have to set up the equations manually.

Thanks.

$\endgroup$
2
$\begingroup$

Probably something like this can be done:

n = 5; (* highest-order term *)
f = Sum[C[k] ρ^k, {k, 1, n}] + O[ρ]^(n + 1); (* C[0] == 0 already imposed *)

Solve[Thread[CoefficientList[D[f, {ρ, 2}] + (2 Exp[-k ρ]/ρ - ε) f, ρ] == 0],
      Array[C, n]] // Simplify
   {{C[2] -> -C[1], C[3] -> 1/6 (2 + 2 k + ε) C[1], 
     C[4] -> -(1/36) (2 + 8 k + 3 k^2 + 4 ε) C[1], 
     C[5] -> 1/360 (2 + 33 k^2 + 6 k^3 + 10 ε + 3 ε^2 + 4 k (5 + 3 ε)) C[1]}}

along with a Solve::svars warning since C[1] was not specified. Increase n as seen fit.

$\endgroup$
2
$\begingroup$

This is your equation:

eq= f''[x]+(2 Exp[-k*x]/x-e) f[x]==0

For a given value of max, you are looking for a series solution sol:

max=5; sol=Function[x, Sum[c[i] x^i, {i,1,max}] + O[x]^max]

Substitute this solution in your equation and use the function SolveAlways for finding the coefficients:

SolveAlways[eq /. f->sol, x]

(* {{c[4]->-(1/36) (2+4 e+8 k+3 k^2) c[1],c[3]->1/6 (2+e+2 k) c[1],c[2]->-c[1]}} *)
$\endgroup$
  • $\begingroup$ Now I don't know which solution should I accept, both of them do the job, but this one is perhaps a bit more clear... If I try this for max at least 7 I get two sets of solutions, one with "correct" (non-zero) coefficients, the second with all coefficients being zero. Why is that so? $\endgroup$ – user16320 Apr 11 '16 at 22:33
  • $\begingroup$ It seems that for max at least 7, the function SolveAlways does not give correct results, A rule for the variable e is returned as well. Therfore, I would accept the other, correct, answer. $\endgroup$ – Fred Simons Apr 12 '16 at 7:34

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.