齐次弦振动方程的MATLAB解法
齐次弦振动方程的MATLAB 解法
【摘要】
弦振动问题是一个典型的波动方程的建立与求解问题。本文通过利用MATLAB 特有的方程求解与画图功能，有效地构造和求解了齐次弦振动方程。并通过图像，可以直观感受方程的解，从而加深对这一问题物理意义的理解。
【关键词】
振动方程 MATLAB 求解 数学物理方法
【正文】
在细弦上任意取微元分析其受力情况，通过Newton 定律建立细弦振动的运动方程，可以求得弦振动的泛定方程为2tt xx u a u =。
要得出振动方程的解，除了泛定方程外，我们还需要知道具体问题的初始条件与边界条件。在弦振动问题里，初始条件可以从初始位移和初始速度考虑，即：
00
()()t t
t u
x u x ?ψ==?=??=??
边界条件是描述物理问题在边界上受约束的状态，在弦振动方程里可以归结为三类边界问题：
（1） 第一类边界问题：0
0,0,x x L
u u
==== 称为固定端。
（2） 第二类边界问题：0()
0,,x x x L F t u
u SY
====
特别的，若()0F t =，0,x
x L
u ==称x L =为自由端。
（3） 第三类边界问题：第一类和第二类边界问题的线性组合。
一、 两端固定的弦振动问题
两端固定的弦振动方程的定解问题可表示如下：
2000,0,0
0,0(),()
tt xx x x L t t t u a u x L t u u u
x u x ?ψ====?=<<>?
==??==? 1、初始位移不为0，初始速度为0 不妨设：
734sin
,()()770,()l l x x x l otherwise π
??≤≤?=???
，()0x ψ=
（1）特征函数求解解 由d ’Alembert 公式：
1
(,)[()()]2
u x t x at x at ??=++-
从而我们可以得到方程的级数解：
n A =
()1
43[sin(7)sin(7)](7)77
1
43
[sin(7)sin(7)],7(7)77
n n n n n n n π
πππ
ππ------=-≠-
而我们知道，弦振动的泛定方程属于本征问题：
''00,0x x L X X X X λ==+=?
?
==? 它在两个边界上都有第一类其次边界条件，它的本征值与本征函数为：
222,sin ,1,2,3n n x n l l
ππλ==
将系数带入方程，级数中每一项都是一个驻波，定义子程序wfun.m 计算不同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
22
t c a x ?=? 于是，初始条件可以表示为：
,1,21,11,1,11
[()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
设初始速度为：
34
1,()
(()77
0,
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: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))...
/(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)
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
function wtx=psi1fun1(N,t)
x=0:0.001: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）
clear
N=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);
for k=2:N
set(h,'XData',x,'YData',u(:,2));
drawnow;
pause(0.01)
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