齐次弦振动方程的解法

齐次弦振动方程的MATLAB解法
【摘要】
弦振动问题是一个典型的波动方程的建立与求解问题。

本文通过利用MATLAB特有的方程求解与画图功能，有效地构造和求解了齐次弦振动方程。

并通过图像，可以直观感受方程的解，从而加深对这一问题物理意义的理解。

【关键词】
振动方程 MATLAB求解数学物理方法
【正文】
在细弦上任意取微元分析其受力情况，通过Newton定律建立细弦振
要得出振动方程的解，除了泛定方程外，我们还需要知道具体问题的初始条件与边界条件。

在弦振动问题里，初始条件可以从初始位移和初始速度考虑，即：
边界条件是描述物理问题在边界上受约束的状态，在弦振动方程里可以归结为三类边界问题：
（1） 称为固定端。

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

一、 两端固定的弦振动问题
两端固定的弦振动方程的定解问题可表示如下：
1、初始位移不为0，初始速度为0
不妨设：
（1）特征函数求解解
由d ’Alembert 公式：
从而我们可以得到方程的级数解：
而我们知道，弦振动的泛定方程属于本征问题：
它在两个边界上都有第一类其次边界条件，它的本征值与本征函数为：
将系数带入方程，级数中每一项都是一个驻波，定义子程序计算不同n 的求和各项，再用主程序jxj 将它们加起来，得到动画图形。

（MATLAB 代码见附录1（1））
（2）差分方程求解
利用差分方程同样可以求出问题的解。

令,
x i x t j t
=∆=
∆，将微分方程改写成差分方程，即有
,11,1,,,1
()2(1)
i l i l i l i j i j
u c u u c u u
++--
=++--
其中，
()
()
2
2
2
t
c a
x
∆
=
∆
于是，初始条件可以表示为：
,1
,21,11,1,1
1
[()2(1)]
2
i
i i i i
u
u c u u c u
ϕ
+-
=
=++-
作图时，先画出()x
ϕ的图形，然后再用x at
+或x at
-代替其中的x，改变at的值，就画出了不同时刻()
x at
ϕ+，()
x at
ϕ-的图形。

（MATLAB代码见附录1（2））
解得的动态图形如下：
2、初始位移为0，初始速度不为0
设初始速度为： 341,()(()770,l l x x otherwise
ψ⎧≤≤⎪=⎨⎪⎩
（1）特征函数求解
通过求本征函数与本征值的方法我们可以得到方程的解析解：
1
sin sin
n
n
n a n
u B t x
l l
ππ
∞
=
=∑
其中系数，22
234
(cos cos)
77
n
l
B na na
n a
π
=
-，类似的，用函数计算级数中的各项，再在主函数中调用便可得解。

（MATLAB代码见附录2（1））
（2）差分方程求解
类似于问题1，我们还可以采用差分方程求解，不过需要注意的是，
题目中的初始条件应表示为：
,2
i
u t
ψ
=∆。

（MATLAB代码见附录2（2））
解得的动画图形如下：
【总结】
通过运用MATLAB构造和求解齐次弦振动方程，绘制了相关图像，直观感受了方程解，加深了对其物理意义的理解。

借助于计算机来做计算和
研究的过程涉及到建立模型，选择方法，语言编程和结果分析。

通过此次问题的探究，培养和训练了自学能力和操作能力，获益匪浅。

【参考文献】
1、李明奇田太心《数学物理方程》电子科技大学出版社 2010
2、彭芳麟《数学物理方程的MATLAB解法与可视化》清华大学出版社
2004
3、彭芳麟《计算物理基础》高等教育出版社 2010
4、谢进李大美《MATLAB与计算方法实验》武汉大学出版社 2009【附录】
附录1
（1）
function jxj
N=50
t=0::;
x=0::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::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))...
/(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=;
dt=; 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
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::; x=0::1;
ww=psi1fun1(N,0);
h=plot(x,ww,'linewidth',2);
axis([0,1,,]);
sy=[];
for n=2:length(t)
ww=psi1fun1(N,t(n));
set(h,'ydata',ww);
drawnow;
pause
sy=[sy,sum(ww)];
end
function wtx=psi1fun1(N,t)
x=0::1; a=1; wtx=0;
for k=1:N
Bk=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）
N=4025; dx=;
dt=; c=dt*dt/dx/dx;
x=linspace(0,1,420);
u(1:420,1)=0;
u(180:240,2)=dt*;
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
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);。