本征函数系与本征振动

第四章本征函数系与本征振动一维本征值问题；正余弦函数系二维本征值问题；贝塞尔函数系三维本征值问题；球谐函数系(勒让德函数系)用偏微分方程工具箱(pdetool)求解本征值问题；4.0 PDE分类2

2Laplacian

算子：

x2

y2拉（elliptic）：斯方程

a2u

f热传导方程(parabolic)：tu

a2u

f波动方程(hyperbolic):a

u

f22t2u

初边值问题边值问题设函数u(x,t)

具有变量分离形式，即它可表示为：或则上述三个方程可以写为：4.1

一维本征值问题对于齐次问题，f=0，在一(二)维情况下，三种方程分别变为：二维拉

斯方程

uxx

uyy

0 (0

x

a;0

y

b)一维热传导方程ut

a

u

0

(0

x

l)2xxutt

a

u

0

(0

x

l)2xx一维波动方程u(x,

y)

X

(x)Y

(

y)

u(x,t)

X

(x)T

(t)这些方程的分离变量中都会出现方程：X

X

0该方程与各种边值条件结合可以组成各种本征函数系二维拉斯方程

Y

Y

0;X

X

0;一维热传导方程2

Ta2ta

T

0

T

e

X

X

0一维波动方程02T

a

T

0

T

cos

at

X

X

04.1.1

四种常见的正、余弦本征函数系X

x

sin

nπx

/

l

X

0

0;

X

l

0;

nπ

/

l

2

;

n

1,2,...X

X

01.

sin(nπx/l)X

0

0;

X

l

0;

nπ

/

l2.

cos(nπx/l)

X

X

0X

x

sin

n

1/

2πx

/

lX

0

0;

X

l

0;

n

1/

22

π2

/

l

2

;

n

0,1,...X

X

03.

sin[(n+1/2)πx/l]X

x

cosn

1/

2πx

/

lX

0

0;

X

l

0;

n

1/

22

π2

/

l

2

;

n

0,1,...X

X

04.

cos[(n+1/2)πx/l]4.1.2

