1、 學(xué) 院 本科畢業(yè)論文（設(shè)計(jì)）題題 目目 波動(dòng)方程的數(shù)值解法研究 院院 系系 物理與電子工程學(xué)院 專專 業(yè)業(yè) 物理學(xué) 姓姓 名名 學(xué)學(xué) 號(hào)號(hào) 學(xué)習(xí)年限學(xué)習(xí)年限 2007 年 9 月至 2011 年 7 月 指導(dǎo)教師指導(dǎo)教師 申請(qǐng)學(xué)位申請(qǐng)學(xué)位 理學(xué) 學(xué)士學(xué)位 二零一一年五月二十三日波動(dòng)方程的數(shù)值解法研究摘 要: 本文利用有限差分法及其顯示差分格式表示了波動(dòng)方程的差分格式以及初始、邊界條件的差分格式, 進(jìn)而討論了不同初始、邊界條件下的波動(dòng)方程的解的情況, 并利用 MATLAB 數(shù)學(xué)軟件編程繪制出了各種解的圖像. 為了方便對(duì)振動(dòng)過(guò)程的進(jìn)一步研究, 還給出了振動(dòng)過(guò)程的動(dòng)畫演示. 對(duì)于方程初始、邊界條件

2、比較復(fù)雜因而無(wú)法求得解析解的波動(dòng)方程, 有限差分方法的數(shù)值解法給出了一種較為普適而精確的計(jì)算方法.關(guān)鍵詞: 波動(dòng)方程; 有限差分法; 顯式差分格式; 數(shù)值解Numerical Solution of Wave EquationABSTRACT: In this paper, using finite difference method and show difference scheme showed that the difference scheme for wave equation and the initial and boundary conditions, then discus

3、sed the different initial and boundary conditions of the wave equation under the circumstances, and using mathematical software program gives MATLAB images of each solution. In order to facilitate the further research of vibration process, also gives a vibration process animation. For more complex b

4、oundary conditions and therefore should not be obtained by analytic solution wave equation, provides a relatively simple method.KEYWORDS：Wave equation; finite difference method; show difference scheme; numerical solution 目 錄引 言.11 有限差分法.21.1 有限差分法.22 顯式差分格式.52.1 波動(dòng)方程的差分格式 .52.2 初始、邊界條件的差分格式 .53 計(jì)算示例

5、.83.1 邊界條件為零, 初始位移、初始速度不為零.83.2 初始條件為零, 邊界條件不為零 .123.3 自由項(xiàng)不為零 .173.4 自由項(xiàng)、定解條件都不為零 .21總 結(jié).26致 謝.26參考文獻(xiàn).26 1引引 言言在實(shí)際工作中,經(jīng)常涉及求解常微分方程邊值問(wèn)題、偏微分方程邊值問(wèn)題. 通常,方程的精確解難以求解,而只能得到求解區(qū)域內(nèi)某些點(diǎn)處解的近似值, 即數(shù)值解. 在數(shù)學(xué)上, 二階線性偏微分方程分為雙曲型、拋物型、橢圓型三種, 其中雙曲型方程的典型代表是波動(dòng)方程. 當(dāng)滿足合適的邊界條件和初始條件時(shí), 波動(dòng)方程是可以求解的. 常用的求解方法有分離變量法1,2,3,5,6,7、行波法1,3,5

6、,6、格林函數(shù)法1,3,4、積分變換法1,2,3等等. 其中最為常用的是分離變量法. 而當(dāng)所討論的問(wèn)題的邊界條件較為復(fù)雜時(shí), 分離變量法將不再便于計(jì)算. 因此, 對(duì)于不易求解解析解的波動(dòng)方程的數(shù)值解的研究和發(fā)展具有重要意義.有限差分法是一種較為成熟的數(shù)值解法, 它的主要思想是將連續(xù)的問(wèn)題離散化. 是將每一處導(dǎo)數(shù)由有限差分近似公式代替, 從而把求解偏微分方程的問(wèn)題轉(zhuǎn)換成求解代數(shù)方程的問(wèn)題的方法. 本文在介紹了有限差分法的基礎(chǔ)知識(shí)后給出了波動(dòng)方程及其初始、邊界條件的差分格式,列舉了一些不同定解條件下的波動(dòng)方程的數(shù)值求解方法及圖像.21 有限差分法有限差分法1.1 有限差分法有限差分法有限差分法是在

