Application Center - Maplesoft

# Eigenvalue problems for ordinary differential operators

You can switch back to the summary page by clicking here.

eigenval.mws

Eigenvalue problems for ordinary differential operators

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

NOTE: In this session we solve eigenvalue problems for second order ordinary differential operators.

1. -y''= lambda*y

Problem

 > restart;

 > L:=-diff(y(x),x,x)=lambda*y(x);

Please  input  number of examples  n (1..4) and     execute Section

 > n:=4;

 > if n=1 then   s1:=y(a)=0; s2:=y(b)=0;

 > elif n=2 then s1:=D(y)(a)=0; s2:=y(b)=0;

 > elif n=3 then s1:=y(a)=0; s2:=D(y)(b)=0;

 > elif n=4 then s1:=D(y)(a)=0; s2:=D(y)(b)=0; end if;

 > assume(a

Solving problem

 > assume(lambda>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > linalg[genmatrix](sist,{_C1,_C2});

 > linalg[det](%);

 > combine(%);

 > if type(%,function) then chl:=% else

 > chl:=select(has,%,[sin,cos]) end if;

 > _EnvAllSolutions := true:

 > solve(chl,lambda): sort(%);

 > subs(1/(a-b)^2=1/(b-a)^2,%):sort(%,[b,a]);

 > indets(%) minus {a,b};

 > subs(%[1]='k',%%):

Eigenvalues:

 > ev:=unapply(%,'k');

 > y(x);

 > isolate(s1,_C1);

 > subs(%,lambda=ev(k),y(x)):

 > simplify(combine(%));

 > select(has,%,x): subs(a-b=b-a,-x+a=x-a,%): sort(%,[x,b,a]);

 > %/sqrt(int(%^2,x=a..b)):sort(%,[x,b,a]);

Eigenfunctions:

 > ef:=unapply(%,x,k);

 > y:='y':

Example

 > cat(n,` EXAMPLE `); L, s1, s2;

 > lambda[k]=ev(k),y[k]=ef(x,k);

Checking the Solution

Checking  equation:

 > L;

 > subs(y(x)=ef(x,k),lambda=ev(k),L):

 > simplify(lhs(%)-rhs(%));

Checking  boundary conditions:

 > s1,s2;

 > y:=unapply(ef(x,k),x):

 > s1,s2;

 > y:='y':

 >

2.  -x^2*y''-y/4 = lambda*y

Problem

 > restart;

 > L:=-x^2*diff(y(x),x,x)-y(x)/4=lambda*y(x);

 > s1:=y(a)=0;s2:=y(b)=0;

 > a:=1;b:=exp(1);

Solving problem

 > assume(lambda>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > linalg[genmatrix](sist,{_C1,_C2});

 > linalg[det](%);

 > combine(%);

 > if type(%,function) then chl:=% else

 > chl:=select(has,%,[sin,cos]) end if;

 > _EnvAllSolutions := true:

 > solve(chl,lambda);

 > indets(%) minus {a,b};

 > subs(%[1]='k',%%):

Eigenvalues:

 > ev:=unapply(%,'k');

 > y:='y':

 > lyg:=subs(lambda=ev(k),L);

 > assume(k,posint);

 > dsolve({lyg,op(0,lhs(s1))(a)=0,op(0,lhs(s2))(b)=0},{y(x)});

 > subs(_C1=1,_C2=1,rhs(%));

 > %/sqrt(int(%^2,x=a..b));

Eigenfunctions:

 > ef:=unapply(%,x,k);

 >

Solution

 > L, s1, s2;

 > lambda[k]=ev(k),y[k]=ef(x,k);

Checking the Solution

Checking  equation:

 > L;

 > subs(y(x)=ef(x,k),lambda=ev(k),L):

 > simplify(lhs(%)-rhs(%));

Checking  boundary conditions:

 > s1,s2;

 > y:=unapply(ef(x,k),x):

 > s1,s2;

 > y:='y':

 >

3.  -x^2*y''-x*y'=lambda*y

Problem

 > restart;

 > L := -x^2*diff(diff(y(x),x),x)-x*diff(y(x),x)=lambda*y(x);

 > s1:=y(a)=0;s2:=y(b)=0;assume(a>0,b>a):

 > #a:=1;b:=2;

Solving problem

 > assume(lambda>0);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > linalg[genmatrix](sist,{_C1,_C2});

 > linalg[det](%);

 > combine(%);

 > if type(%,function) then chl:=% else

 > chl:=select(has,%,[sin,cos]) end if;

 > _EnvAllSolutions := true:

 > solve(chl,lambda);

 > indets(%,symbol) minus {Pi,a,b};

 > subs(%[1]=k,%%):

Eigenvalues:

 > ev:=unapply(%,k);

 > y:='y':

 > lyg:=subs(lambda=ev(k),L);

 > assume(k,posint);

 > tr:=x=expand(exp(t*(ln(b)-ln(a))+ln(a)));

 > _EnvAllSolutions := false:

 > itr:=solve(tr,{t});

 > PDEtools[dchange](tr,lyg,params=[a,b],simplify);

 > dsolve({%,y(0)=0,y(1)=0},y(t));

 > subs(itr,_C1=1,_C2=1,rhs(%));

Eigenfunctions(not normed):

 > ef:=unapply(%,x,k);

 >

Solution

 > L, s1, s2;

 > lambda[k]=ev(k),y[k]=ef(x,k);

Checking the Solution

Checking  equation:

 > L;

 > subs(y(x)=ef(x,k),lambda=ev(k),L):

 > simplify(lhs(%)-rhs(%));

Checking  boundary conditions:

 > s1,s2;

 > y:=unapply(ef(x,k),x):

 > s1,simplify(s2);

 > y:='y':

 >

4.  -y''-2*y'=lambda*y

Problem

 > restart;

 > L := -diff(diff(y(x),x),x)-2*diff(y(x),x)=lambda*y(x);

 > s1:=y(a)=0;s2:=y(b)=0;

 > a:=1;b:=2;

Solving problem

 > assume(lambda>1);

 > dsolve(L,y(x));

 > y :=unapply( rhs(%),x);

 > sist:={s1,s2};

 > linalg[genmatrix](sist,{_C1,_C2});

 > r:=linalg[det](%);

 > chl:=combine(%);

 > _EnvAllSolutions := true:

 > solve(chl,lambda);

 > indets(%,symbol) minus {Pi,a,b};

 > subs(%[1]=k,%%):

Eigenvalues:

 > ev:=unapply(%,k);

 > y:='y':

 > lyg:=subs(lambda=ev(k),L);

 > assume(k,posint);

 > dsolve({lyg,y(a)=0,y(b)=0},y(x));

 > subs(_C1=1,_C2=1,rhs(%));

 > simplify(%/sqrt(int(%^2,x=a..b)));

Eigenfunctions:

 > ef:=unapply(%,x,k);

 >

Solution

 > L, s1, s2;

 > lambda[k]=ev(k),y[k]=ef(x,k);

 >

Checking the Solution

Checking  equation:

 > L;

 > subs(y(x)=ef(x,k),lambda=ev(k),L):

 > simplify(lhs(%)-rhs(%));

Checking  boundary conditions:

 > s1,s2;

 > y:=unapply(ef(x,k),x):

 > simplify(s1),simplify(s2);

 > y:='y':

 >

5.   -y''-y'/x=lambda*y (Bessel functions)

Problem

 > restart;

 > eq:=-diff(y(x),x,x)-diff(y(x),x)/x=lambda*y(x);

 > abs(y(0))

 > assume(R>0);

 >

Solving problem

 > Order:=10:Digits:=20:

 > sol:=dsolve({eq,y(0)

 > cheq:=subs(x=R,diff(rhs(sol),x));

 > solve(cheq,lambda);

 > mu:=[BesselJZeros(1.,0..Order)];

 > #plot(BesselJ(1,x),x=0..20);

Eigenvalues:

 > ev:=[seq(solve(sqrt(lambda)*R=mu[k],lambda),k=1..Order)];

Eigenfunctions:

 > ef:=[seq(simplify(subs(lambda=ev[k],y(0)=1,rhs(sol))),k=1..Order)];

Solution

 > for k to Order do print(lambda[k]=ev[k],y[k]=ef[k]) od;

Checking

Checking  equation:

 > for k to Order do

 > expand(simplify(subs(lambda=ev[k],y(x)=ef[k],lhs(eq)-rhs(eq))));

 > od;

Checking  boundary condition:

 > for k to Order do expand(subs(x=R,diff(ef[k],x))) od;

 >

6. Gegenbauer (ultraspherical) polynomials G(n,a,x)

Problem

 > restart;with(orthopoly):

 > eq := -(1-x^2)*diff(y(x),x,x)+(2*a+1)*x*diff(y(x),x)=lambda*y(x);

 > abs(y(-1))

Solving problem

Eigenfunctions we search in polynomials  f,    .

 > Geg:=proc(n::nonnegint,a) local f,m,id,eq,cf,koef,F,K; eq := (1-x^2)*diff(y(x),x,x)-(2*a+1)*x*diff(y(x),x)+lambda*y(x) = 0; f:=binomial(n+2*a-1,n)+sum('A[k]*x^k','k'=0..n-1)*(x-1); subs(y(x)=f,eq); id:=expand(lhs(%)-rhs(%))=0; K:=indets(id) minus {x}; koef:=solve(identity(id,x),K); m:=nops([koef]); F:=(a,b)->evalb(degree(subs(a,f))>degree(subs(b,f))); cf:=sort([koef],F); RETURN([subs(cf[1],lambda),sort(expand(subs(cf[1],f)))]); end:

Solution

We assign     for example :

 > a:=1/2:

 > for k from 0 to Order do lambda[k]=Geg(k,a)[1], y[k]=Geg(k,a)[2]; end do;

Eigenfunctions is Gegenbauer    polynomials:

 > seq(orthopoly[G](k,a,x),k=0..Order);

and     :

 > seq(n*(2*a+n),n=0..Order);

Checking the Solution

 > seq(simplify(subs(y(x)=orthopoly[G](k,a,x), lambda=k*(2*a+k),lhs(eq)-rhs(eq))),k=0..Order);

 >

Other relations

1. Generating function for Gegenbauer polynomials  is

 > map(simplify,series((1-2*s*x+s^2)^(-a),s));

2.

 > seq(simplify(binomial(n+2*a-1,n)*hypergeom([n+2*a,-n],[a+1/2],(1-x)/2)), n=0..Order-1);

 >

7. Hermite polynomials  H(n,x)

Problem

 > restart;with(orthopoly,H,T):

 > eq := -diff(y(x),x,x)+2*x*diff(y(x),x)=lambda*y(x);

 > abs(y(0))

 >

Solving problem

Eigenfunctions we search in polynomials  f,   + ...    .

 > Her:=proc(n::nonnegint) local f,id,K; global eq; f:=2^n*x^n+sum('A[k]*x^k','k'=0..n-1); subs(y(x)=f,eq); id:=expand(lhs(%)-rhs(%))=0; K:=indets(id) minus {x}; solve(identity(id,x),K); subs(%,f); RETURN([subs(%%,lambda),sort(%)]); end:

 > seq(Her(i),i=0..Order);

Solution

 > for k from 0 to Order do lambda[k]=Her(k)[1], y[k]=Her(k)[2]; end do;

Eigenfunctions is Hermite  polynomials:

 > seq(orthopoly[H](k,x),k=0..Order);

and     :

 > seq(2*n,n=0..Order);

 > k:='k':

Checking the Solution

 > seq(expand(value(subs(y(x)=orthopoly[H](k,x), lambda=2*k,lhs(eq)-rhs(eq)))),k=0..Order);

Other relations

1.

:

Checking :

 > convert(series(exp(-t^2+2*t*x),t),polynom):

 > sum(H(k,x)*t^k/k!,k = 0 .. Order-1):

 > simplify(%-%%);

2.

Checking :

 > simplify([seq((-1)^n*(2*n)!/n!*hypergeom([-n],[1/2],x^2)-H(2*n,x), n=0..Order)]);

 > simplify([seq((-1)^n*(2*n+1)!*2/n!*x*hypergeom([-n],[3/2],x^2)-H(2*n+1,x), n=0..Order)]);

3. If  x>=0 then

Checking :

 > seq(expand(2^n*KummerU(-n/2,1/2,x^2)),n=0..4):

 > plot([%],x=0..3,y=-20..20);

 > seq(H(n,x),n=0..4):

 > plot([%],x=0..3,y=-20..20);

 >

8.   Laguerre polynomials L(n,x)

Problem

 > restart;with(orthopoly):

 > eq := -x*diff(y(x),x,x)+(x-1)*diff(y(x),x)=lambda*y(x);

 > abs(y(0))

Solving problem

Eigenfunctions we search in polynomials  f,   .

 > n:=Order;

 > f:=1-sum(A[k]*x^k,k=0..n-1)*x;

 > subs(y(x)=f,eq);

 > id:=expand(lhs(%)-rhs(%))=0;

 > K:=indets(id) minus {x};

 > koef:=solve(identity(id,x),K);

 > m:=nops([koef]);

 > F:=(a,b)->evalb(degree(subs(a,f))

 > cf:=sort([koef],F);

 >

Solution

 > for k to m do lambda[k-1]=subs(cf[k],lambda), y[k-1]=simplify(subs(cf[k],f)) od;

Eigenfunctions is Laguerre    polynomials:

 > seq(orthopoly[L](k,x),k=0..Order);

and     :

 > seq(n,n=0..Order);

 > k:='k':

Tes t

 > seq(simplify(value(subs(y(x)=orthopoly[L](k,x), lambda=k,lhs(eq)-rhs(eq)))),k=0..Order);

 >

Other relations

1.

Checking :

 > convert(series(exp(-x*t/(1-t))/(1-t),t),polynom):

 > sum(L(k,x)*t^k,k = 0 .. Order-1):

 > simplify(%-%%);

2.

:

Checking :

 > seq(simplify(hypergeom([-n],[1],x)-L(n,x)),n=0..Order);

 >

9. Legendre polynomials P(n,x)

Problem

 > restart;with(orthopoly):

 > eq := -diff((1-x^2)*diff(y(x),x),x)=lambda*y(x);

 > abs(y(-1))

Solving problem

Eigenfunctions we search in polynomials  f,   .

 > n:=Order;

 > f:=1-sum(A[k]*x^k,k=0..n-1)*(x-1);

 > subs(y(x)=f,eq);

 > T:=expand(lhs(%)-rhs(%))=0;

 > K:=indets(T) minus {x};

 > koef:=solve(identity(T,x),K);

 > m:=nops([koef]);

 > F:=(a,b)->evalb(degree(subs(a,f))

 > cf:=sort([koef],F);

 >

Solution

 > for k to m do lambda[k-1]=subs(cf[k],lambda), y[k-1]=simplify(subs(cf[k],f)) od;

Eigenfunctions is Legendre polynomials:

 > seq(orthopoly[P](k,x),k=0..Order);

and     :

 > seq(n*(n+1),n=0..Order);

 >

Checking the Solution

 > seq(simplify(subs(y(x)=orthopoly[P](k,x), lambda=k*(k+1),lhs(eq)-rhs(eq))),k=0..Order);

 >

Other relations

1. Generating function for Legendre polynomials  is    :

 > simplify(series(1/sqrt(1-2*t*x+t^2),t));

2.

:

 > seq(simplify(hypergeom([n+1, -n],[1],(1-x)/2)),n=0..Order-1);

 >

10. Chebyshev polynomials T(n,x)

Problem

 > restart;with(orthopoly):

 > eq := -(1-x^2)*diff(y(x),x,x)+x*diff(y(x),x)=lambda*y(x);

 > abs(y(-1))

Solving problem

Eigenfunctions we search in polynomials  f,   .

 > n:=Order;

 > f:=1-sum(A[k]*x^k,k=0..n-1)*(x-1);

 > subs(y(x)=f,eq);

 > id:=expand(lhs(%)-rhs(%))=0;

 > K:=indets(id) minus {x};

 > koef:=solve(identity(id,x),K);

 > m:=nops([koef]);

 > F:=(a,b)->evalb(degree(subs(a,f))

 > cf:=sort([koef],F);

 >

Solution

 > for k to m do lambda[k-1]=subs(cf[k],lambda), y[k-1]=simplify(subs(cf[k],f)) od;

Eigenfunctions is Chebyshev polynomials:

 > seq(orthopoly[T](k,x),k=0..Order);

and     :

 > seq(n^2,n=0..Order);

 >

Checking the Solution

 > seq(simplify(subs(y(x)=orthopoly[T](k,x), lambda=k^2,lhs(eq)-rhs(eq))), m=0..Order);

Other relations

1. Generating function for Chebyshev polynomials  is  :

 > simplify(series((1-t*x)/(1-2*t*x+t^2),t));

2.

:

 > seq(simplify(hypergeom([n,-n],[1/2],(1-x)/2)),n=0..Order-1);

 >

11. Chebyshev polynomials of the second kind   U(n,x)

Problem

 > restart; with(orthopoly):

 > eq := -(1-x^2)*diff(y(x),x,x)+3*x*diff(y(x),x)=lambda*y(x);

 > abs(y(-1))

Procedure

Eigenfunctions we search in polynomials  f,

 > Ch2:=proc(n::nonnegint) local f,id,eq,cf,F,K; eq := -(1-x^2)*diff(y(x),x,x)+3*x*diff(y(x),x)=lambda*y(x); f:=n+1+sum('A[k]*x^k','k'=0..n-1)*(x-1); subs(y(x)=f,eq); id:=expand(lhs(%)-rhs(%))=0; K:=indets(id) minus {x}; cf:=solve(identity(id,x),K); F:=(a,b)->evalb(degree(subs(a,f))>degree(subs(b,f))); cf:=sort([cf],F); RETURN([subs(cf[1],lambda),expand(subs(cf[1],f))]); end:

Solution

 > seq(Ch2(k),k=0..Order);

 > for k from 0 to Order do lambda[k]=Ch2(k)[1], y[k]=Ch2(k)[2] od;

Eigenfunctions is Chebyshev polynomials of the second kind:

 > seq(orthopoly[U](k,x),k=0..Order);

and     :

 > seq(n*(n+2),n=0..Order);

Checking the Solution

 > seq(simplify(subs(y(x)=orthopoly[U](k,x), lambda=k*(k+2),lhs(eq)-rhs(eq))),k=0..Order);

Other relations

1. Generating function for Chebyshev polynomials of the second kind is  :

 > simplify(series(1/(1-2*t*x+t^2),t));

2.

:

 > seq(simplify((n+1)*hypergeom([n+2,-n],[3/2],(1-x)/2)),n=0..Order-1);

3.

:

 > seq(diff(T(n+1,x),x)/(n+1),n=0..Order-1);

12. Associated Legendre functions

Problem

 > restart;

 > eq := -(1-x^2)*diff(y(x),x,x)+2*x*diff(y(x),x)+m^2/(1-x^2)*y(x)=lambda*y(x);

 > abs(y(-1))

 > assume(m,nonnegint);

 >

Solution

Eigenfunctions are associated Legendre functions

,   , ...

where    Legendre polynomials and

We define  :

 > P:=proc(n::nonnegint,m::nonnegint,x) local i; if n

For m=0     :

 > seq(P(n,0,x),n=0..5);

For m=1

 > seq(P(n,1,x),n=1..5);

For m=2

 > seq(P(n,2,x),n=2..6);

For m=3

 > seq(P(n,3,x),n=3..6);

 >

Checking the Solution

For m=1

 > seq(simplify(subs(y(x)=P(n,1,x),m=1, lambda=n*(n+1),lhs(eq)-rhs(eq))), n=1..Order);

For m=4

 > seq(simplify(subs(y(x)=P(n,4,x),m=4, lambda=n*(n+1),lhs(eq)-rhs(eq))), n=4..2*Order);

 >

Other relations

1.

Generating function

 > gf:=(2*m)!*(1-x^2)^(m/2)/(2^m)/m!/((1-2*x*t+t^2)^(m+1/2));

Checking  for m=1

 > gf1:=subs(m=1,gf);

 > series(gf1,t,5);

 > seq(P(n,1,x),n=1..5);

Checking  for m=2

 > gf2:=subs(m=2,gf);

 > series(gf2,t,5);

 > seq(P(n,2,x),n=2..Order);

2.

Checking :

 > P1:=proc(n::nonnegint,m::nonnegint,x) if n

For m=1

 > seq(simplify(P(n,1,x)-P1(n,1,x)),n=1..Order);

For m=2

 > seq(simplify(P(n,2,x)-P1(n,2,x)),n=2..Order);

 >

Note

In some references associated Legendre functions are defined by

see  in Maple ?LegendreP,  in Mathematica LegendreP[n,m,x] .

Examples with associated Legendre functions see in elliptic3d.mws

 >

While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.