本征函数系的图像及其运动上节中的四个本征函数的图像如下页图。显然，端点为第一类边界条件时，该端点为波节，端点为第二类边界条件时，该端为波腹。程序如下：%ex301(p65)clear;x=pi*(0:0.001:1);A=sin([1:4]'*x);

B=cos([0:3]'*x);C=sin([1/2:7/2]'*x);

D=cos([1/2:7/2]'*x);figure(1);subplot(4,1,1);plot(x,A);subplot(4,1,2);plot(x,B);subplot(4,1,3);plot(x,C);subplot(4,1,4);plot(x,D);实际上对于波动方程和热传导方程，还要考虑时间本征函数：对于波动方程，时间函数为cos(nπa

t/l)或cos[(n+1/2)πa

t/l)对于热传导方程，时间函数为exp(-λa2t)

(文中该处p66有误)考虑时间因子后的动画分别如程序ex302和ex303所示。显然，对于波动方程，将形成驻波，固定端(第一类边界条件)为波节，开口端(第二类边界条件)为波腹。对于热传导方程，若两端温度相同(第一类边界条件)，对于任何初始温度分布，最终都将趋于平衡态，%ex302(p65)function

zbt=0:0.005:2.0;

x=0:0.001:1;w1=wf

(1,0); w2=wf(2,0);

w3=wf

(3,0);

w4=wf(4,0);figure(1);subplot(2,1,1);h1=plot(x,w1,'r');hold

on;

h3=plot(x,w3,'g');grid

on;

legend(‘k=1’,‘k=3’);

title('驻波k=1、3');ym1=max(abs(w1));axis([0,1,-ym1,ym1]);subplot(2,1,2);h2=plot(x,w2,'r');hold

on;

h4=plot(x,w4,'g');grid

on;

legend('k=2','k=4');

title('驻波k=2、4');ym4=max(abs(w4));axis([0,1,-ym4,ym4])

;forn=2:length(t)w1=wf(1,t(n));

w2=wf(2,t(n));

w3=wf(3,t(n));

w4=wf(4,t(n));drawnow;drawnow;set(h1,'ydata',w1);

set(h3,'ydata',w3);set(h2,'ydata',w2);

set(h4,'ydata',w4);end;function

p=wf(k,t)x=0:0.001:1;

a=1;

p=sin(k*pi*x).*cos(k*pi*a*t);%ex303(p67)functionsfbx=0:0.01:1;

t=10^-5*(0:10:2000);A=sf(1,t(1));

B=sf(2,t(1));

C=sf(3,t(1));

D=sf(4,t(1));figure(1);subplot(2,1,1);h1=plot(x,A,'r');hold

on;

h3=plot(x,B,'g');grid

on;y1=max(abs(A));axis([0,1,-y1,y1]);subplot(2,1,2);h2=plot(x,B);hold

on;

h4=plot(x,D,'g');gridon;y2=max(abs(B));axis([0,1,-y2,y2]);forn=2:length(t)A=sf(1,t(n));

set(h1,'ydata',A);drawnow;C=sf(3,t(n));

set(h3,'ydata',C);B=sf(2,t(n));

set(h2,'ydata',B);D=sf(4,t(n));

set(h4,'ydata',D);drawnow;endfunction

wtx=sf(k,t)x=0:0.01:1;

wtx=sin(k*pi*x)*exp(-(2*pi*k)^2*t);2

22

ux2

y2

u

u

ux,y;22yxX

Y

kkxyXYx,yuX

Y

令

x

X

x

sin

nπx

/

aX

0

X

a

0X

k

2

Xu0,

y

ua,

y

ux,0

ux,

b

0;可通过分离变量求得本征函数：y

Y

y

sin

mπy

/

b

0

Y

b

0YY

k

2

y二维本征值问题矩形区域的本征模与本征振动具有固定边界的矩形区域拉 斯方程的本征值问题：ux,

y

sin

nπx

/

asin

mπy

/

b;

nm

nm

k

2

nπ

/

a2

mπ

/

b2

;其前四个模式如下页图所示。程序如下：%ex304(p70)a=2;b=1;

[m,n]=meshgrid(1:3);L=((n*pi/a).^2+(m*pi/b).^2)x=0:0.05:2;

y=0:0.025:1;

[X,Y]=meshgrid(x,y);w1=sin(pi*X/a).*sin(pi*Y/b);

w2=sin(pi*X/a).*sin(2*pi*Y/b);w3=sin(2*pi*X/a).*sin(pi*Y/b);

w4=sin(3*pi*X/a).*sin(pi*Y/b);figure(1);subplot(2,2,1);mesh(X,Y,w1);

subplot(2,2,2);mesh(X,Y,w2);subplot(2,2,3);mesh(X,Y,w3);

subplot(2,2,4);mesh(X,Y,w4);在波动问题中，还要考虑时间因子：T

t

sin

kmnat

0

相应的动画演示程序为ex305ex304%ex305(p71)a0=1;

a=2;b=1;x=0:0.05:a;

y=0:0.025:b;

[X,Y]=meshgrid(x,y);form=1forn=1:3fori=1:50k=sqrt((m*pi/a)^2+(n*pi/b)^2);

t=i*0.02;Z=sin(k*a0*t)*sin(n*pi*X/a).*sin(m*pi*Y/b);T=['本征振动:','m=',int2str(m),',n=',int2str(n)];mesh(X,Y,Z);

axis([0,a,0,b,-1,1]);

title('T');p(:,i)=getframe;movie(p);end;end;end;0222

u

0;

k

u;1

u

1

2u

2

u

的第n个零点022Rmnm,n

m

mmnm

mim其中k

n

xm

/

;

xm为m阶贝塞尔函数m

n

0

n

u

,

J

k r

e

R

J

k r

2

m2

0

Rk

0

2

R

R

k

2

m2

R

0

k

R

R

4.2.2

园形区域的本征模与本征振动具有固定边界的圆形区域拉 斯方程的本征值问题：可通过分离变量求得本征函数：令u,

R

ann)函数%ex106(p32-33)

