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

2、比較復(fù)雜因而無法求得解析解的波動方程, 有限差分方法的數(shù)值解法給出了一種較為普適而精確的計算方法.關(guān)鍵詞: 波動方程; 有限差分法; 顯式差分格式; 數(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 波動方程的差分格式 .52.2 初始、邊界條件的差分格式 .53 計算示例

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

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

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

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

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

10、得到(1.6)式, 它們都截斷 df xf xhf xdxh于項, 而把項以及更高冪次的項全部略去. (1.7)式相當(dāng)于把相應(yīng)的泰勒公式 df xhdx2h(1.8)式及(1.9)式聯(lián)立起來, 得到: (1.10) 333223!df xd f xf xhf xhhhdxdx當(dāng)略去項以及更高冪次的項以后, 便得到3h 2df xf xhf xhdxh即(1.7)式. 很明顯, 以上三種差商表達式中, 以(1.7)式所示的差商的截斷誤差最小, 因為(1.7)式略去的是三階以上的無窮小, 而其他的是二階以上的無窮小.對于二階導(dǎo)數(shù), 同樣可用一階差商的差商來近似表示, 即 (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截斷到項, 略去了項以及更高冪次的項. 222d f xhdx4h偏導(dǎo)數(shù)也可以仿照上述方法表示為差商, 它用各離散點上函數(shù)的差商來近似代替該點的偏導(dǎo)數(shù), 將需求解的邊值問題轉(zhuǎn)化為一組差分方程, 而后根據(jù)差分方程組(代數(shù)方程組), 解出位于各離散點上的待求函數(shù)值, 便可得到相應(yīng)的數(shù)值解.由上可見, 有限差分法是以差分原理為基礎(chǔ)的一種數(shù)值方法, 它實質(zhì)上是將物理上連續(xù)

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

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

14、, 這種格式是收斂和穩(wěn)定的.22OOh1P 2.2 初始、邊界條件的差分格式初始、邊界條件的差分格式初始條件用向前差商代替, 有,0u xt (2.5)101,0iiu xuut初始條件變?yōu)? (2.6) 01112,iiiiiuR xuR xRx為找到誤差為邊界條件的差分式, 取等分坐標(biāo)軸, 結(jié)點坐標(biāo)為, 相2Ohxhix應(yīng)的有, 則有泰勒展開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 計算示例計算示例3.1 邊界條件為零邊界條件為零, 初始位移、初始速度不為零初始位移、初始速度不為零計算一維波動方程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é)點坐標(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計算程序如下:%計算程序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;%計算矩陣 ufor k=2:K-1 for ii=2:I-1圖 3.1 此振動圖像是時間每間隔為 0.05秒時的位移隨坐標(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計算結(jié)果如圖 3.1 所示, 圖中曲線至上而下時間間隔為 0.05 秒. 將計算程序稍作修改, 物理過程圖像如下:圖 3.2 此圖像為圖 3.1 所表示的振動過程中初始時刻(用點劃線表示)及時間為0.25s 時刻(用實線表示)的波動圖像.圖 3.3 此圖像

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

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

23、點劃線形狀, 經(jīng)過后, 各個質(zhì)點波1s0.75s0.25s動圖像變?yōu)閷嵕€形狀.為觀察實際振動時的情況, 可將上面程序修改為如下動畫程序: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 此圖像為振動過程中的截圖.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;%計算矩陣 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 初始條件為零初始條件為零, 邊界條件不為零邊界條件不為零計算一維波動方程:20ttxxuc u0000sin,0,00 xx ltttuAtuuxlu令 (其中單位為,單位為 , 單位為, 單位為), 1,1,2 ,1cAlhms cm slm編寫計算程序如下.%計算程序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%計算矩陣 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 此振動圖像是時間每間隔為 0.05秒時的位移隨坐標(biāo)的變化圖像.起始時為一橫軸上的一條直線, 最后為圖像中最左邊的曲線.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計算結(jié)果如圖 3.7 所示, 起始時間曲線的位置與橫軸吻合, 最終的振動位移如圖中最左面的曲線. 其中為間隔 0.05 秒時的振動曲線. 將計算程序稍作修改, 物理過程圖像如下:圖 3.8 此圖像為圖 3.7 所表示的振動過程中初始時刻(用點劃線表示)及時間為0.25s 時刻(用實線表示)的波動圖像.圖 3.9 此圖像為圖 3.7 所表示的振動過程中 0.25s 時刻(用點劃線表示)及時間為0.5s 時刻(用實線表示)的波動圖像.圖 3.10 此圖像為圖 3.7

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

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

30、面程序修改為如下動畫程序: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 此圖為振動過程中的截圖.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%計算矩陣 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 自由項不為零自由項不為零, 初始、邊界條件為零初始、邊界條件為零計算一維波動方程2cossinttxxxuc uAtl00000;0,00 xx l

32、tttuuuxlu令 (其中單位為,單位為 , 單位為, 單位為), 1,1,2 ,1cAlhms cm slm編寫計算程序如下.%計算程序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 此振動圖像是時間每間隔為 0.05秒時的位移隨坐標(biāo)的變化圖像.起始時為一橫軸上的一條直線, 最后為圖像中左部最下面的曲線,右部為最上面的曲線.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;%計算矩陣 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計算結(jié)果如圖 3.13 所示, 起始時為一橫軸上的一條直線, 最后為圖像中左部最下面的曲線, 右部為最上面的曲線.將計算程序稍作修改, 物理過程圖像如下:此振動過程可以看做是兩端固定, 初始時刻的位移為零, 初始速度也為零的弦在驅(qū)動力的作用下的振動, 具體過程分析如下:cossinxAtl、如圖 3.14 所示, 點劃線表示振動過程的初始時刻的狀態(tài), 實線表示振動過程1的時刻的狀態(tài). 初始時刻各個質(zhì)點波動圖像為點劃線形狀, 經(jīng)過后, 各個0.25s0.25s質(zhì)點波動圖像變?yōu)閷嵕€形狀.、如圖 3.15 所示, 點劃線表示振動過程的時刻的狀態(tài)

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

