一類弦振動(dòng)問(wèn)題的計(jì)算機(jī)輔助分析*

孫麗男，張馥菊

(1.黑河學(xué)院;2.大連理工大學(xué);3.哈爾濱師范大學(xué))

0 引言

弦振動(dòng)問(wèn)題是數(shù)學(xué)物理方程中的經(jīng)典問(wèn)題，弦振動(dòng)理論在生產(chǎn)、生活中有著較為廣泛的應(yīng)用.如對(duì)于斜拉橋、高空觀光纜車運(yùn)行安全和游客人身安全而言，檢測(cè)并控制鋼索工作過(guò)程中的張力是至關(guān)重要的，而此類問(wèn)題最終歸結(jié)為弦振動(dòng)問(wèn)題的研究.另外，工廠車間進(jìn)行絡(luò)紗時(shí)，可以觀察到紗會(huì)以相當(dāng)快的速度沿著中心軸進(jìn)行退繞，達(dá)到一定條件時(shí)便會(huì)形成多節(jié)氣圈，這種現(xiàn)象也可以通過(guò)弦振動(dòng)理論加以說(shuō)明.弦振動(dòng)是波動(dòng)的一種特殊形式，因此弦振動(dòng)理論在聲學(xué)中也占有舉足輕重的地位，弦振動(dòng)理論與現(xiàn)代先進(jìn)的軟件和硬件相結(jié)合，在許多領(lǐng)域諸如音樂(lè)物理學(xué)、材料學(xué)和系統(tǒng)分析中也都得到了廣泛的應(yīng)用［1］.弦振動(dòng)問(wèn)題是由弦振動(dòng)方程附加初邊值條件構(gòu)成的，弦振動(dòng)方程是一類重要的偏微分方程，它關(guān)于時(shí)間和空間變量的導(dǎo)數(shù)均是二階的，在解決實(shí)際問(wèn)題的過(guò)程中，遇到的弦振動(dòng)問(wèn)題通常既有初始條件又有邊界條件，是一類混合問(wèn)題.對(duì)于一些簡(jiǎn)單的混合弦振動(dòng)問(wèn)題可以通過(guò)分離變量法、疊加原理、齊次化原理及傅里葉變換等方法求其解析解或者是廣義解，但是對(duì)于一些復(fù)雜的弦振動(dòng)問(wèn)題，傳統(tǒng)的方法求解起來(lái)就比較困難，或者求出的解的形式相對(duì)復(fù)雜，與實(shí)際相結(jié)合較為困難.這時(shí)，就可以通過(guò)數(shù)值方法進(jìn)行求解，借助計(jì)算機(jī)進(jìn)行輔助分析，從而使求解過(guò)程變的輕松，并且能夠?qū)⑺媒Y(jié)果以圖形等形式直觀地顯示出來(lái)［3］.

文中考慮了以下兩端固定的有界弦振動(dòng)問(wèn)題

其中，上述問(wèn)題的解析解可以利用分離變量等方法進(jìn)行求解［2］，但為了更直觀地展示上述問(wèn)題解的形式，以下將采用數(shù)值方法進(jìn)行求解并利用功能強(qiáng)大的Matlab軟件進(jìn)行輔助分析，最后將解析解和數(shù)值解進(jìn)行比較.

1 數(shù)值分析

該節(jié)首先給出求解問(wèn)題(1)的有效的數(shù)值方法，然后利用這些方法及Matlab軟件求出數(shù)值解，并且比較幾種數(shù)值方法的優(yōu)劣.

1.1 數(shù)值方法

求解偏微分方程的數(shù)值方法有很多，比較常見(jiàn)的方法就是有限元法和差分法，文獻(xiàn)［4］中利用有限元法對(duì)(1)進(jìn)行了求解.文中將利用差分法對(duì)(1)進(jìn)行求解，主要采用顯式、隱式兩種差分格式.

1.1.1 區(qū)域剖分

