1、齊次弦振動方程的MATLAB解法【摘要】弦振動問題是一個典型的波動方程的建立與求解問題。本文通過利用MATLAB特有的方程求解與畫圖功能，有效地構造和求解了齊次弦振動方程。并通過圖像，可以直觀感受方程的解，從而加深對這一問題物理意義的理解?！娟P鍵詞】振動方程 MATLAB求解 數學物理方法【正文】在細弦上任意取微元分析其受力情況，通過Newton定律建立細弦振動的運動方程，可以求得弦振動的泛定方程為。要得出振動方程的解，除了泛定方程外，我們還需要知道具體問題的初始條件與邊界條件。在弦振動問題里，初始條件可以從初始位移和初始速度考慮，即：邊界條件是描述物理問題在邊界上受約束的狀態，在弦振動方程里

2、可以歸結為三類邊界問題：（1） 第一類邊界問題： 稱為固定端。（2） 第二類邊界問題：特別的，若，稱為自由端。（3） 第三類邊界問題：第一類和第二類邊界問題的線性組合。一、 兩端固定的弦振動問題兩端固定的弦振動方程的定解問題可表示如下：1、初始位移不為0，初始速度為0不妨設：，（1）特征函數求解解由dAlembert公式：從而我們可以得到方程的級數解：而我們知道，弦振動的泛定方程屬于本征問題：它在兩個邊界上都有第一類其次邊界條件，它的本征值與本征函數為：將系數帶入方程，級數中每一項都是一個駐波，定義子程序wfun.m計算不同n的求和各項，再用主程序jxj將它們加起來，得到動畫圖形。（MATLA

3、B代碼見附錄1（1）（2）差分方程求解利用差分方程同樣可以求出問題的解。令，將微分方程改寫成差分方程，即有其中，于是，初始條件可以表示為：作圖時，先畫出的圖形，然后再用或代替其中的，改變的值，就畫出了不同時刻，的圖形。（MATLAB代碼見附錄1（2）解得的動態圖形如下：2、初始位移為0，初始速度不為0設初始速度為：（1） 特征函數求解通過求本征函數與本征值的方法我們可以得到方程的解析解：其中系數，類似的，用函數計算級數中的各項，再在主函數中調用便可得解。（MATLAB代碼見附錄2（1）（2） 差分方程求解類似于問題1，我們還可以采用差分方程求解，不過需要注意的是，題目中的初始條件應表示為：。（

4、MATLAB代碼見附錄2（2）解得的動畫圖形如下：【總結】通過運用MATLAB構造和求解齊次弦振動方程，繪制了相關圖像，直觀感受了方程解，加深了對其物理意義的理解。借助于計算機來做計算和研究的過程涉及到建立模型，選擇方法，語言編程和結果分析。通過此次問題的探究，培養和訓練了自學能力和操作能力，獲益匪淺?！緟⒖嘉墨I】1、 李明奇 田太心 數學物理方程 電子科技大學出版社 20102、 彭芳麟 數學物理方程的MATLAB解法與可視化 清華大學出版社 20043、 彭芳麟 計算物理基礎 高等教育出版社 20104、 謝進 李大美 MATLAB與計算方法實驗 武漢大學出版社 2009【附錄】附錄1（1

5、）function jxjN=50t=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);endfunction 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-

6、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); 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(

7、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 psiN=50;t=0:0.0

8、05: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);endfunction 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）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);for k=2:N set(h,XData,x,YData,u(:,2); drawno