利用傅里叶变换可以得到： 2 2 U tt a ik U U k , t Ak cosakt Bk sin akt U k , t 0 k ; U t k , t 0 k ; 再用初始条件得： U k , t k cosakt k sin akt / ak 通过反变换可得 1 1 x at u x, t x at x at x at d 2 2a 达朗贝尔公式：
x at 0 0; 1 x at d x at / 2a; 0 x at 1 2a 1 / 2a; x at 1
%ex602; (p159) 无限长弦波动的解析解(初位移为0, 初速不为0) clear; M=100; N=80; a=1.0; L=10; T1=8; dx=L/M; dt=T1/N; x=-L:dx:L; t=0:dt:T1;[X,T]=meshgrid(x,t); xp=X+a*T; xp(find(xp<=0))=0; xp(find(xp>=1))=1; xm=X-a*T; xm(find(xm<=0))=0; xm(find(xm>=1))=1; S=(xp-xm)/(2*a); figure(1); h=plot(x,S(1,:),'linewidth',3); axis([-L L 0 .6]); set(h,'erasemode','xor'); for k=2:N+1; pause(0.01); set(h,'ydata',S(k,:)); drawnow; end;
2l Bn 2 2 cos3nπ / 7 cos4nπ / 7 nπa An 0;
%ex603; (p161) % 两端固定的弦波动的解析解 clear; N=100; M=500; K=100;a=1.0;L=1;T=.4; dx=L/M; dt=T/K; x=dx*(0:M); t=dt*(0:K); [X,T]=meshgrid(x,t); w=0; for n=1:N; p=(7+n)*pi/7; q=(7-n)*pi/7+eps; r=n*pi/L; s=n*pi; A(n)=(sin(4*q)-sin(3*q))/q/7-(sin(4*p)-sin(3*p))/p/7; B(n)=0; %B(n)=2*a*L/(s*a)^2*(cos(3*s/7)-cos(4*s/7));A(n)=0; w=w+(A(n)*cos(r*a*T)+B(n)*sin(r*a*T)).*sin(r*X); end; ymax=1.1*max(max(abs(w))); figure(1);h=plot(x,w(1,:),'linewidth',3); axis([0 M*dx -ymax ymax]); set(h,'erasemode','xor'); for k=2:K; pause(0.02); set(h,'ydata',w(k,:)); drawnow; end; figure(2);mesh(X,T,w);

%ex6012; (p157-161) 无穷长弦波动的差分解 clear; II=500; N=240; L=10; T=4; a=1.0; A=1; dx=L/II; dt=T/N; x=dx*(0:II); t=dt*(0:N); K=1:II+1;fai(K)=0; psi(K)=0; C=(a*dt/dx)^2; I=2:II; %fai=A*sin(7*pi*x/L);fai(find(x>4*L/7|x<3*L/7))=0;% (初位移) psi(find(x<4*L/7&x>3*L/7))=1; %(初位移不为0，初速为0)； u(1,:)=fai; u(2,:)=fai+dt*psi; figure(1); h=plot(x,u(1,:),'linewidth',3); %画动画; axis([0,II*dx,-A,A]); set(h,'erasemode','xor'); set(h,'ydata',u(2,:)); drawnow; pause(0.01); for k=2:N; u(k+1,1)=0; u(k+1,I+1)=0; u(k+1,I)=2*u(k,I)-u(k-1,I)+C*(u(k,I+1)-2*u(k,I)+u(k,I-1)); set(h,'ydata',u(k+1,:)); drawnow; pause(0.01); end; figure(2); mesh(x,t,u);






%ex6081; (p171)两端固定弦振动的差分解(有阻尼) clear;II=500; N=300; L=10; T=4; a=1.0; A=1; lamda=5; dx=L/II; dt=T/N; x=dx*(0:II); t=dt*(0:N); K=1:II+1;fai(K)=0; psi(K)=0; C=(a*dt/dx)^2; I=2:II; fai=A*sin(7*pi*x/L);fai(find(x>4*L/7|x<3*L/7))=0; %(初位移) u(1,:)=fai; u(2,I)=u(1,I)+0.5*C*(u(1,I+1)-2*u(1,I)+u(1,I-1)); figure(1); h=plot(x,u(1,:),'linewidth',3); %画动画; axis([0,L,-1.1*A,1.1*A]); set(h,'erasemode','xor'); set(h,'ydata',u(2,:)); drawnow; pause(0.01); for n=2:N; u(n+1,1)=0; u(n+1,II+1)=0; u(n+1,I)=2*u(n,I)-u(n-1,I)+C*(u(n,I+1)-2*u(n,I)+u(n,I-1)); u(n+1,I)=(u(n+1,I)+lamda*dt*u(n,I))/(1+lamda*dt); set(h,'ydata',u(n+1,:)); drawnow; pause(0.01); end; figure(2); mesh(x,t,u);title('lamda')
sin 7x / l ; 3l / 7 x 4l / 7 1. 取 l=1, a=1， x 0; x 其他; 0
1 sin 47 nπ / 7 sin 37 n π / l An 7 n π 1 sin 47 nπ / 7 sin 37 nπ / l ; n 7 7 n π A7 1 / 7; Bn 0; 1; 3l / 7 x 4l / 7 x 0; x 2. 取 l=1, a=1， 其他; 0
A sin 7x / l ; 3l / 7 x 4l / 7 x 1.若初始条件为： 0; 其他 x 0 1 则解析解为： u x, t x at x at
一. 解析解
2
%ex601; (p157) 一维无限长波动的解析解(初速为0)； clear; N=140; M=60; L=10; a=1; A=1; x1= 3*L/7;x2=4*L/7; x=L*(0:N)/N; dt=0.01; t=dt*(0:M); u0=A*sin(7*pi*x/L);u0(find(x<x1))=0;u0(find(x>x2))=0; figure(1); h=plot(x,u0,'linewidth',3);axis([0,L,-A,A]); set(h,'erasemode','xor'); b=7*pi/L; for k=2:length(t); xl=x+a*t(k); xr=x-a*t(k); ul=A*sin(b*xl); ul(find(xl<x1))=0; ul(find(xl>x2))=0; ur=A*sin(b*xr);ur(find(xr<x1))=0;ur(find(xr>x2))=0; set(h,'ydata',(ul+ur)/2 ); drawnow; pause(0.1); end;
差分解程序ex6012 (无穷长，有限长)
差分解程序ex6012 (无穷长，有限长)
7.1.2~3 两端固定的弦的自由振动
两端固定的弦的自由振动的定解问题为：
utt a 2u xx (0 x l ) u 0, t u l , t 0 u x, t 0 x ; u x, t 0 x ; t
2. 初始条件为： x 0
1; 0 x 1 x 0; 其他
则解析解为： 1 x at 1 x at 1 x at u x, t d d d 2a x at 2a 2a x at 0 0; 1 x at 其中： d x at / 2a; 0 x at 1 2a 1 / 2a; x at 1
解析解 程序ex603
解析解 程序ex603
7.1.4 两端固定的弦的振动(有阻尼)
当存在阻尼时，弦振动的振幅会不断减小，定解问题为：
utt a 2u xx 2ut ; (0 x l ) u 0, t 0; u l , t 0; sin 7 πx / l ; 3l/ 7 x 4l/ 7 u x, t 0 x 0; (其他) ut x, t 0 x 0; 2 2 n 1 a t n n 1 n n n n 1 n 1 u 2 u u u 2 u u t u u i i i i 1 i i 1 i i 2 x 1 2 0 u x ; u u i i i 2 t xi 0; i