齐次弦振动方程的MATLAB解法

1齐次弦振动方程的MATLAB解法【摘要】弦振动问题是一个典型的波动方程的建立与求解问题。本文通过利用MATLAB特有的方程求解与画图功能，有效地构造和求解了齐次弦振动方程。并通过图像，可以直观感受方程的解，从而加深对这一问题物理意义的理解。【关键词】振动方程MATLAB求解数学物理方法【正文】在细弦上任意取微元分析其受力情况，通过Newton定律建立细弦振动的运动方程，可以求得弦振动的泛定方程为2ttxxuau。要得出振动方程的解，除了泛定方程外，我们还需要知道具体问题的初始条件与边界条件。在弦振动问题里，初始条件可以从初始位移和初始速度考虑，即：00()()tttuxux边界条件是描述物理问题在边界上受约束的状态，在弦振动方程里可以归结为三类边界问题：2（1）第一类边界问题：00,0,xxLuu称为固定端。（2）第二类边界问题：0()0,,xxxLFtuuSY特别的，若()0Ft，0,xxLu称xL为自由端。（3）第三类边界问题：第一类和第二类边界问题的线性组合。一、两端固定的弦振动问题两端固定的弦振动方程的定解问题可表示如下：2000,0,00,0(),()ttxxxxLtttuauxLtuuuxux1、初始位移不为0，初始速度为0不妨设：734sin,()()770,()llxxxlotherwise，()0x（1）特征函数求解解由d’Alembert公式：1(,)[()()]2uxtxatxat从而我们可以得到方程的级数解：nA143[sin(7)sin(7)](7)77143[sin(7)sin(7)],7(7)77nnnnnnn3而我们知道，弦振动的泛定方程属于本征问题：0''00,0xxLXXXX它在两个边界上都有第一类其次边界条件，它的本征值与本征函数为：222,sin,1,2,3nnxnll将系数带入方程，级数中每一项都是一个驻波，定义子程序wfun.m计算不同n的求和各项，再用主程序jxj将它们加起来，得到动画图形。（MATLAB代码见附录1（1））（2）差分方程求解利用差分方程同样可以求出问题的解。令,xixtjt，将微分方程改写成差分方程，即有,11,1,,,1()2(1)ilililijijucuucuu其中，222tcax于是，初始条件可以表示为：,1,21,11,1,11[()2(1)]2iiiiiuucuucu作图时，先画出x的图形，然后再用xat或xat代替其中的x，改变at的值，就画出了不同时刻()xat，()xat4的图形。（MATLAB代码见附录1（2））解得的动态图形如下：52、初始位移为0，初始速度不为0设初始速度为：341,()(()770,llxxotherwise（1）特征函数求解通过求本征函数与本征值的方法我们可以得到方程的解析解：1sinsinnnnanuBtxll其中系数，22234(coscos)77nlBnanana，类似的，用函数计算级数中的各项，再在主函数中调用便可得解。（MATLAB代码见附录2（1））（2）差分方程求解类似于问题1，我们还可以采用差分方程求解，不过需要注意的是，题目中的初始条件应表示为：,2iut。6（MATLAB代码见附录2（2））解得的动画图形如下：7【总结】通过运用MATLAB构造和求解齐次弦振动方程，绘制了相关图像，直观感受了方程解，加深了对其物理意义的理解。借助于计算机来做计算和研究的过程涉及到建立模型，选择方法，语言编程和结果分析。通过此次问题的探究，培养和训练了自学能力和操作能力，获益匪浅。【参考文献】1、李明奇田太心《数学物理方程》电子科技大学出版社20102、彭芳麟《数学物理方程的MATLAB解法与可视化》清华大学出版社20043、彭芳麟《计算物理基础》高等教育出版社20104、谢进李大美《MATLAB与计算方法实验》武汉大学出版社2009【附录】附录1（1）functionjxjN=508t=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=[];forn=2:length(t)ww=wfun(N,t(n));set(h,'ydata',ww);drawnow;sy=[sy,sum(ww)];endfunctionwtx=wfun(N,t)x=0:0.001:1;a=1;wtx=0;forI=1:NifI~=7wtx=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))...9/(7+I)/pi)*cos(I*pi*a*t).*sin(I*pi*x);elsewtx=wtx+1/7*cos(I*pi*a*t).*sin(I*pi*x);endend（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);fork=2:Nset(h,'XData',x,'YData',u(:,2));drawnow;pause(0.1)10u(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）functionpsiN=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=[];forn=2:length(t)ww=psi1fun1(N,t(n));set(h,'ydata',ww);drawnow;pause(1.5)sy=[sy,sum(ww)];end11functionwtx=psi1fun1(N,t)x=0:0.001:1;a=1;wtx=0;fork=1:NBk=2/(k*k*pi*pi)*(cos(3*k*pi/7)-cos(4*k*pi/7));wtx=wtx+Bk*sin(k*pi*t)*sin(k*pi*x);end（2）clearN=4025;dx=0.0024;dt=0.0005;c=dt*dt/dx/dx;x=linspace(0,1,420);u(1:420,1)=0;u(180:240,2)=dt*0.5;h=plot(x,u(:,1),'linewidth',2);axis([0,1,-1,1]);set(h,'EraseMode','xor','MarkerSize',18);fork=2:Nset(h,'XData',x,'YData',u(:,2));drawnow;pause(0.01)12u(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