470.html

Series Solutions near an Ordinary Point

Example: Solve y''+y=0 .

Solution: We look for a solution of the form .

> y:=sum(a[k]*x^k,k=0..infinity);y1:=sum(k*a[k]*x^(k-1),k=1..infinity);y2:=sum(k*(k-1)*a[k]*x^(k-2),k=2..infinity);

> y2+y;

Using copy and pasn=left> Using copy and paste change the index of summation. In the first sum replace k-2 with k, and don't forget to correctly adjust the summation limits

> sum((k+2)*(k+1)*a[k+2]*x^k,k = 0 .. infinity)+sum(a[k]*x^k,k = 0 .. infinity);

Write this as a single sum by collecting the powers of x.

> sum(((k+2)*(k+1)*a[k+2]+a[k])*x^k,k = 0 .. infinity);

Hence the recurrence relation is , ... .

We solve this recurence relation as follows:

> restart;

> l:=[seq(a[i],i=0..10)];

> a[0]:=1;a[1]:=1;

> for k from 0 to 10 do a[k+2]:=-a[k]/((k+2)^2-(k+2)) od:

> l;

Hence the first 11 terms of the series solution is

> c[1]*add(a[2*k]*x^(2*k),k=0..5)+c[2]*add(a[2*k+1]*x^(2*k+1),k=0..4);

We recognize that the sums are Taylor expensions of and respectively. Thus the solution is .

Example 2: Find a series solution in powers of x of Airy's equation y''-xy=0.

Solution:

> y2-x*y;

Using copy and paste multiply the second sum by x

> sum(k*(k-1)*a[k]*x^(k-2),k = 2 .. infinity)-sum(a[k]*x^(k+1),k = 0 .. infinity);

In the first sum replace k-2 with k+1, that is replace k with k+3. Don't forget to correctly adjust the summation limits

> sum((k+3)*(k+2)*a[k+3]*x^(k+1),k = -1 .. infinity)-sum(a[k]*x^(k+1),k = 0 .. infinity);

Write this as a single sum. Notice that the sum has to start from k=0, so there will be an extra term from the first sum that corresponds to k=-1.

> sum(((k+3)*(k+2)*a[k+3]-a[k])*x^(k+1),k = 0 .. infinity)+2*1*a[2];

Hence the recurrence relation is for ..., together with . Notice that to solve this recurrence relation we can set and arbitrarily while has to be zero. We find first 20 terms of the series solution in a similar way as we did in the preway as we did in the previous example.

> restart;

> l:=[seq(a[i],i=0..20)]:

> a[0]:=1;a[1]:=1;a[2]:=0;

> for k from 0 to 20 do a[k+3]:=a[k]/((k+3)*(k+2)) od:

Hence the first 21 terms of the series solution is

> c[1]*add(a[2*k]*x^(2*k),k=0..10)+c[2]*add(a[2*k+1]*x^(2*k+1),k=0..9);

Example 3: Solve Airy's equation in powers of x-1.

Solution:

> restart;

> y:=sum(a[k]*(x-1)^k,k=0..infinity);y1:=sum(k*a[k]*(x-1)^(k-1),k=1..infinity);y2:=sum(k*(k-1)*a[k]*(x-1)^(k-2),k=2..infinity);

>

> y2-x*y;

Since all the sums have to be in powers of (x-1), we first write x as (x-1)+1

> y2-(x-1)*y-y;

Bring (x-1) inside the second sum

> sum(k*(k-1)*a[k]*(x-1)^(k-2),k = 2 .. infinity)-sum(a[k]*(x-1)^(k+1),k = 0 .. i 2 .. infinity)-sum(a[k]*(x-1)^(k+1),k = 0 .. infinity)-sum(a[k]*(x-1)^k,k = 0 .. infinity);

To write this as a single sum we have to fix powers of (x-1). One way to do this is to replace (k-2) with k in the first sum, and to replace (k+1) with k in the second sum. This amounts to replacing k with k+2 in the first sum, and with replacing k with (k-1) in the second sum.

> sum((k+2)*(k+1)*a[k+2]*(x-1)^k,k = 0 .. infinity)-sum(a[k-1]*(x-1)^k,k = 1 .. infinity)-sum(a[k]*(x-1)^k,k = 0 .. infinity);

We write these three sums as a single sum. The common sum has to start at k=1, so there will be to terms that correspond to k=0 from the first and the third sum.

> sum(((k+2)*(k+1)*a[k+2]-a[k-1]-a[k])*(x-1)^k,k = 1 .. infinity)+2*1*a[2]-a[0];

Hence the recurrence relation is for ..., together with .

> restart;

> l:=[seq(a[i],i=0..20)]:

> a[0]:=c[1];a[1]:=c[2];a[2]:=c[1]/2;