7、采用數(shù)值計(jì)算方法求解偏微分方程時(shí), 將每一處導(dǎo)數(shù)由有限差分近似公式代替, 從而把求解偏微分方程的問(wèn)題轉(zhuǎn)換成求解代數(shù)方程的問(wèn)題的方法. 有限差分法首先用差分代替微分方程中的微分, 將連續(xù)變化的變量離散化, 從而得到差分方程組的數(shù)學(xué)形式; 繼而再求解差分方程組.設(shè)函數(shù)為, 其獨(dú)立變量有一很小的增量, 則該函數(shù)的增量為: f xxxh f x (1.1) f xf xhf x稱為函數(shù)的一階差分, 它與微分不同, 因是有限量的差, 故稱為有限差分. 而一 f x階差分除以增量的商, 即為一階差商: f xh (1.2) f xf xhf xxh一階差分仍是獨(dú)立變量的函數(shù), 類似地,可以得到二階差分,

8、即x 2f x (1.3) 2f xf xhf x顯然, 只要上述增量很小, 差分與微分之間的h fx差異將很小. 由于一階導(dǎo)數(shù) (1.4) 0limxdf xf xdxx是無(wú)限小的微分除以無(wú)限小微分 0limxdf xf x商. 因此, 如圖 2.1.1 所示, 當(dāng)用 2、3 點(diǎn)的量來(lái)0limxdxx表示函數(shù)的導(dǎo)數(shù)時(shí), 可近似地用差分表示為: f x (向前差分) (1.5) df xf xf xhf xdxxh即用有限差分除以增量的商來(lái)表示, 也稱為差商. 同理, 當(dāng)用 f xx df xdx f xx1、2 點(diǎn)的量來(lái)表示函數(shù)的一階導(dǎo)數(shù)時(shí), 可近似表達(dá)為: f x (向后差分) (1.6)

9、 df xf xf xf xhdxxh而用 1、3 點(diǎn)的量來(lái)表示時(shí), 則有: fx (中心差分) (1.7) 2df xf xf xhf xhdxxhf xh f xf xhOxhxxh123x f x圖 2.1.1 差分運(yùn)算原理圖3可見(jiàn), 函數(shù)在處的導(dǎo)數(shù)可用三種差商近似表示, 但哪一種精確度更高, 則 f xx可以通過(guò)泰勒公式的展開(kāi)式來(lái)分析. (1.8) 22212!df xd f xf xhf xhhdxdx和 (1.9) 22212!df xd f xf xhf xhhdxdx可見(jiàn), 將(1.8)式中的二階以及二階以上項(xiàng)略去后, 便可得到, 這就是(1.5)式. 同理, 由(1.9)式可