36、75s 時刻(用實線表示)的波動圖像.圖 3.17 此圖像為圖 3.13 所表示的振動過程中 0.75s 時刻(用點劃線表示)及時間為1s 時刻(用實線表示)的波動圖像.20的時刻的狀態(tài). 時刻各個質(zhì)點波動圖像為點劃線形狀, 經(jīng)過后, 各個質(zhì)0.75s0.5s0.25s點波動圖像變?yōu)閷嵕€形狀.、如圖 3.17 所示, 點劃線表示振動過程的時刻的狀態(tài), 實線表示振動過程40.75s的時刻的狀態(tài). 時刻各個質(zhì)點波動圖像為點劃線形狀, 經(jīng)過后, 各個質(zhì)點1s0.75s0.25s波動圖像變?yōu)閷嵕€形狀.為觀察實際振動時的情況, 可將上面程序修改為如下動畫程序:%計算 3 動畫程序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 此圖為振動過程中的截圖.21%設(shè)置左端第一類邊界條件u(1,:)=0;%設(shè)置右端第一類邊界條件u(I,:)=0;%計算矩陣 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 自由項、定解條件都不為零自由項、定解條件都不為零考慮如下波動方程： 22210,0,sin,00,0ttxxtxuc uxutu

39、 l ttu xuxx lx令(其中單位為,單位為 , 單位為,1,1,1,0.05,1,2cPLhThms cm s單位為), 則編寫計算程序如下:lm22%計算程序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%計算矩陣 ufor k=2:K-1圖 3.19 此振動圖像是時間每間隔 0.05 秒位移隨坐標(biāo)的變化圖像.初始時刻為橫軸上一條直線,終止時為圖像左下方的曲線.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將計算程序稍作修改, 物理過程圖像如下:圖 3.20 此圖像為圖 3.19 所表示的振動過程中初始時刻(用點劃線表示)及時間為0.25s 時刻(用實線表示)的波動圖像.圖 3.21 此圖像為圖 3.19 所表示的振動過程中 0.25s 時刻(用點劃線表示)及時間為0.5s 時刻(用實線表示)的波動圖像.圖 3.22 此圖像為圖 3.19 所表示的振動過程中 0.5s 時刻(用點劃線表示)及時間為0.75s 時刻(用實線表示)的波動圖像.圖 3.23 此圖像為圖 3.19 所表示的振動過

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