本征函數(shù)系與本征振動

第四章本征函數(shù)系與本征振動一維本征值問題；正余弦函數(shù)系二維本征值問題；貝塞爾函數(shù)系三維本征值問題；球諧函數(shù)系(勒讓德函數(shù)系)用偏微分方程工具箱(pdetool)求解本征值問題；4.0 PDE分類2

2Laplacian

算子：

x2

y2拉（elliptic）：斯方程

a2u

f熱傳導(dǎo)方程(parabolic)：tu

a2u

f波動方程(hyperbolic):a

u

f22t2u

初邊值問題邊值問題設(shè)函數(shù)u(x,t)

具有變量分離形式，即它可表示為：或則上述三個方程可以寫為：4.1

一維本征值問題對于齊次問題，f=0，在一(二)維情況下，三種方程分別變?yōu)椋憾S拉

斯方程

uxx

uyy

0 (0

x

a;0

y

b)一維熱傳導(dǎo)方程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)這些方程的分離變量中都會出現(xiàn)方程：X

X

0該方程與各種邊值條件結(jié)合可以組成各種本征函數(shù)系二維拉斯方程

Y

Y

0;X

X

0;一維熱傳導(dǎo)方程2

Ta2ta

T

0

T

e

X

X

0一維波動方程02T

a

T

0

T

cos

at

X

X

04.1.1

四種常見的正、余弦本征函數(shù)系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

本征函數(shù)系的圖像及其運動上節(jié)中的四個本征函數(shù)的圖像如下頁圖。顯然，端點為第一類邊界條件時，該端點為波節(jié)，端點為第二類邊界條件時，該端為波腹。程序如下：%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);實際上對于波動方程和熱傳導(dǎo)方程，還要考慮時間本征函數(shù)：對于波動方程，時間函數(shù)為cos(nπa

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

t/l)對于熱傳導(dǎo)方程，時間函數(shù)為exp(-λa2t)

(文中該處p66有誤)考慮時間因子后的動畫分別如程序ex302和ex303所示。顯然，對于波動方程，將形成駐波，固定端(第一類邊界條件)為波節(jié)，開口端(第二類邊界條件)為波腹。對于熱傳導(dǎo)方程，若兩端溫度相同(第一類邊界條件)，對于任何初始溫度分布，最終都將趨于平衡態(tài)，%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;可通過分離變量求得本征函數(shù)：y

Y

y

sin

mπy

/

b

0

Y

b

0YY

k

2

y二維本征值問題矩形區(qū)域的本征模與本征振動具有固定邊界的矩形區(qū)域拉 斯方程的本征值問題：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

相應(yīng)的動畫演示程序為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階貝塞爾函數(shù)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

園形區(qū)域的本征模與本征振動具有固定邊界的圓形區(qū)域拉 斯方程的本征值問題：可通過分離變量求得本征函數(shù)：令u,

R

ann)函數(shù)%ex106(p32-33)

;%貝塞耳(Bessel)函數(shù)\諾伊曼(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函數(shù)');subplot(2,1,2);y=bessely(k,x);plot(x,y);grid

on;

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

ann函數(shù)');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三維本征值問題(球貝塞耳函數(shù)和勒讓德函數(shù))球坐標系中的拉 斯方程在球坐標系中拉 斯方程的本征值問題：可通過分離變量求得本征函數(shù)：

令ur,,

其中稱為l

階連帶(或締合)勒讓德函數(shù)，在中，用legendre(l,x)表示。它們的圖像分別如下幾頁所示。jl(r)與nl(r)分別稱為l階球貝塞耳函數(shù)和l

階球諾伊曼函數(shù)；它們與半數(shù)階貝塞耳函數(shù)Jl+1/2(r)和半數(shù)階諾伊曼函數(shù)Yl+1/2(r)的關(guān)系為：

immlmmlPYcos

el

m

!

4πl(wèi)

m

!2l

1

,

1mlP

cos

l而Y

m(θ,φ)稱為球諧函數(shù)，可以寫為：

nl

π

πY

r2r

l

1/

2J2r

l

1/

2r

;

r

jl

r

在 中，l

階球貝塞耳和

l

階球諾伊曼函數(shù)分別用：sqrt(pi/2/r)*besselj(l+1/2,r)

與sqrt(pi/2/r)*

bessely(l+1/2,r)表示。%ex106(p32-33);Bessel\諾伊曼函數(shù)與球Bessel\諾伊曼函數(shù)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函數(shù)');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函數(shù)');figure(1);subplot(2,1,1);y=besselj(k,x');plot(x,y);grid

on;

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

on;

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

ann函數(shù)');ex106%ex311,勒讓德函數(shù)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θ)