> for k from 1 to 20 do a[k+2]:=(a[k-1]+a[k])/((k+1)*(k+2)) od:

Hence the first 21 terms of the series solution is

> <>
> add(a[k]*(x-1)^k,k=0..20);

>

Example 4: Solve Airy's equation in powers of x-1.

Solution: To solve the differential equation y''-xy=0 in powers of (x-1), is the same as to solve the differential equation y''-(t+1)y=0 in powers of t. Notice that all we did is to replace x-1 by t, that is to replace x by t.

> y:=sum(a[k]*t^k,k=0..infinity);y1:=sum(k*a[k]*t^(k-1),k=1..infinity);y2:=sum(k*(k-1)*a[k]*t^(k-2),k=2..infinity);

> y2-t*y-y;

and this is equal to . Hence we obtained exactly the same recurrence relation as before.

Finding Radius of Convergence of Series Solution

The radius of convergence is at least as big as the radius of convergence of p and q. We ilustrate this on the next example.

Example 1. What is the radius of convergence of the Taylor series for about x=0? About x=1?

Solution: First we solve where the fraction is not defined, that is where the denominator is zero.

> solve(x^2-2*x+2=0,x);

Next we find the distance of the solutions from x=0.

> sqrt((1-0)^2+(1-0)^2); sqrt((1-0)^2+(-1-0)^2);

Hence the radius is .

To answer the second part we have to find the distance from x=1.

> sqrt((1-1)^2+(1-0)^2); sqrt((1-1)^2+(-1-0)^2);

Hence the radius of convergence for the Taylor series around x=1 is 1.

Example 2: Determine a lower bound for the radius of convergence of series solutions about x=0 and x=-1/2 for the differential equation .

Solution: Here and . Both p and q are not defined at

> solve(1+x^2=0);

idth=43 height=20 alt="[Maple Math]">

The distance of these points to x=0 are clearly 1 so a lower bound for the radius of convergence of series solution about x=0 is 1.

> sqrt((0-(-1/2))^2+(1-0)^2);sqrt((0-(-1/2))^2+(-1-0)^2);

On the other hand the distance of these points to is so a lower bound for the radius of convergence of series solution about is .

Regular Singular Points

In this section we are looking for solutions of Second Order Linear Equation near singular points. Singular points are those points Singular points are those points where p and q fail to be analytic.

Example: Determine the singular points of the differential equation .

Solution: Here and . Since p and q are rational functions the singular points are those were the denominator is zero

> solve( (x+2)^2*(x-1)=0);

Hence the singular points are x=-2 and x=1. We will consider only the so called regular singular points. A point is a regular singular point if the limits and

are finite.

> limit((x+2)/(x+2)^2,x=-2);

We see that x=-2 is not a regular singular point. On the other hand the both limits

> limit((x-1)/(x+2)^2,x=1);limit((x-1)^2*2*(x+2)/((x+2)^2*(x-1)),x=1);

are finite so x=1 is a regular singular point.

Example: Determine the singular points of .

Solution: The only singular point is . We evaluate the following limits

> limit((x-Pi/2)*cos(x)/(x-Pi/2)^2,x=Pi/2);limit((x-Pi/2)^2*sin(x)/(x-Pi/2)^2,x=Pi/2);

Both limits are finite so is a regular singular point.

Series Solutions near a Regular Singular Point

Theorem:

Consider the differential equation , where is a regular singular point. Then alt="[Maple Math]" align=middle> and are analytic at with convergent power series expansions , for , where is the minimum of the radii of convergence of the power series for and .Let and and be the roots of the indicial equation

,

with if and are real. Then in either of the interval ( ) or ( ), there is a solution of the form

,

where the are given by the recurrence relation

+ [ ]=0, with and . (**)

If is not zero or a positive integer, then in either of intervals ( ) or ( ), there exists a second linearly independent solution of the form

.

The are also determined b00000> are also determined by the recurrence relation ** ,with and . The two power series converge at least for .

If , then the second solution is

.

If , a positive integer, then the second solution is

.

The coefficients , , , and the constant can be determined by substituting the form of the series solution for in the differential equation. The constant may turn out to be zero, in which case there is no logarithmic term in the solution. The two power series converge at least for and defines a function that is analytic in some neighborhood of .

Example:

Determine the indicial equation and the exponents at the singularity for each regular rity for each regular singular point of .

Solution: Clearly x=0 and x=1 are the only singular points, and . Next we check for regularity of these two points.

> limit(x*(-(1+x))/(x^2*(1-x)),x=0);

> limit((x-1)*(-(1+x))/(x^2*(1-x)),x=1);

> limit((x-1)^2*2*x/(x^2*(1-x)),x=1);

Hence the only regular singular point is x=1. The indicial equation is

. The exponents of singularity are

> solve(r*(r-1)+2*r+0 = 0);

