Topic 5 Series Solutions

5.1 Power series solutions

Recall that a point \(x_0\) is a regular point of the equation \(y''+p(x)y'+q(x)y=0\), where \(p\) and \(q\) are rational functions, if \(x_0\) is not a discontinuity of both \(p(x)\) and \(q(x)\). At a regular point \(x_0\), we can expect a power series solution \(y(x)=\sum\limits_{n=0}^\infty a_nx^n\). The coefficients are be determined recurrence relations that can be obtained by the method of undetermined coefficients. In maple, we can find a power series solution by the method described above, or using the command powsolve supported by the package powerseries, or using dsolve with the option series. The powsolve returns a procedure which can be expressed explicitly but truncated using the command tpsform.

Example 5.1 Find the general solution in power series of the equation \[y''-xy'-y=0.\]

Solution. We first defined the differential equation.

ode41:=diff(y(x), x, x)-x*diff(y(x), x)-y(x)=0;

Method 1: Follow the manual procedure.

Suppose \(y(x)=\sum\limits_{n=0}^\infty a_nx^n\). In Maple, we define \(y\) using the command sum.

y:=sum(a[n]*x^n, n=0..infinity);

Now plug \(y\) in to the equation and simplify.

pode41 := combine(subs[eval](y(x) = y, ode41));

The output is \[\color{blue}{\sum\limits_{n=0}^\infty\left(-n a_{n} x^{n -2}+a_{n} x^{n -2} n^{2}-n a_{n} x^{n}-a_{n} x^{n}\right) = 0}.\]

For polynomials, Maple can combine like terms. Unfortunately, for series, Maple couldn’t do the shifting or combine like terms for us. Since we are interested in the coefficient of the general term in \(x^n\), we can extract the terms that contains \(x^n\). The trick it to substitute n=0..infinity by n=k..k+2. Because, only when \(n=k\), \(k+2\), there are terms in \(x^k\).

kterms:=simplify(subs((n=0..infinity)=(n=k..k+2), pode41));

Now we can extract the coefficient of \(x^k\) using the command coeff. To keep the equation form, we use map to apply the command coeff to \(x^k\) in both sides of the equation.

coeffk:=map(coeff, kterms, x^k);

The output is \[\color{blue}{\left(1+k \right) \left(k a_{k +2}-a_{k}+2 a_{k +2}\right) = 0}.\]

We can then get the recurrence relations by isolate \(a_{k+2}\) and replace \(k\) by \(n-2\).

recurrence_relations:=subs[eval](k=n-2, isolate(coeffk, a[k+2]));

The output is \[\color{blue}{a_{n} = \frac{a_{n -2}}{n}}.\]

To get an explicit formula for \(a_n\), we use the Maple command rsolve

coeffn:=rsolve(recurrence_relation, a[n]);

We can get a polynomial approximation up to \(x^{10}\) for the solution \(y\).

y:=a[0] + a[1]*x + sum(eval(coeffn, n = k)*x^k, k = 2 .. 10);

Method 2: Using the command powsolve.

Load the package powseries

with(powseries):

Solve the equation using powsolve(equation). Note that the output is a procedure that that produces in each step the coefficient of the degree _k term.

ypowsol:=powsolve(ode41);

The coefficient of the \(n\)-th term is given by

an:= eval(ypowsol(_k), _k=n);

To get a truncated power series, we use the Maple command tpsform. For example, the truncated power series for ypowsol up to 10th order can be obtained by the following command.

tpsform(ypowsol, x, 10);

Method 3: Using `dsolve(equation, y(x), type=‘series’).

Using this method, we get a power series solution directly.

dsolve({ode41, y(0) = a[0], D(y)(0) = a[1]}, y(x), 'series');

Running this command gives the solution \[y\left(x \right) = a_{0}+a_{1} x +\frac{1}{2} a_{0} x^{2}+\frac{1}{3} a_{1} x^{3}+\frac{1}{8} a_{0} x^{4}+\frac{1}{15} a_{1} x^{5}+\mathrm{O}\left(x^{6}\right)\]

Each method has advantages and disadvantages.

The first method works in general, but not efficient.
The powsolve command can only work with polynomial coefficient equations and the power series solution is always at 0.
The dsolve command cannot produce recurrence relations.

In all three methods, initial conditions can be added.

In the first method, initial conditions can be imposed with rsolve({recurrence_relation, a[0]=a0, a[1]=a1}).
In the second and the third method, initial conditions can be imposed by {equation, y(0)=a0, y'(0)=a1}.

Exercise 5.1 Find the general solution in power series form of the equation \[y''+ xy'+y=0.\]

5.2 Euler Equations

For differential equations with regular singular point, the first method and the the third method still work. If the equation is an Euler equation, that is, \(ax^2y''+bxy'+c=0\). The solution is easier to find even by hands. Since once can always use dsolve to solve a differential equation. We will focus on how to solve Euler equations manually. Recall the following theorem.

Theorem 5.1 Suppose the solutions of the indicial equation \[ar(r-1)+br+c=0\] are \(r_1\) and \(r_2\). Then the general solution of the Euler equation \[ax^2y''+bxy'+cy=0\] on \((0,\infty)\) is

\[y= c_1x^{r_1}+c_2x^{r_2}\] if \(r_1\) and \(r_2\) are distinct real numbers;
\[y= x^{r}(c_1+c_2\ln x)\] if \(r_1=r_2=r\);
\[y=x^{\lambda}\left[c_1\cos\left(\omega\ln x\right)+ c_2\sin\left(\omega\ln x \right)\right]\] if \(r_1,r_2=\lambda\pm \mathrm{i}\omega\) with \(\omega>0\).

From the theorem, we see that the key to find the general solution is to solve the indicial equation \(ar(r-1)+br+c=0\). This equation can be obtained substituting \(y''\) by \(r(r-1)\), \(y'\) by \(r\), and \(x\) and \(y\) by \(1\). It can also be obtained by the Maple command indicialeq(ode, indepvar, regsp, depvar) supported by DETools.

Example 5.2 Find the general solution of the Euler equations \[x^2y''-xy'+5y=0\]

Solution. Define the equation

ode42:=x^2*diff(y(x), x$2)-x*diff(y(x), x)+5*y(x)=0;

The indicial equation can be obtained using one of the following command

indeq1:=subs[eval]({diff(y(x), x, x)=r*(r-1), diff(y(x), x)=r, y(x)=1, x=1}, ode42);
indeq2:=DETools[indicialeq](ode42, x, 0, y(x));

Solve the indicial equation, we will use the equation indeq1,

rcomplex:=solve(indeq1, r);

The solutions are \[\color{blue}{1+2I, 1-2I}\]

Construct the solution using the theorem. First find the real and imaginary parts.

alpha:=Re(rcomplex[1]);
beta:=Im(rcomplex[2]);

So the general solution is \[y=e^{alpha x}(c_1\cos(\beta\ln x)+c[2]\sin(\beta\ln x)).\]

You can check that using the Maple command odetest.

yg := x -> x^alpha*(c[1]*cos(beta*ln(x)) + c[2]*sin(beta*ln(x))):
odetest(y(x)=yg(x), ode42);

The output \(\color{blue}{0}\) means that yg is a solution of the equation ode42.

Exercise 5.2 Find the general solution of the Euler equation \[x^2y''-2y'-4y=0.\]

Exercise 5.3 Find the general solution of the Euler equation \[x^2y''-3y'+5y=0.\]

Exercise 5.4 Find the general solution of the Euler equation \[x^2y''+3y'+y=0.\]