数学物理建模与计算机辅助设计数学物理建模与计算机辅助设计上机练习专题上机练习专题6 6：数学物理方程的计：数学物理方程的计 算机求解和可视化算机求解和可视化数学物理建模与计算机辅助设计Page 2一维一维无界无界弦自由振动（弦自由振动（无强迫力无强迫力）齐次齐次定解问题为：定解问题为：直接利用直接利用达朗贝尔公式达朗贝尔公式求解：求解：（ ） 利用利用行波法行波法求解数学物理定解问题求解数学物理定解问题数学物理建模与计算机辅助设计达朗贝尔解的物理意义达朗贝尔解的物理意义：正向波（或右行波）：正向波（或右行波）：反向波（或左行波）：反向波（或左行波）传播速度为传播速度为 a a数学物理建模与计算机辅助设计实例实例1 1：一维：一维无界弦无界弦振动定解问题振动定解问题: :根据达朗贝尔公式，位移为根据达朗贝尔公式，位移为 ：数学物理建模与计算机辅助设计clearclear u0=2; u(1:560)=0;u0=2; u(1:560)=0; x1=-1;x2=1;x1=-1;x2=1; x=linspace(-2,2,560);x=linspace(-2,2,560); u(141:280)=2*u0*(x(141:280)-x1)/(x2-x1);u(141:280)=2*u0*(x(141:280)-x1)/(x2-x1); u(281:420)=2*u0*(x2-x(281:420)/(x2-x1);u(281:420)=2*u0*(x2-x(281:420)/(x2-x1); uu=u;uu=u; h=plot(x,u,linewidth,3);h=plot(x,u,linewidth,3); axis(2*x1,2*x2,-4,4);axis(2*x1,2*x2,-4,4); set(h,EraseMode,xor)set(h,EraseMode,xor) for at=2:140for at=2:140lu(1:560)=0; ru(1:560)=0; lu(1:560)=0; ru(1:560)=0;lx=141:420-at; rx=141:420+at; lx=141:420-at; rx=141:420+at;lu(lx)=0.5*uu(141:420); ru(rx)=0.5*uu(141:420); lu(lx)=0.5*uu(141:420); ru(rx)=0.5*uu(141:420);u=lu+ru; u=lu+ru;set(h,Xdata,x,YData,u); set(h,Xdata,x,YData,u);drawnow; drawnow;pause(0.1) pause(0.1) endend数学物理建模与计算机辅助设计某时刻无界弦振动图型某时刻无界弦振动图型u u( (x x, ,t t) )x x数学物理建模与计算机辅助设计实例实例2 2：一维：一维无界弦无界弦振动定解问题振动定解问题: :由由达朗达朗贝贝贝贝尔尔公式公式得得: :数学物理建模与计算机辅助设计其中其中: :数学物理建模与计算机辅助设计clearclear x1=0;x2=1;t=0:0.005:8;x=-10:0.1:10;a=1;x1=0;x2=1;t=0:0.005:8;x=-10:0.1:10;a=1; X,T=meshgrid(x,t); X,T=meshgrid(x,t); xpat=X+a*T;xpat=X+a*T; xpat(find(xpat=x2)=1; xpat(find(xpat=x2)=1; xmat=X-a*T;xmat=X-a*T; xmat(find(xmat=x2)=1; xmat(find(xmat=x2)=1; jf1=1/2/a*(xpat-x1);jf2=1/2/a*(xmat-x1);jf3=1/2/a*(xpat-x1)-(xmat-x1);jf1=1/2/a*(xpat-x1);jf2=1/2/a*(xmat-x1);jf3=1/2/a*(xpat-x1)-(xmat-x1); subplot(3,1,1)subplot(3,1,1) h1=plot(x,jf1(1,:),linewidth,3);h1=plot(x,jf1(1,:),linewidth,3); set(h1,erasemode,xor);set(h1,erasemode,xor); axis(-10 10 -1 1)axis(-10 10 -1 1) hold onhold on subplot(3,1,2)subplot(3,1,2) h2=plot(x,jf2(1,:),linewidth,3);h2=plot(x,jf2(1,:),linewidth,3); set(h2,erasemode,xor);set(h2,erasemode,xor); axis(-10 10 -1 1)axis(-10 10 -1 1) hold onhold on subplot(3,1,3)subplot(3,1,3) h3=plot(x,jf3(1,:),linewidth,3);h3=plot(x,jf3(1,:),linewidth,3); set(h3,erasemode,xor);set(h3,erasemode,xor); axis(-10 10 -1 1)axis(-10 10 -1 1) hold onhold on for j=2:length(t)for j=2:length(t)pause(0.01) pause(0.01)set(h1,ydata,jf1(j,:); set(h2,ydata,jf2(j,:); set(h3,ydata,jf3(j,:); set(h1,ydata,jf1(j,:); set(h2,ydata,jf2(j,:); set(h3,ydata,jf3(j,:);drawnow; drawnow; endend数学物理建模与计算机辅助设计左行波左行波: : 右行波右行波: : 叠加后叠加后 的波形的波形: : 数学物理建模与计算机辅助设计利用利用分离变量法分离变量法求解数学物理定解问题求解数学物理定解问题长为长为 l l ，两端固定的均匀弦的自，两端固定的均匀弦的自由振动定解问题：由振动定解问题：通解：通解：数学物理建模与计算机辅助设计特解特解改写为：改写为： 驻波驻波叠加叠加角频率角频率初位相初位相数学物理建模与计算机辅助设计振幅振幅: : 波节波节: : 波腹波腹: :数学物理建模与计算机辅助设计function zbfunction zb t=0:0.005:2.0; x=0:0.001:1;t=0:0.005:2.0; x=0:0.001:1; ww1=ww1=wfunwfun(1,0);ww2=(1,0);ww2=wfunwfun(2,0);(2,0); ww3=ww3=wfunwfun(3,0);ww4=(3,0);ww4=wfunwfun(4,0);(4,0); ymax1=max(abs(ww1);ymax1=max(abs(ww1); figurefigure subplot(4,1,1)subplot(4,1,1) h1=plot(x,ww1,r,linewidth,5);h1=plot(x,ww1,r,linewidth,5); axis(0,1,-ymax1,ymax1)axis(0,1,-ymax1,ymax1) subplot(4,1,2)subplot(4,1,2) ymax2=max(abs(ww2);ymax2=max(abs(ww2); h2=plot(x,ww2,r,linewidth,5);h2=plot(x,ww2,r,linewidth,5); axis(0,1,-ymax2,ymax2)axis(0,1,-ymax2,ymax2) subplot(4,1,3)subplot(4,1,3) ymax3=max(abs(ww3);ymax3=max(abs(ww3); h3=plot(x,ww3,r,linewidth,5);h3=plot(x,ww3,r,linewidth,5); axis(0,1,-ymax3,ymax3)axis(0,1,-ymax3,ymax3) subplot(4,1,4)subplot(4,1,4) ymax4=max(abs(ww4);ymax4=max(abs(ww4); h4=plot(x,ww4,r,linewidth,5);h4=plot(x,ww4,r,linewidth,5); axis(0,1,-ymax4,ymax4)axis(0,1,-ymax4,ymax4)for n=2:length(t)for n=2:length(t)ww1=wfun(1,t(n); ww1=wfun(1,t(n);set(h1,ydata,ww1); set(h1,ydata,ww1);ww2=wfun(2,t(n); ww2=wfun(2,t(n);set(h2,ydata,ww2); set(h2,ydata,ww2); drawnow; drawnow;ww3=wfun(3,t(n); ww3=wfun(3,t(n);set(h3,ydata,ww3); set(h3,ydata,ww3); ww4=wfun(4,t(n); ww4=wfun(4,t(n);set(h4,ydata,ww4) set(h4,ydata,ww4)drawnow; drawnow; endendfunction wtx= function wtx=wfunwfun(k,t)(k,t)x=0:0.001:1; a=1; x=0:0.001:1; a=1;wtx=cos(k*pi*a*t).*sin(k*pi*x); wtx=cos(k*pi*a*t).*sin(k*pi*x);数学物理建模与计算机辅助设计n=1n=1n=2n=2n=3n=3n=4n=4数学物理建模与计算机辅助设计Page 16对代数表达式进行可视化二维本征值问题二维本征值问题 矩形区域的本征模与本征振动矩形区域的本征模与本征振动边长为b和c的的四周固定的矩形膜振动的本征值问题为采用分离变量法可以得到本征模和本征值为数学物理建模与计算机辅助设计Page 17 求前9个本征值 绘制前4个本征函数的图形b=2; c=1; m,n=meshgrid(1:3); L=(m*pi./c).2+(n*pi./b).2) L =12.3370 41.9458 91.293819.7392 49.3480 98.696032.0762 61.6850 111.0330x=0:0.01:b; y=0:0.01:c; X,Y=meshgrid(x,y); w11=sin(pi*Y./c).*sin(pi*X./b); w12=sin(2*pi*Y./c).*sin(pi*X./b); w21=sin(pi*Y./c).