Since both roots of the indicial equation are real there is a solution of the form . Since the difference is 1 the second solution is

. We now proceed in finding .

> y:=1+sum(a[n]*(x-1)^n,n = 1 .. infinity);

> y1:=sum(n*a[n]*(x-1)^(n-1),n = 1 .. infinity);

> y2:=sum(n*(n-1)*a[n]*(x-1)^(n-2),n = 2 .. infinity);

> x^2*(1-x)*y2-(1+x)*y1+2*x*y=0;

Since all the powers should be in (x-1) we write as

, (1+x) as (x-1)+2 and x as (x-1)+1.

. Multiply all the sums

. At this point we have to write the sums in the same powers, for example insums in the same powers, for example in the powers of . Don't forget to adjust the indicies correctly.

. Now all the sums are of the same power so we can put them together as a single sum. Notice that the single sum has to start at n=3.

and here I give up.

Bessel's Equation

Bessel's equation is an equation of the form , where is a constant. x=0 is a regular singular point. We consider several possibilites for the valveral possibilites for the value of .

Bessel's Equation of Order Zero,

This is the equation . One can find a solution that is of the form . We will not embark on tedious calculations, we will just say that the solutions are Bessel functions of the first and second kinds, respectively. Maple has built in these functions and calling sequence is BesselJ and BesselY. For example two independent solutions to Bessel's equation of order 0 are BesselJ(0,x) and BesselY(0,x). We check that this is indeed the case.