首先，對(duì)定解問(wèn)題(1)的求解區(qū)域W={(x，t)|0＜x＜l，t＞0}進(jìn)行剖分.在區(qū)間［0，1］內(nèi)取n-1個(gè)等分點(diǎn)xi=ih，i=1，…，n-1將區(qū)間n等分，其中x0=0，xn=l，h為空間步長(zhǎng).同樣，在時(shí)間區(qū)間內(nèi)插入若干等分點(diǎn)tj=jτ，j=1，2，…，其中，t0=0，τ為時(shí)間步長(zhǎng).用兩族平行直線x=xi，1＜i＜n;t=tj，j=1，2，… 將區(qū)域W剖分成矩形網(wǎng)格，其中(xi，tj)為結(jié)點(diǎn)，以下將討論問(wèn)題(1)在結(jié)點(diǎn)處的數(shù)值解.

1.1.2 方程離散

用差分方法求解偏微分方程的關(guān)鍵是對(duì)方程的離散，下面采用顯示、隱式差分格式對(duì)問(wèn)題(1)中方程進(jìn)行離散，得到如下形式的差分方程.

(i)顯示差分方程

(i)隱式差分格式

其中i=1，2，…，n-1;j=1，2，…，令則上式可改寫為

1.1.3 初值條件的處理

問(wèn)題(1)對(duì)應(yīng)的初值條件可進(jìn)行如下處理

1.1.4 邊值條件的處理

問(wèn)題(1)對(duì)應(yīng)的邊值可進(jìn)行如下處理

1.2 數(shù)值實(shí)驗(yàn)

下面將利用計(jì)算機(jī)實(shí)驗(yàn)對(duì)問(wèn)題(1)的數(shù)值解進(jìn)行輔助分析，并將其可視化，實(shí)驗(yàn)主要采用的軟件是MATLAB軟件.

首先，根據(jù)顯示差分格式編制MATLAB程序模擬問(wèn)題(1)的數(shù)值解，獲得動(dòng)態(tài)結(jié)果如圖1所示.

圖1

在圖1實(shí)驗(yàn)中取初始位移為

初始速度ψ(x)取為零，時(shí)間步長(zhǎng)τ=0.0005，空間步長(zhǎng)h=0.0024.通過(guò)對(duì)上述結(jié)果進(jìn)行觀察，發(fā)現(xiàn)初始狀態(tài)的波平均分成兩個(gè)波向相反方向傳播，在它們到達(dá)弦的端點(diǎn)之前，和波在無(wú)限長(zhǎng)弦上的傳播是一致的，到達(dá)弦的端點(diǎn)以后，兩個(gè)波同時(shí)發(fā)生反射，出現(xiàn)了半波損失，這是有別于無(wú)限長(zhǎng)弦的振動(dòng)的，接下來(lái)，被反射回來(lái)的兩個(gè)波進(jìn)行相向運(yùn)動(dòng)并在弦的中間相遇，相遇后兩個(gè)波疊加繼續(xù)傳播，整個(gè)波的傳播過(guò)程清晰明了，實(shí)現(xiàn)了問(wèn)題(1)的數(shù)值解的可視化，有利于進(jìn)一步研究更為復(fù)雜的弦振動(dòng)問(wèn)題.

對(duì)于問(wèn)題(1)還可以根據(jù)隱式差分格式編制MATLAB程序模擬它的數(shù)值解，相應(yīng)程序的編碼要比顯示差分格式復(fù)雜的多，涉及到三對(duì)角方程組的求解等技術(shù)性問(wèn)題，具體程序如下

clear;clc;

format short e

a=1;l=1;n=420;=0.0005;

N=input(＇請(qǐng)輸入運(yùn)行次數(shù) N的值:＇);

h = l/n;x = zeros(n+1，1);c =a^2*τ^2/h^2;

for i=1:n

x(i+1)=i*h;

end

u=zeros(n+1，N);

foc i=181:240

u(i，1)=0.05*sin(pi*x(i)*7/l);

u(i，2)=0.05*sin(pi*x(i)*7/l)+1/2*c*0.05*(sin(pi*x(i+1)*7/l)-2*sin(pi*x(i)*7/l) + sin(pi*x(i -1)*7/l));

end

h=plot(x，u(:，1)，＇linewidth＇，3);axis(［0，1， - 0.05，0.05］);set(h，＇EraseMode＇，＇xor，＇MackecSize＇，18)

