齊次弦振動(dòng)方程的MATLAB解法

1、真誠(chéng)為您提供優(yōu)質(zhì)參考資料，若有不當(dāng)之處，請(qǐng)指正。齊次弦振動(dòng)方程的MATLAB解法【摘要】弦振動(dòng)問題是一個(gè)典型的波動(dòng)方程的建立與求解問題。本文通過利用MATLAB特有的方程求解與畫圖功能，有效地構(gòu)造和求解了齊次弦振動(dòng)方程。并通過圖像，可以直觀感受方程的解，從而加深對(duì)這一問題物理意義的理解?！娟P(guān)鍵詞】振動(dòng)方程 MATLAB求解 數(shù)學(xué)物理方法【正文】在細(xì)弦上任意取微元分析其受力情況，通過Newton定律建立細(xì)弦振動(dòng)的運(yùn)動(dòng)方程，可以求得弦振動(dòng)的泛定方程為。要得出振動(dòng)方程的解，除了泛定方程外，我們還需要知道具體問題的初始條件與邊界條件。在弦振動(dòng)問題里，初始條件可以從初始位移和初始速度考慮，即：邊界條件是

2、描述物理問題在邊界上受約束的狀態(tài)，在弦振動(dòng)方程里可以歸結(jié)為三類邊界問題：（1） 第一類邊界問題： 稱為固定端。（2） 第二類邊界問題：特別的，若，稱為自由端。（3） 第三類邊界問題：第一類和第二類邊界問題的線性組合。一、 兩端固定的弦振動(dòng)問題兩端固定的弦振動(dòng)方程的定解問題可表示如下：1、初始位移不為0，初始速度為0不妨設(shè)：，（1）特征函數(shù)求解解由dAlembert公式：從而我們可以得到方程的級(jí)數(shù)解：而我們知道，弦振動(dòng)的泛定方程屬于本征問題：它在兩個(gè)邊界上都有第一類其次邊界條件，它的本征值與本征函數(shù)為：將系數(shù)帶入方程，級(jí)數(shù)中每一項(xiàng)都是一個(gè)駐波，定義子程序wfun.m計(jì)算不同n的求和各項(xiàng)，再用主程

3、序jxj將它們加起來(lái)，得到動(dòng)畫圖形。（MATLAB代碼見附錄1（1）（2）差分方程求解利用差分方程同樣可以求出問題的解。令，將微分方程改寫成差分方程，即有其中，于是，初始條件可以表示為：作圖時(shí)，先畫出的圖形，然后再用或代替其中的，改變的值，就畫出了不同時(shí)刻，的圖形。（MATLAB代碼見附錄1（2）解得的動(dòng)態(tài)圖形如下：2、初始位移為0，初始速度不為0設(shè)初始速度為：（1） 特征函數(shù)求解通過求本征函數(shù)與本征值的方法我們可以得到方程的解析解：其中系數(shù)，類似的，用函數(shù)計(jì)算級(jí)數(shù)中的各項(xiàng)，再在主函數(shù)中調(diào)用便可得解。（MATLAB代碼見附錄2（1）（2） 差分方程求解類似于問題1，我們還可以采用差分方程求解，

4、不過需要注意的是，題目中的初始條件應(yīng)表示為：。（MATLAB代碼見附錄2（2）解得的動(dòng)畫圖形如下：【總結(jié)】通過運(yùn)用MATLAB構(gòu)造和求解齊次弦振動(dòng)方程，繪制了相關(guān)圖像，直觀感受了方程解，加深了對(duì)其物理意義的理解。借助于計(jì)算機(jī)來(lái)做計(jì)算和研究的過程涉及到建立模型，選擇方法，語(yǔ)言編程和結(jié)果分析。通過此次問題的探究，培養(yǎng)和訓(xùn)練了自學(xué)能力和操作能力，獲益匪淺?！緟⒖嘉墨I(xiàn)】1、 李明奇 田太心 數(shù)學(xué)物理方程 電子科技大學(xué)出版社 20102、 彭芳麟 數(shù)學(xué)物理方程的MATLAB解法與可視化 清華大學(xué)出版社 20043、 彭芳麟 計(jì)算物理基礎(chǔ) 高等教育出版社 20104、 謝進(jìn) 李大美 MATLAB與計(jì)算方法

5、實(shí)驗(yàn) 武漢大學(xué)出版社 2009【附錄】附錄1（1）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)

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

7、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）

8、function psiN=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);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); drawnow; pau