;%贝塞耳(Bessel)函数\诺伊曼(clear;k=(0:3)';

m=20;

x=.05:0.1:m;figure(1);subplot(2,1,1);y=besselj(k,x);plot(x,y);grid

on;

axis([0,m,-.5,1]);title('0-3阶Bessel函数');subplot(2,1,2);y=bessely(k,x);plot(x,y);grid

on;

axis([0,m,-1,.6]);title('0-3阶

ann函数');ex306图3-10

(m=0,n=0~3)ex307图3-14

(m=1,n=0~3)ex308图3-16

(m=2,n=0~3)ex30922

22r

k

u;r

r1

2

u

1

u

1

2u

2u

sin

sin

r

sin

sin

1

sin

1

222222m

l

l

Y

1

,

Pl

,Y

1

YY

sin

sin

Rr

c1

jl

kr

c2n

l

krr

2

R

2rR

kr2

ll

1R

0

krY

1

YRY

sin

sinr

R

2rR三维本征值问题(球贝塞耳函数和勒让德函数)球坐标系中的拉 斯方程在球坐标系中拉 斯方程的本征值问题：可通过分离变量求得本征函数：

令ur,,

其中称为l

阶连带(或缔合)勒让德函数，在中，用legendre(l,x)表示。它们的图像分别如下几页所示。jl(r)与nl(r)分别称为l阶球贝塞耳函数和l

阶球诺伊曼函数；它们与半数阶贝塞耳函数Jl+1/2(r)和半数阶诺伊曼函数Yl+1/2(r)的关系为：

immlmmlPYcos

el

m

!

4πl

m

!2l

1

,

1mlP

cos

l而Y

m(θ,φ)称为球谐函数，可以写为：

nl

π

πY

r2r

l

1/

2J2r

l

1/

2r

;

r

jl

r

在 中，l

阶球贝塞耳和

l

阶球诺伊曼函数分别用：sqrt(pi/2/r)*besselj(l+1/2,r)

与sqrt(pi/2/r)*

bessely(l+1/2,r)表示。%ex106(p32-33);Bessel\诺伊曼函数与球Bessel\诺伊曼函数clear;

k=0:3;m=20;x=.05:0.1:m;sqrtx=ones(4,1)*sqrt(pi/2./x);figure(2);subplot(2,1,1);y=sqrtx'.*besselj(k+1/2,x');plot(x,y);grid

on;

axis([0,m,-.3,1]);title('0-3阶球Bessel函数');subplot(2,1,2);y=sqrtx'.*bessely(k+1/2,x');plot(x,y);grid

on;

axis([0,m,-1,.5]);title('0-3阶球

ann函数');figure(1);subplot(2,1,1);y=besselj(k,x');plot(x,y);grid

on;

axis([0,m,-.5,1]);title('0-3阶Bessel函数');subplot(2,1,2);y=bessely(k,x');plot(x,y);grid

on;

axis([0,m,-1,.6]);title('0-3阶

ann函数');ex106%ex311,勒让德函数clear;

L=5;

theta=pi*(0:1/100:1);fork=0:L;P=legendre(L,cos(theta));A=sqrt(prod(1:L-k)*(2*L+1)/(4*pi*prod(L+k)))Y(k+1,:)=A*P(k+1,:)';figure(1);

subplot(3,2,k+1);

plot(theta/pi,Y(k+1,:));xlabel('\theta

/

pi');

box

on; grid

on;P=legendre(k+1,cos(theta));P0(k+1,:)=P(1,:)';figure(3);

subplot(3,2,k+1);

plot(theta/pi,P0(k+1,:))xlabel('\theta

/

pi');

box

on; gridon;end;Y50(cosθ)Y51(cosθ)Y53(cosθ)Y55(cosθ)Y52(cosθ)