第二章 常微分方程的初值問題

1、第二章 常微分方程的初值問題本章要研究的物理問題：經(jīng)典動力學(xué)中的有序和混沌本章內(nèi)容4簡單方法123多步法和隱式法龍格庫塔法穩(wěn)定性問題4動力學(xué)中的有序和混沌52.0 引子常微分方程是物理學(xué)中經(jīng)常碰到，以一維運(yùn)動粒子為例一般形式為M個藕合的一階方程實(shí)際問題往往涉及不止一類問題，例如對偏微分方程分離變量時(shí)，同時(shí)得到分離開來的初值問題與本征問題。本章只討論初值問題 常微分方程定解問題的分類本征值問題本征值問題邊界值問題邊界值問題初值問題給定待求函數(shù)在某個初始點(diǎn)上的值在自變量的兩個端點(diǎn)上對待求函數(shù)施加約束含有待定參數(shù)，只有在參數(shù)取特定值時(shí)，方程才有非零解2.1 簡單方法1). 離散化：將區(qū)間 0, 1

2、分為 N 個等間隔的子區(qū)間，每個區(qū)間寬度為 h=1/N待求問題求 x=1 處 y 的值策略2). 尋找一個遞推關(guān)系，把 yn 同 yn-1, yn-2, 聯(lián)系起來。局部誤差為 O(h2)在 xn 點(diǎn)，將微分方程左邊的微分用向前差分公式替代歐拉法得到遞推關(guān)系全局誤差為 NO(h2)O(h)精度太低！例子：牛頓動力學(xué)方程一維運(yùn)動的粒子方程為用歐拉法寫出具體的遞推關(guān)系歐拉法的Matlab實(shí)現(xiàn)標(biāo)量形式矢量形式for n = 1:N x(n+1) = x(n) + h * p(n) /m p(n+1) = p(n) + h * g(x(n), t(n)endf=(x, t)x(2)/m, g(x(1),

3、 t)for n = 1:N z(n+1,:) = z(n,:) + h*f(z(n,:), t(n)end這個公式在已知 f 的解析形式時(shí)很好用。階數(shù)增大時(shí)，變得復(fù)雜。泰勒級數(shù)法將 yn+1 在 yn 附近做泰勒展開而又已知從而有練習(xí)：推導(dǎo)三階泰勒級數(shù)法2.2 多步法和隱式法在子區(qū)間 xn, xn+1 將微分方程寫成積分形式關(guān)鍵在于對積分號下的 f(x, y) 取合理的近似歐拉法實(shí)際上將 f(x, y) 近似為 f(xn, yn) 利用 f 在 xn-1 和 xn 處的值通過線型插值得到 f 在積分區(qū)域的值多步法為了獲得更高的精度，可以將 yn+1 不僅僅是同 yn，而且還同更早的點(diǎn)，如 y

4、n-1, yn-2 等點(diǎn)相聯(lián)系例如:代入積分公式得Adams-Bashforth 二步積分法類似的，Adams-Bashforth 四步積分法為多步法有個小麻煩：最初幾個格點(diǎn)的啟動值需要通過其它的積分法如歐拉法、泰勒級數(shù)法來獲得隱式法在區(qū)間 xn , xn+1 上過 fn， fn+1 兩點(diǎn)對 f 插值，得帶入積分公式得隱式法意味著在每一個積分步都必須解一個方程，會非常耗時(shí)。一個特別簡單的情況是，如果 f 對于 y 是線性的，即 f(x,y)=g(x) y，則上面的方程可以為得到其顯式解利用 fn-1, fn, fn+1 三個值在區(qū)間 xn , xn+1 對 f 進(jìn)行二次多項(xiàng)式插值，得到隱式遞推

5、關(guān)系A(chǔ)dams-Moulton方法多步隱式法帶入積分公式得先通過顯式法“預(yù)報(bào)”一個 yn+1的值，然后利用隱式法來校正得到一個更精確的值。這樣的算法有個優(yōu)點(diǎn)，它可以持續(xù)監(jiān)控積分的精度。用三次多項(xiàng)式插值，得到相應(yīng)的三步法為隱式法很少直接使用，它通常用在預(yù)報(bào)校正算法中一個常常使用的具有局域誤差 O(h5)的預(yù)報(bào)校正算法：顯式的 Adams-Bashforth 四步法Adams-Moulton三步法練習(xí)：簡諧振子方程2.3 Runge-Kutta 法思想：從歐拉法說起在子區(qū)間 xn, xn+1 將微分方程寫成積分形式中值定理平均斜率歐拉法就是用 xn 點(diǎn)的斜率近似 xn, xn+1 區(qū)間的平均斜率

6、Kave改進(jìn)的歐拉法K1、K2 分別是 xn xn+1 點(diǎn)處的斜率 但是由于 yn+1 待定，因此需要做“預(yù)報(bào)”用 xn xn+1 兩點(diǎn)斜率的平均值來近似 xn, xn+1 平均斜率 思想：為了提高精度，多取幾點(diǎn)的斜率值作為加權(quán)平均當(dāng)作平均斜率其中 i (i = 1,2,.,m) 和 ij (i = 2,3,.,m 且 j i) 是待定參數(shù)Runge-Kutta 法將上式展開到 O(hm)，得到 m 個方程，而有 m+m(m1)/2個待定參數(shù) i , ij. 所以還有靈活選擇的空間以 m=2 為例 , Runge-Kutta 法公式為其中將 K2 做泰勒展開到 O(h2) 項(xiàng)，得利用將 yn+

7、1在 xn 點(diǎn)附近做泰勒展開至 O(h2) 項(xiàng)，得兩式相比，有可選解為或若取得二階Runge-Kutta法二階Runge-Kutta法與二階泰勒級數(shù)法比較可以看出，相比二階泰勒級數(shù)法而言，二階Runge-Kutta法適用性更廣，使用也更為方便三階Runge-Kutta法最為常用的是四階Runge-Kutta法通常認(rèn)為四階Runge-Kutta法在效率和精度間達(dá)到了最好的平衡練習(xí)：四階Runge-Kutta法處理一階常微分方程f = (x, y)(-x*y); y(1) = 1; for i=1:n K1 =f(x(i), y(i); K2 =f(x(i) + h/2, y(i) + h*K1/

8、2 ); K3 =f(x(i) + h/2, y(i) + h*K2/2 ); K4(i)=f(x(i) + h, y(i) + h*K3 ); y(i+1)=y(i) + 1/6*h*(K1 + 2*K2 + 2*K3 + K4 ); end練習(xí)：二階Runge-Kutta法處理二階常微分方程2.4 算法的穩(wěn)定性以歐拉法的一個簡單擴(kuò)展為例對微分方程進(jìn)行積分時(shí)，一個首要的考慮是所用算法的數(shù)值穩(wěn)定性，也就是說，舍入誤差或數(shù)值計(jì)算中的其它誤差能被放大的程度。為了啟動上面的遞推關(guān)系，還需要 y1 的值，這可以由泰勒展開獲得當(dāng)把這一方法應(yīng)用到下述問題這個問題的解析解是 y=e-x設(shè)解為指數(shù)形式 y=A

9、rn，代入遞推關(guān)系得正根略小于1，它對應(yīng)于我們要找的按指數(shù)減小的解，但是負(fù)根略小于-1，因而它對應(yīng)于一個虛假的解上述方程的解為其大小隨n增大并且在格點(diǎn)上逐點(diǎn)發(fā)生振蕩。線性差分方程的通解正是這兩個指數(shù)解的一個線性組合。雖然可以精心安排初值 y0 和 y1， 使得當(dāng) x 值小時(shí)只呈現(xiàn)指數(shù)衰減的解，但是在遞推過程中的舍入誤差將會摻入一個小的“壞”解，它最后將增長到在解中占壓倒地位。一個好的經(jīng)驗(yàn)規(guī)則是，每當(dāng)積分一個隨著迭代過程急劇衰減的解時(shí)，就應(yīng)當(dāng)小心不穩(wěn)定性和舍入誤差。例 1. 強(qiáng)迫鐘擺一根長度為 l 鐘擺被限制在一個垂直的平面內(nèi)，在強(qiáng)迫外力 fd 和阻力 fr 的作用下振蕩運(yùn)動。鐘擺的運(yùn)動可以通過

10、 Newton 方程來描述 其中 fg = -mg sin 是重力在運(yùn)動方向的分力，a=l d2 /dt2 是沿切線方向的加速度， 是桿和垂線的夾角。2.5 動力學(xué)中的有序和混沌 那么 Newton 運(yùn)動方程就可以寫成如下形式假設(shè)強(qiáng)迫外 力為 , 阻力為并且取 (l/g)1/2 為時(shí)間單位。 其中 該運(yùn)動方程就可以化為一階方程組設(shè) q=0.5,b=0.9, 0 =2/3, t=100。 在這 種情況下,鐘擺運(yùn)動是一個有序的周期運(yùn)動 桿和垂線的夾角、角速度隨時(shí)間演化過程角速度與夾角的軌跡q=0.5, b=1.15, 0 =2/3, t=1000，在這 種情況下鐘擺運(yùn)動是一個無序運(yùn)動,呈現(xiàn)出分形的

11、特征桿和垂線的夾角、角速度隨時(shí)間演化過程角速度與夾角的軌跡例 2. 自激振動范德波爾方程范德波爾描述范德波爾描述非線性有阻尼的自激動振動系統(tǒng)其中 是一個正的小量自激系統(tǒng)能將非振動的能源通過系統(tǒng)本身的反饋調(diào)節(jié)吸收進(jìn)來，以補(bǔ)充被損耗的能量。自激系統(tǒng)的例子包括心臟等。VDP方程不能準(zhǔn)確的表示心臟振蕩圖形的細(xì)節(jié)，為了更準(zhǔn)確的表示心臟的波動，F(xiàn)itzHugh對VDP方程進(jìn)行了修改，提出了所謂的BVP方程，形式為例 3. 化學(xué)振蕩BZ反應(yīng)CA, CB， CC， CD 分別為四種化合物濃度。例 4. 二維粒子的運(yùn)動考慮一個單位質(zhì)量的粒子,它在一個位勢 V 中作二維運(yùn)動。設(shè)在 t 時(shí)刻,粒子的平面坐標(biāo)為(x,

12、 y), 它的共軛動量為(px , py), 則 Hamilton 量的形式為 粒子的軌跡就由坐標(biāo)和動量隨時(shí)間的演化由如下 Hamilton 方程規(guī)定該約束條件把粒子運(yùn)動的軌道限制在四維相空間的三維流形上. 它是四個耦合的一階微分方程組.這些方程使能量 E 守恒,即滿足約束條件 在可分離變量的系統(tǒng)中,位勢 V 具有 (x,y) 的可分離形式,即 V(x,y)=Vx(x)+Vy(y)，此時(shí) Hamilton 方程可寫為可分離變量的系統(tǒng)和由此可見, x 方向和 y 方向的運(yùn)動是相互不耦合的，每一個 Hamilton 量單獨(dú)都是一個運(yùn)動常數(shù)，其中 作為一個例子,考慮兩個坐標(biāo)方向的運(yùn)動都是簡諧運(yùn)動的情況,即位勢對應(yīng)的 Hamilton 方程為和該問題是精確可解的。 我們可以通過直接積分,求得滿足初始條件 (x0, y0 , px0 , py0 ) 的解為(x, px)平面和 (y, py) 平面的軌跡圖(x, y) 平面的運(yùn)動軌跡在 Henon-Heiles 位勢系統(tǒng)中，位勢的形式為 Henon-Heiles