Re: Power series solution to differential equations
- To: mathgroup at smc.vnet.net
- Subject: [mg41606] Re: Power series solution to differential equations
- From: "Peltio" <peltio at twilight.zone>
- Date: Wed, 28 May 2003 04:57:35 -0400 (EDT)
- References: <basnpa$kpe$1@smc.vnet.net>
- Reply-to: "Peltio" <peltioNOSP at Miname.com.invalid>
- Sender: owner-wri-mathgroup at wolfram.com
"Dr. Wolfgang Hintze" wrote
>Given a differential equation of the form
>diffeq = a[x] u''[x] + b[x] u'[x] + f[x, u[x]] == 0
>where ' means d/dx we assume that u[x] has a power series expansion
>about x0 of the form (t = x-x0)
>u[t] = Sum[ c[k] t^(k+z) , {k, 0, Infinity }]
>We have to determine z and the coefficients c[k].
>
>Question: what is the best way to tackle this problem in Mathematica?
>Any hint is greatly appreciated.
You could let Mathematica solve the equation with the power method 'by
hand'.
The package Summa.m (it somewhere in the mathsource, I think it has been
classified under "typesetting" : | ) can do this by using an object Summa
that looks like an ordinary Sum (the newer version - which I'll upload to
the mathsource as soon as I find a little bit of time to take care of it -
can switch on and off the evaluation of sums and is totally transparent to
the user). Anyway, here's how it works with one of the equations shown in
this thread
<<Summa.m
eq = (x + 1) y''[x] + y'[x] + x y[x] == 0;
We look for a solution in the form of power series. The function and its
derivatives are defined in terms of a Summa (version 2 of the package allows
for the use of Sum)
ys = Summa[ a[n] x^n, {n, 0, Infinity}];
y1s = DSum[ys, x];
y2s = DSum[y1s, x];
The equation becomes
eqs = eq /. {y[x] -> ys, y'[x] -> y1s, y''[x] -> y2s}
(1 + x) Sum[( n-1) n x^(n-2) a[n], {n, 0, Infinity}] +
Sum[n x^(n-1)a[n], {n, 0,Infinity}] +
x Sum[x^n a[n], {n, 0, Infinity}] == 0
Let's expand the products and bring the factors inside the sums, simplifying
the result:
eqs = ExpandSum[BringIn[eqs]];
eqs = AlignExponent[eqs, x, n -> n];
res = BreakSum[eqs, 1]//SimplifySum
a[1] + 2 a[2] +
Sum[x^n (a[n-1]+(n+1)((n+1)a[n+1]+(n+2)a[n+2])), {n,1,Infinity}] == 0
Now all we have to do is to apply the principle of identity of polynomials.
With a slight abuse of notation we suppress the Sum[_, it] structure
eqRicors = GetElement[res]
a[1] + 2 a[2] +
x^n (a[n-1]+(n+1)((n+1)a[n+1]+(n+2)a[n+2])) == 0
So that we can easily equate the coefficients of x^n and x^0 to zero,
solving for the unknown coefficients:
eq0=GetCoefficient[eqRicors[[1]], x,0]==0;
{ci} = Solve[eq0, a[2]]
{{ a[2] -> -a[1]/2 }}
eqn=GetCoefficient[eqRicors[[1]], x^n]==0;
{cn} = Simplify[Solve[eqn, a[n+2]] /. n->n-2]
{{ a[n] -> -(a[n-3]+(n-1)^2 a[n-1])/(n(n-1)) }}
A particular solution with a0 = 1 and a1 = 0 can be found by fiddling with
the recursion relation we have just found.
c[0] = 1;
c[1] = 0;
c[2] = -c[1]/2;
c[n_ /; n > 2] := c[n] = -((c[n - 1] (n - 1)^2 + c[n - 3])/((n - 1) n))
Table[c[k], {k, 0, 9}]
{1,0,0,-1/6,1/8,-1/10,4/45,-19/240,191/2688,191/2688,-2921/45360}
All in all, the method is not as elegant as the others already proposed, but
if you need a symbolic solution you can have mathematica do what you would
have done with paper and pencil. It could be used, I suppose, to tackle that
problem of yours with the symbolic value z that need to be determined.
(I said "I think" : )) )
Hope this helps,
Peltio
Invalid address in reply-to. Crafty demunging required to mail me.