(MATLAB代码见附录2(2)) 解得的动画图形如下:
6
【总结】
通过运用MATLAB构造和求解齐次弦振动方程,绘制了相关图像,直观感受了方程解,加深了对其物理意义的理解。借助于计算机来做计算和研究的过程涉及到建立模型,选择方法,语言编程和结果分析。通过此次问题的探究,培养和训练了自学能力和操作能力,获益匪浅。
【参考文献】
1、 李明奇 田太心 《数学物理方程》 电子科技大学出版社
2010
2、 彭芳麟 《数学物理方程的MATLAB解法与可视化》 清华大
学出版社 2004
3、 彭芳麟 《计算物理基础》 高等教育出版社 2010 4、 谢进 李大美 《MATLAB与计算方法实验》 武汉大学出
版社 2009 【附录】 附录1 (1)
function jxj N=50
7
t=0:0.005:2.0; x=0:0.001:1; ww=wfun(N,0); ymax=max(abs(ww)); h=plot(x,ww);
axis([0,1,-ymax,ymax]) sy=[];
for n=2:length(t) ww=wfun(N,t(n)); set(h,'ydata',ww); drawnow;
sy=[sy,sum(ww)]; end
function wtx=wfun(N,t) x=0:0.001:1; a=1; wtx=0; for I=1:N if I~=7
wtx=wtx+((sin(pi*(7-I)*4/7)-sin(pi*(7-I)*3/7))...
/(7-I)/pi-(sin(pi*(7+I)*4/7)-sin(pi*(7+I)*3/7))...
8
/(7+I)/pi)*cos(I*pi*a*t).*sin(I*pi*x); else
wtx=wtx+1/7*cos(I*pi*a*t).*sin(I*pi*x); end end (2)
N=4010; dx=0.0024; dt=0.0005; c=dt*dt/dx/dx; x=linspace(0,1,420); u(1:420,1)=0;
u(181:240,1)=sin(pi*x(181:240)*7);
u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1));
h=plot(x,u(:,1),'linewidth',2); axis([0,1,-1,1]);
set(h,'EraseMode','xor','MarkerSize',18); for k=2:N
set(h,'XData',x,'YData',u(:,2)); drawnow; pause(0.1)
9
u(2:419,3)=2*u(2:419,2)-u(2:419,1)+c*(u(3:420,2)...
-2*u(2:419,2)+u(1:418,2)); u(2:419,1)=u(2:419,2); u(2:419,2)=u(2:419,3); end 2 1)
function psi N=50;
t=0:0.005:2.0; x=0:0.001:1; ww=psi1fun1(N,0);
h=plot(x,ww,'linewidth',2); axis([0,1,-0.1,0.1]); sy=[];
for n=2:length(t) ww=psi1fun1(N,t(n)); set(h,'ydata',ww); drawnow; pause(1.5) sy=[sy,sum(ww)]; end
10
附录(