10、得到(1.6)式, 它們都截?cái)?df xf xhf xdxh于項(xiàng), 而把項(xiàng)以及更高冪次的項(xiàng)全部略去. (1.7)式相當(dāng)于把相應(yīng)的泰勒公式 df xhdx2h(1.8)式及(1.9)式聯(lián)立起來(lái), 得到: (1.10) 333223!df xd f xf xhf xhhhdxdx當(dāng)略去項(xiàng)以及更高冪次的項(xiàng)以后, 便得到3h 2df xf xhf xhdxh即(1.7)式. 很明顯, 以上三種差商表達(dá)式中, 以(1.7)式所示的差商的截?cái)嗾`差最小, 因?yàn)?1.7)式略去的是三階以上的無(wú)窮小, 而其他的是二階以上的無(wú)窮小.對(duì)于二階導(dǎo)數(shù), 同樣可用一階差商的差商來(lái)近似表示, 即 (1.11) 222111

11、2xxd f xdf xdf xf xhf xdxxdxdxhhhf xhf xf xf xhhhhf xhf xf xhh上式相當(dāng)于把泰勒公式 (1.12) 242424224!d f xd f xf xhf xhf xhhdxdx截?cái)嗟巾?xiàng), 略去了項(xiàng)以及更高冪次的項(xiàng). 222d f xhdx4h偏導(dǎo)數(shù)也可以仿照上述方法表示為差商, 它用各離散點(diǎn)上函數(shù)的差商來(lái)近似代替該點(diǎn)的偏導(dǎo)數(shù), 將需求解的邊值問(wèn)題轉(zhuǎn)化為一組差分方程, 而后根據(jù)差分方程組(代數(shù)方程組), 解出位于各離散點(diǎn)上的待求函數(shù)值, 便可得到相應(yīng)的數(shù)值解.由上可見(jiàn), 有限差分法是以差分原理為基礎(chǔ)的一種數(shù)值方法, 它實(shí)質(zhì)上是將物理上連續(xù)

12、域的問(wèn)題變換為離散系統(tǒng)的問(wèn)題來(lái)求解, 也就是用網(wǎng)格狀離散化模型上各離散4點(diǎn)的數(shù)值解來(lái)逼近連續(xù)場(chǎng)域的真實(shí)解的.有限差分法的應(yīng)用范圍很廣, 不但能求解均勻媒質(zhì)和不均勻線性媒質(zhì)中的位場(chǎng), 而且還能求解非線性媒質(zhì)中的場(chǎng); 它不僅能求解恒定場(chǎng)或似穩(wěn)場(chǎng), 還能求解時(shí)變場(chǎng). 在邊值問(wèn)題的數(shù)值方法中, 此法是相當(dāng)簡(jiǎn)便的, 在計(jì)算機(jī)存儲(chǔ)容量容許的情況下, 有可能采用較精細(xì)的網(wǎng)格, 使離散化模型能較精細(xì)地逼近真實(shí)問(wèn)題, 獲得具有足夠精度的數(shù)值解.52 顯式差分格式顯式差分格式2.1 波動(dòng)方程的差分格式波動(dòng)方程的差分格式考慮如下一維波動(dòng)方程的定解問(wèn)題 (2.1) 2220max22120011122200, 0,0

13、ttuucfxattxtxu x tuRRx txxtuaubcxnua ubcxan 取進(jìn)行離散化, 節(jié)點(diǎn)坐標(biāo)為,xht (2.2)max1,2,11,2,ikxhiNitkkk節(jié)點(diǎn)處的函數(shù)為. 在點(diǎn), 、用中心差商代替, 則式(1.1)中,kikiu x tu, i k22ut22ux的偏微分方程變?yōu)? (2. 3)2111122122kkkkkkkiiiiiiicfuuuuuuh引入, 上式變?yōu)?Pc h (2. 4)1211121kkkkkkiiiiiiuuPfuuuP由和時(shí)刻的可直接求時(shí)刻的, 不必解聯(lián)立方程組, 故這種差分格式是顯1k ku1k u式的. 誤差為, 可以證明, 當(dāng)時(shí)

14、, 這種格式是收斂和穩(wěn)定的.22OOh1P 2.2 初始、邊界條件的差分格式初始、邊界條件的差分格式初始條件用向前差商代替, 有,0u xt (2.5)101,0iiu xuut初始條件變?yōu)? (2.6) 01112,iiiiiuR xuR xRx為找到誤差為邊界條件的差分式, 取等分坐標(biāo)軸, 結(jié)點(diǎn)坐標(biāo)為, 相2Ohxhix應(yīng)的有, 則有泰勒展開(kāi)1, 2,iu i (2.7)2312!3!iiiiihhuuhuuu (2.8)2312!3!iiiiihhuuhuuu (2.9)2322222!3!iiiiihhuuhuuu (2.10)2322222!3!iiiiihhuuhuuu由(2.7)

15、可得一階向前差商公式 (2.11) 1iiiduuuOhdxh由(2.8)可得一階向后差商公式 (2.12) 1iiiduuuOhdxh式(2.9)式(2.7), 可得二階向前差商公式2 (2.13) 221222iiiid uuuuOhdxh式(2.10)式(2.8), 可得二階向后差商公式2 (2.14) 212222iiiid uuuuOhdxh將二階向前差商公式(2.13)代入(2.7), 則:232112212121222!3!2222432iiiiiiiiiiiiiiiiih uuuhuuhuuhuuhuuuuOhuuhuuOh 由此可得一階向前差商 (2.15)212342iii

16、iduuuuOhdxh同理將二階向前差商公式(2.14)代入(2.8), 可得一階向后差商 (2.16)212432iiiiduuuuOhdxh7將(2.15)代入邊界條件后得111uaubcn21111112111111112111111211342243223242423iiiiiiiiiiiiiiiuuuaubchhaubububuhchaubuhcbubuhcbubuuhab 令得到左邊界條件為:1i (2.17)112131112423hcbubuuhab將(2.16)代入邊界條件同理可得222ua ubcn212222222122221222243224322423iiiiiiii

17、iiiuuua ubchha ub ub ub uhchcb ub uuhab 令得到右邊界條件為:in (2.18)22122222423nnnhcb ub uuhab83 計(jì)算示例計(jì)算示例3.1 邊界條件為零邊界條件為零, 初始位移、初始速度不為零初始位移、初始速度不為零計(jì)算一維波動(dòng)方程2222200001,0sin,0;1,100,1,0ttuucxttxuxx txl lu x txxtuuttt微分方程的差分格式為1211121kkkkkkiiiiiiuuPfuuuP其中令, , 因此有 (其中單位為,單位0,1,fc21Pc h0.05h 0.05hm為 , 單位為, 單位為),

18、上式變?yōu)閟 cm slm1111kkkkiiiiuuuu其節(jié)點(diǎn)坐標(biāo)為max1,2,11,2,ikxhiNitkkk其中, 并選.1 0.05 121N max21k其初始條件 01112,iiiiiuR xuR xRx可寫為01sinsin,2,3,211sin2,3,211111iiiuxihiuhihhiii邊界條件1121 31112423hcbubuuhab22122222423nnnhcb ub uuhab由于, 則邊界條件可寫為12120ccbb10kknuu9計(jì)算程序如下:%計(jì)算程序clear%設(shè)置參數(shù)c=1;P=1;L=1;h=0.05; T=1;tao=P*h/c;N=fix

19、(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1)=sin(ii-1).*h.*pi);u(ii,2)=sin(ii-1).*h.*pi)+(ii-1).*h.*tao.*(1-(ii-1).*h);end%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件u(I,:)=0;%計(jì)算矩陣 ufor k=2:K-1 for ii=2:I-1圖 3.1 此振動(dòng)圖像是時(shí)間每間隔為 0.05秒時(shí)的位移隨坐標(biāo)的變化圖像.1

20、0 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1:K plot(x,u(:,k),-k,LineWidth,2) hold onend%axis(0,1,0,1)xlabel (fontsize14bfx/cm)ylabel (fontsize14bfu/cm)grid on計(jì)算結(jié)果如圖 3.1 所示, 圖中曲線至上而下時(shí)間間隔為 0.05 秒. 將計(jì)算程序稍作修改, 物理過(guò)程圖像如下:圖 3.2 此圖像為圖 3.1 所表示的振動(dòng)過(guò)程中初始時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.25s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.3 此圖像

21、為圖 3.1 所表示的振動(dòng)過(guò)程中 0.25s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.5s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.4 此圖像為圖 3.1 所表示的振動(dòng)過(guò)程中 0.5s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.75s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.5 此圖像為圖 3.1 所表示的振動(dòng)過(guò)程中 0.75s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為 1s時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.11此過(guò)程可以看做是兩端固定初始時(shí)刻位移為, 初始速度為的弦的振sinx1xx動(dòng),具體過(guò)程分析如下:、如圖 3.2 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的初始時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程的1時(shí)刻的狀態(tài).初始時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線

22、形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)點(diǎn)0.25s0.25s波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.3 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程20.25s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.5s0.25s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.4 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程的30.5s時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)點(diǎn)0.75s0.5s0.25s波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.5 所示,點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程的40.75s時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為

23、點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)點(diǎn)波1s0.75s0.25s動(dòng)圖像變?yōu)閷?shí)線形狀.為觀察實(shí)際振動(dòng)時(shí)的情況, 可將上面程序修改為如下動(dòng)畫程序:clear%設(shè)置參數(shù)c=1;P=1;L=1;h=0.05; T=1;tao=P*h/c;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1)=sin(ii-1).*h.*pi);圖 3.6 此圖像為振動(dòng)過(guò)程中的截圖.12u(ii,2)=sin(ii-1).*h.*pi)+(ii-1)

24、.*h.*tao.*(1-(ii-1).*h);end%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件u(I,:)=0;%計(jì)算矩陣 ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1: K plot(x,u(:,k),-k,LineWidth,2) hold on comet(x,u(:,k) pause(0.5)end%axis(0,1,0,1)xlabel (fontsize14bfx/cm)ylabel (fontsize14bfu/cm)grid on3.

25、2 初始條件為零初始條件為零, 邊界條件不為零邊界條件不為零計(jì)算一維波動(dòng)方程:20ttxxuc u0000sin,0,00 xx ltttuAtuuxlu令 (其中單位為,單位為 , 單位為, 單位為), 1,1,2 ,1cAlhms cm slm編寫計(jì)算程序如下.%計(jì)算程序clear13%設(shè)置參數(shù)c=1;P=1;L=1;h=0.05; T=1;tao=P.*h/c;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1

26、)=0;u(ii,2)=0;end%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件for ii=1:Ku(I,ii)=1.*sin(2.*pi.*(ii-1).*tao);end%計(jì)算矩陣 ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); end圖 3.7 此振動(dòng)圖像是時(shí)間每間隔為 0.05秒時(shí)的位移隨坐標(biāo)的變化圖像.起始時(shí)為一橫軸上的一條直線, 最后為圖像中最左邊的曲線.14endu;for k=1:K plot(x,u(:,k),-k,LineWidth,2) hold onend%axi

27、s(0,1,0,1)xlabel (fontsize14bfx/cm)ylabel (fontsize14bfu/cm)grid on計(jì)算結(jié)果如圖 3.7 所示, 起始時(shí)間曲線的位置與橫軸吻合, 最終的振動(dòng)位移如圖中最左面的曲線. 其中為間隔 0.05 秒時(shí)的振動(dòng)曲線. 將計(jì)算程序稍作修改, 物理過(guò)程圖像如下:圖 3.8 此圖像為圖 3.7 所表示的振動(dòng)過(guò)程中初始時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.25s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.9 此圖像為圖 3.7 所表示的振動(dòng)過(guò)程中 0.25s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.5s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.10 此圖像為圖 3.7

28、 所表示的振動(dòng)過(guò)程中 0.5s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.75s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.11 此圖像為圖 3.7 所表示的振動(dòng)過(guò)程中 0.75s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為1s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.15此振動(dòng)過(guò)程可以看做是始時(shí)刻的位移為零, 初始速度為零, 左端固定, 右端位移為的弦的振動(dòng), 具體過(guò)程分析如下:sinAt、如圖 3.8 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的初始時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程的1時(shí)刻的狀態(tài). 初始時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.25s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.9 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻

29、的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程20.25s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.5s0.25s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.10 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程30.5s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.75s0.5s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.11 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程40.75s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)1s0.75s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.為觀察實(shí)際振動(dòng)時(shí)的情況, 可將上

30、面程序修改為如下動(dòng)畫程序:clear%設(shè)置參數(shù)v=1;P=1;L=1;h=0.02; T=1;tao=P*h/v;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1; 圖 3.12 此圖為振動(dòng)過(guò)程中的截圖.16for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件for ii=1:Ku(I,ii)=1.*sin(2*pi.*(ii-1).*tao);end%計(jì)算矩陣 ufo

31、r k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1); endendu;for k=1:K axis(0,1,-1,1) plot(x,u(:,k),-k,LineWidth,2) hold on comet(x,u(:,k) pause(0.5)17endxlabel (fontsize14bfx/cm)ylabel (fontsize14bfu/cm)grid on3.3 自由項(xiàng)不為零自由項(xiàng)不為零, 初始、邊界條件為零初始、邊界條件為零計(jì)算一維波動(dòng)方程2cossinttxxxuc uAtl00000;0,00 xx l

32、tttuuuxlu令 (其中單位為,單位為 , 單位為, 單位為), 1,1,2 ,1cAlhms cm slm編寫計(jì)算程序如下.%計(jì)算程序clear%設(shè)置參數(shù)A=1;omega=2*pi;c=1;P=1;L=1;h=0.05; T=1;tao=P*h/c;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值圖 3.13 此振動(dòng)圖像是時(shí)間每間隔為 0.05秒時(shí)的位移隨坐標(biāo)的變化圖像.起始時(shí)為一橫軸上的一條直線, 最后為圖像中左部最下面的曲線,右部為最上面的曲線.18I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros

33、;%設(shè)置初始條件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件u(I,:)=0;%計(jì)算矩陣 ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1) .+tao2.*A.*cos(ii-1)*h.*pi/L).*sin(omega.*(k-1).*tao); endendu;for k=1: K plot(x,u(:,k),-k,LineWidth,2) hold onendxlabel (fontsize14bfx/cm)ylabel

34、 (fontsize14bfu/cm)grid on19計(jì)算結(jié)果如圖 3.13 所示, 起始時(shí)為一橫軸上的一條直線, 最后為圖像中左部最下面的曲線, 右部為最上面的曲線.將計(jì)算程序稍作修改, 物理過(guò)程圖像如下:此振動(dòng)過(guò)程可以看做是兩端固定, 初始時(shí)刻的位移為零, 初始速度也為零的弦在驅(qū)動(dòng)力的作用下的振動(dòng), 具體過(guò)程分析如下:cossinxAtl、如圖 3.14 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的初始時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程1的時(shí)刻的狀態(tài). 初始時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)0.25s0.25s質(zhì)點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.15 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài)

35、, 實(shí)線表示振動(dòng)過(guò)程20.25s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.5s0.25s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.16 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程30.5s圖 3.14 此圖像為圖 3.13 所表示的振動(dòng)過(guò)程中初始時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.25s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.15 此圖像為圖 3.13 所表示的振動(dòng)過(guò)程中 0.25s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.5s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.16 此圖像為圖 3.13 所表示的振動(dòng)過(guò)程中 0.5s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.

36、75s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.17 此圖像為圖 3.13 所表示的振動(dòng)過(guò)程中 0.75s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為1s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.20的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.75s0.5s0.25s點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.17 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程40.75s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)點(diǎn)1s0.75s0.25s波動(dòng)圖像變?yōu)閷?shí)線形狀.為觀察實(shí)際振動(dòng)時(shí)的情況, 可將上面程序修改為如下動(dòng)畫程序:%計(jì)算 3 動(dòng)畫程序clear%設(shè)置參數(shù)A=

37、1;omega=2*pi;c=1;P=1;L=1;h=0.02; T=1;tao=P*h/c;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1;for ii=1:I x(ii)=(ii-1)*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1)=0;u(ii,2)=0;end 圖 3.18 此圖為振動(dòng)過(guò)程中的截圖.21%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件u(I,:)=0;%計(jì)算矩陣 ufor k=2:K-1 for ii=2:I-1 u(ii,k+1)=u(ii+1,k)+u(ii-

38、1,k)-u(ii,k-1)+. tao2.*A.*cos(ii-1)*h.*pi/L).*sin(omega.*(k-1).*tao); endendu;for k=1:K axis(0,1,-0.08,0.08) plot(x,u(:,k),-k,LineWidth,2) hold on comet(x,u(:,k) pause(0.1)endxlabel (fontsize14bfx/cm)ylabel (fontsize14bfu/cm)grid on3.4 自由項(xiàng)、定解條件都不為零自由項(xiàng)、定解條件都不為零考慮如下波動(dòng)方程： 22210,0,sin,00,0ttxxtxuc uxutu

39、 l ttu xuxx lx令(其中單位為,單位為 , 單位為,1,1,1,0.05,1,2cPLhThms cm s單位為), 則編寫計(jì)算程序如下:lm22%計(jì)算程序clear%設(shè)置參數(shù)c=1;P=1;L=1;h=0.05;T=1;omega=2*pi;tao =P*h/c;N=fix(L/h);M=fix(T/tao);%設(shè)置 u 矩陣及 x 的值I=N+1;K=M+1;for ii=1:Ix(ii)=(ii-1).*h;endu(I,K)=zeros;%設(shè)置初始條件for ii=1:Iu(ii,1)=0;u(ii,2)=(ii-1).*h.*(1-(ii-1).*h);endfor ii

40、=1:K%設(shè)置左端第一類邊界條件u(1,ii)=0;%設(shè)置右端第一類邊界條件u(I,ii)=sin(omega.*(ii-1).*tao);end%計(jì)算矩陣 ufor k=2:K-1圖 3.19 此振動(dòng)圖像是時(shí)間每間隔 0.05 秒位移隨坐標(biāo)的變化圖像.初始時(shí)刻為橫軸上一條直線,終止時(shí)為圖像左下方的曲線.23 for ii=2:I-1u(ii,k+1)=u(ii+1,k)+u(ii-1,k)-u(ii,k-1)+ tao2.*(ii-1).*h/(1+(ii-1).*h)2)2;endendu;for k=1:Kplot(x,u(:,k),-k,linewidth,2)hold onendxl

41、abel(fontsize14bfx/cm) ylabel(fontsize14bfu/cm) grid on將計(jì)算程序稍作修改, 物理過(guò)程圖像如下:圖 3.20 此圖像為圖 3.19 所表示的振動(dòng)過(guò)程中初始時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.25s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.21 此圖像為圖 3.19 所表示的振動(dòng)過(guò)程中 0.25s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.5s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.22 此圖像為圖 3.19 所表示的振動(dòng)過(guò)程中 0.5s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為0.75s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.圖 3.23 此圖像為圖 3.19 所表示的振動(dòng)過(guò)

42、程中 0.75s 時(shí)刻(用點(diǎn)劃線表示)及時(shí)間為1s 時(shí)刻(用實(shí)線表示)的波動(dòng)圖像.24此振動(dòng)過(guò)程可以看做是左端固定, 右端位移為, 初始時(shí)刻的位移為零, 初始sint速度為的弦在驅(qū)動(dòng)力的作用下的振動(dòng), 具體過(guò)程分析如下:1xx221xx、如圖 3.20 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的初始時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程1的時(shí)刻的狀態(tài). 初始時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)0.25s0.25s質(zhì)點(diǎn)波動(dòng)圖像變?yōu)閷?shí)線形狀.、如圖 3.21 所示, 點(diǎn)劃線表示振動(dòng)過(guò)程的時(shí)刻的狀態(tài), 實(shí)線表示振動(dòng)過(guò)程20.25s的時(shí)刻的狀態(tài). 時(shí)刻各個(gè)質(zhì)點(diǎn)波動(dòng)圖像為點(diǎn)劃線形狀, 經(jīng)過(guò)后, 各個(gè)質(zhì)0.5s0.25