for i=1:N

set(h，＇XData＇，x，＇YData＇，u(:，1));drawnow;

B=zeros(n-1，1);A=zeros(n-2，1);R=zeros(n-2，1);S=zeros(n-1，1);E=zeros(n-1，1);

F=zeros(n-2，1);G=zeros(n-2，1);I=zeros(n-1，1);H=zeros(［n-1，n-1］);

for i=1:n-2

E(i)=-(1+c);F(i)=1/2*c;G(i)=1/2*c;I(i)=u(i，1);

end

E(n-1)=-(1+c);H=diag(E)+diag(F，1)+diag(G，-1);I(n-1)=u(n-1，1);J=H*I;

foc i=1:n-2

B(i)=1+c;A(i)=-1/2*c;C(i)=-1/2*c;S(i)=2*u(i，2)+J(i);

end

B(n-1)=1+c;S(n-1)=2*u(n-1，2)+J(n-1);u(1，2)=0;u(n+1，2)=0;S(1)=u(1，2)+J(1);

%追趕法

Y=zecos(n-1，1);X=Y;v=zecos(1，n-1);g=v;v(1)=B(1);Y(1)=S(1);A(2:n-1)=A;

foc i=2:n-1

g(i)=A(i)/v(i-1);v(i)=B(i)-g(i)*C(i-1);Y(i)=S(i)-g(i)*Y(i-1);

end

U=zecos(n-1);L=eye(n-1);

foc i=1:n-2

L(i+1，i)=g(i+1);

end

foc i=1:n-2

U(i，i)=v(i);

U(i，i+1)=C(i);

end

U(n-1，n-1)=v(n-1);X(n-1)=Y(n-1)/v(n-1);foc i=n-2:-1:1

X(i)=(Y(i)-C(i)*X(i+1))/v(i);

end

u(1:n-1，3)=X;u(:，1)=u(:，2);u(:，2)=u(:，3);

end

運(yùn)行結(jié)果的動(dòng)態(tài)展示如圖2所示.

圖2

利用隱式差分格式得到的數(shù)值結(jié)果與顯式差分格式總體上是一致的，二者均較好地逼近了解析解［2］.

其中的系數(shù)是:當(dāng)n≠7時(shí)，

當(dāng)n=7時(shí)另外，B=0.利用MATLABn軟件也可以清晰的比較解析解中級(jí)數(shù)項(xiàng)數(shù)的多少對(duì)解的精確度的影響，當(dāng)級(jí)數(shù)達(dá)到50項(xiàng)以上才可以有較高的精確度，另外，根據(jù)模擬結(jié)果可以觀察到解析解中級(jí)數(shù)的每一項(xiàng)是一個(gè)駐波，不同頻率的駐波疊加成行波，生動(dòng)形象.

2 結(jié)束語(yǔ)

以上利用MATLAB軟件及有限差分方法求解兩端固定的有界弦振動(dòng)問(wèn)題，首先利用區(qū)域轉(zhuǎn)化的思想對(duì)方程進(jìn)行離散，然后編制MATLAB程序進(jìn)行求解，并將數(shù)值解可視化，使對(duì)問(wèn)題(1)的解的物理意義有了更深刻的認(rèn)識(shí)，另外，通過(guò)與解析解比較，也充分說(shuō)明了所選擇數(shù)值方法的可靠性和計(jì)算機(jī)輔助分析的重要性.

［1］ 李韻，郭怡文，呂郁文.基于 Matlab環(huán)境的弦振動(dòng)方程的圖像與音效模擬.科協(xié)論壇［J］，2009(7):74-75.

［2］ 谷超豪，李大潛，等.數(shù)學(xué)物理方程:第三版［M］.北京:高等教育出版社，2012.

［3］ 孫麗男，張馥菊.基于MATLAB的一類輸運(yùn)問(wèn)題的數(shù)值分析.哈爾濱師范大學(xué)自然科學(xué)學(xué)報(bào)，2012，28(1):14-17.

［4］ 高峰，陳君若，等.有界弦振動(dòng)方程的有限元方法求解［J］.機(jī)械設(shè)計(jì)，2008，25(11):15-18.