齐次弦振动方程的MATLAB解法

齐次弦振动方程的MATLAB解法

【摘要】

弦振动问题是一个典型的波动方程的建立与求解问题。本文通过利用MATLAB特有的方程求解与画图功能，有效地构造和求解了齐次弦振动方程。并通过图像，可以直观感受方程的解，从而加深对这一问题物理意义的理解。

【关键词】

振动方程 MATLAB求解 数学物理方法

【正文】

在细弦上任意取微元分析其受力情况，通过Newton定律建立细弦振动的运动方程，可以求得弦振动的泛定方程为utt?a2uxx。

要得出振动方程的解，除了泛定方程外，我们还需要知道具体问题的初始条件与边界条件。在弦振动问题里，初始条件可以从初始位移和初始速度考虑，即：

??u???utt?0t?0??(x)??(x)

边界条件是描述物理问题在边界上受约束的状态，在弦振动方程里可以归结为三类边界问题：

1

（1） 第一类边界问题：u（2） 第二类边界问题：uuxx?Lx?0?0,ux?L?0, 称为固定端。

F(t),特别的，若F(t)?0，SYx?0?0,uxx?L??0,称x?L为自由端。

（3） 第三类边界问题：第一类和第二类边界问题的线性组合。

一、 两端固定的弦振动问题

两端固定的弦振动方程的定解问题可表示如下：

??u2tt?auxx,0?x?L,t?0?ux?0?0,ux?L?0?x),u ?ut?0??(tt?0??(x)1、初始位移不为0，初始速度为0 不妨设：

?7?3l4l?(x)???sin?lx,(7?x?7)，?(x)??0,(otherwise)0

（1）特征函数求解解 由d’Alembert公式：

u(x,t)?12[?(x?at)??(x?at)]

从而我们可以得到方程的级数解：

A4n?1(7?n)?[sin(7?n)?7?sin(7?n)?37]?14

(7?n)?[sin(7?n)?7?sin(7?n)?37],?n?7?

2

而我们知道，弦振动的泛定方程属于本征问题：

X''??X?0? ?X?0,X?0x?L?x?0它在两个边界上都有第一类其次边界条件，它的本征值与本征函数为：

n2?2n???2,sinx,n?1,2,3?

ll将系数带入方程，级数中每一项都是一个驻波，定义子程序wfun.m计算不同n的求和各项，再用主程序jxj将它们加起来，得到动画图形。

（MATLAB代码见附录1（1））

（2）差分方程求解

利用差分方程同样可以求出问题的解。令x?i?x,t?j?t，将微分方程改写成差分方程，即有

ui,l?1?c(ui?1,l?ui?1,l)?2(1?c)ui,j?ui,j?1

其中，

2??t? c?a2??x?2于是，初始条件可以表示为：

ui,1?? 1ui,2?[c(ui?1,1?ui?1,1)?2(1?c)ui,1]2作图时，先画出??x?的图形，然后再用x?at或x?at代替其中的x，改变at的值，就画出了不同时刻?(x?at)，?(x?at)3

的图形。

（MATLAB代码见附录1（2））

解得的动态图形如下：

4

2、初始位移为0，初始速度不为0 设初始速度为：

4l?3l?1,(?x?)?((x)??77

??0,otherwise（1） 特征函数求解

通过求本征函数与本征值的方法我们可以得到方程的解析解：

n?an?u??Bnsintsinx

lln?1?其中系数，Bn?34(cosna?cosna)，类似的，用函数计算级22n?a772l数中的各项，再在主函数中调用便可得解。 （MATLAB代码见附录2（1））

（2） 差分方程求解

类似于问题1，我们还可以采用差分方程求解，不过需要注意的是，题目中的初始条件应表示为：ui,2???t。

5

（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

附录（

本文来源：https://www.bwwdw.com/article/iec2.html