> x^2*diff(BesselJ(000>x^2*diff(BesselJ(0,x),x,x)+x*diff(BesselJ(0,x),x)+x^2*BesselJ(0,x);

> simplify(%);

Similarly we check that BesselY(0,x) is also a solution to the equation above.

> x^2*diff(BesselY(0,x),x,x)+x*diff(BesselY(0,x),x)+x^2*BesselY(0,x);

> simplify(%);

Hence the general solution is . To get better idea how these functions look like we will graph them.

> p1:=plot(BesselJ(0,x),x=0..14,color=red):

> p2:=plot(BesselY(0,x),x=0..14,color=green):

>
> with(plots):

> display({p1,p2});

Bessel's Equation of Order Zero,

Solve . We recognize this equation as Bessel's equation of order 1/2. Thus two independent solutions are BesselJ(1/2,x) and BesselY(1/2,x). We check that these are solutions as follows.

> x^2*diff(BesselJ(1/2,x),x,x)+x*diff(BesselJ(1/2,x),x)+(x^2-1/4)*BesselJ(1/2,x);

> simplify(%);

> x^2*diff(BesselY(1/2,x),x,x)+x*diff(BesselY(1/2,x),x)+(x^2-1/4)*BesselY(1/2,x);

> simplify(%);

These two Bessel functions have closed form, as we can see below

> BesselJ(1/2,x);

> BesselY(1/2,x);

Bessel's equation of order 1

Example: Solve . We recognize this equation as Bessel's equation of order 1. Thus two independent solutions are BesselJ(1,x) and BesselY(1,x). We check that these are solutions as follows.

> > x^2*diff(BesselJ(1,x),x,x)+x*diff(BesselJ(1,x),x)+(x^2-1)*BesselJ(1,x);

> simplify(%);

> x^2*diff(BesselY(1,x),x,x)+x*diff(BesselY(1,x),x)+(x^2-1)*BesselY(1,x);

> simplify(%);

Laplace Transform

Example 1: Solve with initial conditions and .

Solution: We use Laplace transform to solve this differential equation. First we define initial conditions:

> y(0):=2;D(y)(0):=1;

Next we take Laplace transform of the differential equation

> laplace(diff(y(x),x,x)+y(x)=sin(2*x),x,s);

We solve this equation as follows

> solve(%,laplace(y(x),x,s));

Finally we use invlaplace to find the inverse Laplace

> invlaplace(%,s,x);

Example 2. Solve where with initial conditions and .

Solution: We use Laplace transform to solve this differential equation. First we define initial conditions:

> y(0):=0;D(y)(0):=0;g:=t-> Heaviside(t-5)-Heaviside(t-20);

Next we take Laplace transform of the differential equation

> laplace(2*diff(y(x),x,x)+diff(y(x),x)+2*y(x)=g(x),x,s);

We solve this equation as follows

> solve(%,laplace(y(x),x,s));

Finally we use invlaplace to find the inverse Laplace

>

> invlaplace(%,s,x);

> p:=unapply(%,x);

> plot(p,0..40);

>

Systems of First Order Linear Equations

In order to solve a system, of first order linear equations with constant coefficients we have to find eigenvalues and eigenvectors of the matrix . We will use Maple to do this. To access Maple commands we first need to load Linear Algebra commands.

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

Example 1: Solve the system .

> A:=matrix(2,2,[-3,sqrt(2),sqrt(2),-2]);

> eigenvals(A);

To find eigenvectors we need the column vector and the identity matrix . Notice below how we defined .

> t:=matrix(2,1,[t1,t2]);t:=matrix(2,1,[t1,t2]);

> E:=array(1..2,1..2,identity);

We obtain eigenvectors for the egenvalue -1 by solving the system

.

> multiply(A+E,t);

> solve(-2*t1+2^(1/2)*t2=0,t2);

Hence all eigenvectors of the eigenvalue -1 are of the form .We repeat the same procedure for the eigenvalue -4.

> multiply(A+4*E,t);

> > solve(t1+2^(1/2)*t2=0,t2);

Hence all eigenvectors of the eigenvalue -4 are of the form .

Since the eigenvalues are disitinct by setting t1=1 we obtain two independent eigenvectors and .Hence two independent solutions of the system are and .

Example 2: Solve the system .

Solution:

> A:=matrix(3,3,[0,1,1,1,0,1,1,1,0]);

> E:=array(1..3,1..3,identity);

> t:=matrix(3,1,[t1,t2,t3]);

> eigenvals(A);

> multiply(A-2*E,t);

We know that the eigenvalue 2 has only one eigenvector. We set t1=t2=1 and we find t3 as follows

> solve(-2+1+t3=0,t3);

Hence one eigenvector is and the corresponding solution is

. On the other hand eigenvalue -1 has two independent eigenvectors.

> multiply(A+E,t);

> solve({t1+t2+t3=0,t1+t2+t3=0,t1+t2+t3=0},{t1,t2,t3});

Hence there are two independetn eigenvectors which we can obtain by setting first t2=0 and t3=-1, and tehn by setting t2=1 and t3=-1..

Hence two independent eigenvectors that correspond to -1 are and and thus two independent solutions of the differntial equation are

and .

In the case when eigenvalues are complex we have to find real and imaginary parts of the complex-valued solutions which we obtain in the same way as in the real case.

Example 3. Solve .

Solution:

> A:=matrix(2,2,[-1/2,1,-1,-1/2]);

> E:=array(1..2,1..2,identity);

> t:=matrix(2,1,[t1,t2]);

> eigenvalues(A);

> multiply(A-(-1/2+I)*E,t);

We set t1=1, and solve -It1+t2=0 in t2.

> solve(-I+t2=0,t2);

Hence one eigenvector is and the corresponding solution is

. This is a complex-valued solution so we have to find the real and imaginary part. Rewrite as . Therefore we can write the complex solution in the form

. Hence two real and independent solutions are

and .

Example 4: Solve the system .

Solution: As before, we define the matrix , the vector and the identity matrix .

>
> A:=matrix(2,2,[1,- 1,1, 3]);
> t:=matrix(2,1,[t1,t2]); > E:=array(1..2,1..2,identity); Next we find eigenvalues of . > eigenvalues(A); Notice that the eigenvalue 2 is a repeated eigenvalue, which will yield only one independent solution. We find this eigenvector in the same way we did before. > multiply(A-2*E,t); > solve(-t1-t2=0,t2); Hence one eigenvector is which gives one independent solution . We find a second solution as follows: First solve the matrix equation . > multiply(A-2*E,t); > solve({-t1-t2=1,t1+t2=-1},{t1,t2}); > Set t1=-1 (It is the same value of t2 as in the eigenvector .), and find t1. Thus one solution of the matrix equation is . Now the second solution of the differential equation is .Example 5: Solve the system . Solution: > A:=matrix(3,3,[5,-3,-2,8,-5,-4,-4,3,3]); > E:=array(1..3,1..3,identity); > t:=matrix(3,1,[t1,t2,t3]); > <> > eigenvalues(A); Here, 1 one is a triple eigenvalue. > multiply(A-1*E,t); > solve({4*t1-3*t2-2*t3=0,8*t1-6*t2-4*t3=0,-4*t1+3*t2+2*t3=0},{t1,t2,t3}); The eigenvalue 1 has two independent eigenvectors, and we could obtain them by setting first t1=1 and t2=0, and then by setting t1=0 and t2=1. Hence two independent egienvectors corresponding to the eigenvalue 1 are and ., and the corresponding two independent solutions of the differential equation are and . We find a third independent solution in the similar way as we did it in Example 4. We have to solve the matrix equation . This system does not have a solution for all choices of and . The second row of is a multiply of the first row, so on the right hand side of the equation we should also have that the second row be twice the first, that is the system will have a solution only if . Now set . > solve({4*t1-3*t2-2*t3=1,8*t1-6*t2-4*t3=2,-4*t1+3*t2+2*t3=-1},{t1,t2,t3}); Thus one solution of the matrix equation that corresponds to and is obtained by setting t1=0 and t2=0 to get t3=-1/2. Hence a third independent solution of the differential equation is that is .