模擬物理 第二章

1、模擬物理 第二章第1頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四主要內(nèi)容簡單方法：Euler方法多步法和隱式法Runge-Kutta方法例子：二維運(yùn)動(dòng)中的有序和混沌第2頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四最常遇到的任務(wù)許多物理定律通過微分方程表述。比如求微分方程的數(shù)值解是模擬物理系統(tǒng)時(shí)最常遇到的一項(xiàng)任務(wù)。第3頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四方程的表達(dá)形式常微分方程最一般的形式是一組M個(gè)耦合的一階方程X是自變量，y是M個(gè)因變量，f有M個(gè)分量。高階微分方程通過引入若干輔助函數(shù)可以寫成這種一階形式。第4頁，共81頁，2022年，5月20日，6點(diǎn)

2、35分，星期四牛頓方程-哈密頓方程例如，一個(gè)質(zhì)量為m的粒子在力場F(z)作用下的一維運(yùn)動(dòng)由二階方程描述若定義動(dòng)量方程可以改為第5頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四只要詳細(xì)討論求解一階方程組的方法就夠了。只有一個(gè)因變量的情況，很容易推廣到有多個(gè)因變量的情況。因此只討論單個(gè)因變量的情況。本章重點(diǎn)討論初值問題。即給定y(x=0)=y0，求y(x)第6頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四2.1 簡單方法設(shè)我們想求解帶有初始條件y(x=0)=y0的微分方程更具體地說，我們通常感興趣的是某一特定x值（比如說x=1）上的y值。第7頁，共81頁，2022年，5月20日

3、，6點(diǎn)35分，星期四解法概述總的策略是把區(qū)間0,1分成N個(gè)等間隔的子區(qū)間，每個(gè)子區(qū)間的寬度h=1/N。然后找出一個(gè)遞推公式，把yn同yn-1,yn-2,聯(lián)系起來。其中yn是對y(xn=nh)的近似。這樣一個(gè)遞推關(guān)系將允許對這個(gè)微分方程進(jìn)行從x=0到x=1的逐步積分。第8頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四Euler方法一個(gè)最簡單的算法是Euler法?？紤]xn點(diǎn)上的情況，并且把微商換成前向差分近似。遞推關(guān)系第9頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四誤差局部誤差O(h2)整體誤差：經(jīng)過N步迭代，y(1)上的誤差為NO(h2)O(h)這一誤差線性地隨h減小。為使

4、最后結(jié)果的不精確度減半，需要h減半，因而步數(shù)增加一倍。第10頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四例子考慮微分方程和初條件它的解析解為下面我們用Euler方法從x=0積分到x=3，步長由外界輸入。第11頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四開始輸入步長：h計(jì)算總步數(shù) N=3/h初始條件 x=0,Y=1.0迭代一步: y=y+h*fnf(x,y)X=x+h循環(huán)計(jì)數(shù)器i=i+1輸出誤差iNN結(jié)束Yfnf(x, y)=-x*y第12頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四精度一般來說Euler法精度太低。這使我們不能通過采用較大的h來減少步數(shù)，從而

5、減少計(jì)算時(shí)間。在上面的例子中，當(dāng)我們試圖積分到更大的x值時(shí)，這一缺點(diǎn)就更明顯。當(dāng)x1/h是，y=0通常我們采用更高階精度的算法第13頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四簡單的高階方法：Taylor級(jí)數(shù)方法一類簡單的高階方法可以由Taylor級(jí)數(shù)展開式導(dǎo)出對其中的導(dǎo)數(shù)可以進(jìn)一步處理第14頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四局部誤差O(h3)整體誤差O(h2)，比Euler方法的精度高一階。缺點(diǎn)：只有當(dāng)f的解析形式已知，并且足夠簡單可以求微商時(shí)，才可以使用。第15頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四2.2 多步法和隱式法達(dá)到更高精度的另一

6、種方法是：使yn+1不僅同yn相聯(lián)系而且同更早的點(diǎn)比如yn-1, yn-2,相聯(lián)系，構(gòu)造一個(gè)包含多步的遞推關(guān)系。第16頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四多步法推導(dǎo)多步法公式：對微分方程做一步積分f(x,y)未知。取一個(gè)可以解析積分的近似。用xn和xn-1點(diǎn)上的f來做一個(gè)線性近似。帶入積分二步法第17頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四更高階的方法可以通過用更高次的多項(xiàng)式外插而得出。例如，若f用一個(gè)與fn,fn-1,fn-2,fn-3擬合的多項(xiàng)式來外插，就得到“四步法”第18頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四由于多步法的遞推關(guān)系式包

7、含前面的好幾步，單單關(guān)于y0的信息不能使它啟動(dòng)。因此必須通過別的方法，比如Euler方法、Taylor級(jí)數(shù)方法，或者下面討論的Runge-Kutta方法，先得出前幾個(gè)格點(diǎn)上的y值。第19頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四上述方法都是“顯式的”。它意味著：yn+1是用已經(jīng)知道的yn直接給出的（迭代）。隱式法：求解一個(gè)方程來決定yn+1。它提供了達(dá)到更高的精度的另一個(gè)手段。第20頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四我們設(shè)兩個(gè)格點(diǎn)的中點(diǎn)xn+1/2=(n+1/2)h?？紤]方程使用對稱差分近似得到遞推公式但是yn+1出現(xiàn)在兩邊。必須解方程。第21頁，共81頁，

8、2022年，5月20日，6點(diǎn)35分，星期四解方程可能很花時(shí)間?；仡櫍呵蟾?，搜索法，牛頓法，弦割法。如果f對于y是線性的，比方說f(x,y)=g(x)y，那么方程可以解出：第22頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四多步隱式法二步法：對方程 ，使用通過fn-1,fn,fn+1的二次多項(xiàng)式擬合f。解析積分得到隱式遞推公式第23頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四用三次多項(xiàng)式內(nèi)插可以推出對應(yīng)的三步公式第24頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四隱式法真的要解方程嗎？很少以通過解隱式方程的方式使用。用于“預(yù)估-校正”算法的基礎(chǔ)。先以顯式法得出

9、yn+1的預(yù)估值，再通過隱式法對它加以校正，給出更好的近似值。第25頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四2.3 Runge-Kutta方法對微分方程進(jìn)行積分的算法是有些自由的。實(shí)際上，的確存在許多算法。每種算法都有其特點(diǎn)和優(yōu)點(diǎn)。一種非常方便和廣泛使用的方式是Runge-Kutta算法。它有不同階的精度。第26頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四推導(dǎo)二階公式對于 ，f用它在積分區(qū)間中點(diǎn)附近的Taylor級(jí)數(shù)展開式逼近。做解析積分其中hf來自零次項(xiàng)，一次項(xiàng)積分為0，誤差來自二次項(xiàng)。第27頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四怎樣處理yn+1

10、/2呢？用Euler方法產(chǎn)生它。誤差為O(h2)二階RK算法：第28頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四特點(diǎn)：它體現(xiàn)了把y的近似值帶入隱式表達(dá)式右邊的想法。優(yōu)點(diǎn)：它同Taylor級(jí)數(shù)法或隱式法同樣精確，但是并不對f加特殊的約束。不要求容易求微商或者f關(guān)于y是線性的。僅使用y在前面的一個(gè)格點(diǎn)上的值。第29頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四高階算法用高階算法計(jì)算xn到xn+1的積分，用高階算法估計(jì)積分區(qū)間中格點(diǎn)上的y值可以推出高階算法。第30頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四三階算法第31頁，共81頁，2022年，5月20日，6點(diǎn)35

11、分，星期四四階算法第32頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四經(jīng)驗(yàn)發(fā)現(xiàn)， 四階算法在精度和計(jì)算量之間給出最佳的折中。第33頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四使用RK方法的例子用4階Runge-Kutta方法積分方程：求x=3時(shí)y的值。第34頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第35頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第36頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四2.5 二維運(yùn)動(dòng)中的有序和混沌在物理學(xué)中使用計(jì)算機(jī)帶來的一個(gè)基本好處是：能夠處理不能解析求解的系統(tǒng)。通常，數(shù)值結(jié)果同我們通過研究可解模型所發(fā)

12、展出的直觀圖象在定性上一致。而定量的數(shù)值結(jié)果又實(shí)際興趣。在少數(shù)情況下，計(jì)算機(jī)結(jié)果否定了我們的直觀。這時(shí)數(shù)值結(jié)果對正確理解這種現(xiàn)象有根本的重要性。第37頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四2.5 二維運(yùn)動(dòng)中的有序和混沌基本好處： 能夠處理不能解析求解的系統(tǒng)。通常在少數(shù)情況下數(shù)值結(jié)果直觀圖象定性實(shí)際興趣定量結(jié)果數(shù)值結(jié)果直觀圖象第38頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四否定經(jīng)典圖象令人驚奇的是，數(shù)值結(jié)果也否定了一些經(jīng)典系統(tǒng)動(dòng)力學(xué)的直觀圖象。牛頓以后，一切運(yùn)轉(zhuǎn)正常相對論指出高速運(yùn)動(dòng)下經(jīng)典的直觀圖象不可行量子力學(xué)指出微觀下經(jīng)典的直觀圖象不可行計(jì)算機(jī)的應(yīng)用指出經(jīng)典情

13、形下，一些經(jīng)典的直觀圖象不可行第39頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四做二維運(yùn)動(dòng)的粒子本例題中，我們將通過用數(shù)值積分計(jì)算一個(gè)做二維運(yùn)動(dòng)的粒子的軌道以此為例研究這種令人驚奇的行為。考慮一個(gè)單位質(zhì)量的粒子。它在一個(gè)位勢V中作二維運(yùn)動(dòng)。假定如果粒子的能量足夠低，它將永遠(yuǎn)被約束在V中。第40頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四若與兩個(gè)坐標(biāo)(x, y)共軛的動(dòng)量是(px, py)，則Hamilton量的形式為粒子的軌道由坐標(biāo)和動(dòng)量隨時(shí)間的演化所確定。注意質(zhì)量第41頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四Hamilton方程坐標(biāo)和動(dòng)量的隨時(shí)間的演化

14、由四個(gè)耦合的一階微分方程規(guī)定，即Hamilton方程第42頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四約束對于任何V，這些方程都能使能量E守恒，約束條件：因此，約束條件把軌道限制在嵌在四維相空間中的一個(gè)三維流形上。除此之外，關(guān)于系統(tǒng)的演化很難作出別的什么一般性判斷。第43頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四可積系統(tǒng)能夠?qū)壍雷鞒鲞M(jìn)一步判斷的一類重要的二維Hamilton量是可積的Hamilton量。對于這種情況，除了能量以外，還有第二個(gè)坐標(biāo)和動(dòng)量的函數(shù)是運(yùn)動(dòng)常數(shù)。于是軌道被約束在相空間的一個(gè)二維流形上。兩種熟悉的可積系統(tǒng)是可分離變量的系統(tǒng)和中心位勢。第44頁，共

15、81頁，2022年，5月20日，6點(diǎn)35分，星期四可積系統(tǒng)1：可分離變量可分離變量的情況下其中Vx, Vy是兩個(gè)獨(dú)立的函數(shù)。因此Hamilton量分離為兩部分。兩個(gè)方向的運(yùn)動(dòng)沒有耦合。每一個(gè)方向上的Hamilton單獨(dú)是一個(gè)運(yùn)動(dòng)常數(shù)。（等價(jià)于，Hx-Hx和E=Hx+Hy）H=Hx+HyHx=px2/2+Vx , Hy=py2/2+Vy第45頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四可積系統(tǒng)2：中心位勢在中心位勢的場合向心力導(dǎo)致角動(dòng)量是第二個(gè)運(yùn)動(dòng)常數(shù)Hamilton量可以寫為其中pr是r的共軛動(dòng)量。第46頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四可積性的困難可積系統(tǒng)的

16、動(dòng)力學(xué)是簡單的。但是要找出這種簡單性常常很不容易。沒有一個(gè)普遍的解析方法可以判斷，在任意一個(gè)位勢中是否存在第二個(gè)運(yùn)動(dòng)常數(shù)，以及如果有的話如何求出它。第47頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四數(shù)值計(jì)算的固有局限數(shù)值計(jì)算看起來也不會(huì)改善情況。因?yàn)閿?shù)值計(jì)算只給出給定初始條件下的軌跡。而這條軌跡即使在我們熟悉的情況下也可以是相當(dāng)復(fù)雜的。比如Lissajous圖形。（示例）第48頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第49頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第50頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四分析相空間通過對相空間的分析

17、，可以得出一個(gè)根據(jù)軌跡來檢測可積性的辦法。考慮一個(gè)可分離變量的位勢的情況。因?yàn)閮蓚€(gè)坐標(biāo)上的運(yùn)動(dòng)是獨(dú)立的，在(x, px)平面和(y, py)平面上畫出的軌跡看起來可能是第51頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四可積性的標(biāo)志一個(gè)粒子在一個(gè)可分離變量的二維位勢中的軌跡在(x, px)和(y, py)平面上的樣子。這些閉合回路的存在是系統(tǒng)可積性的標(biāo)志。第52頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四數(shù)值做法：截面根據(jù)軌跡可以得到一幅低維空間上的圖。假設(shè)每當(dāng)我們觀察到一個(gè)坐標(biāo)例如x通過零時(shí)，我們就在(y, py)平面上畫出粒子的位置。如果x運(yùn)動(dòng)和y運(yùn)動(dòng)的周期是不可通約

18、的（即它們的比值是一個(gè)無理數(shù)），那么隨著軌跡前進(jìn)，這些觀測值將逐漸描繪出完整的(y, py)回路。如果它們的周期是可通約的（即比值是有理數(shù)）那么將得到沿回路的一系列分立點(diǎn)。第53頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四利用截面方法，可以只從軌跡出發(fā)，來研究同任何給定的Hamilton量相聯(lián)系的相空間的拓?fù)湫再|(zhì)。第54頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四在一個(gè)可積的Hamilton量的情形下，相空間的一般的拓?fù)湫再|(zhì)，可以通過考慮中心位勢中的運(yùn)動(dòng)來說明。對于固定的能量值和角動(dòng)量值，徑向運(yùn)動(dòng)被限制在兩個(gè)轉(zhuǎn)折點(diǎn)rin和rout之間。這兩個(gè)r值是方程的兩個(gè)解： pr=

19、0 例如：衛(wèi)星軌道第55頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四這兩個(gè)半徑在(x, y)平面中定義了一個(gè)圓環(huán)面。對于給定的r值，能量守恒允許的徑向運(yùn)動(dòng)的動(dòng)量只能是這兩個(gè)動(dòng)量在(x, y, pr)空間中確定了包含全部軌道的二維流形第56頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第57頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四如果通過考慮x=0平面來畫出一幅(y, py)截面圖，我們將會(huì)得到兩個(gè)閉合回路。畫截面比畫流形簡單的多。x=0第58頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四可以證明，中心位勢的相空間的環(huán)面拓?fù)浣Y(jié)構(gòu)，是一切可積系統(tǒng)所共

20、有的。對于給定的運(yùn)動(dòng)常數(shù)，軌跡所在的流形叫做一個(gè)“不變環(huán)面”。第59頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四環(huán)面的一般的截面圖的樣子第60頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四對可積性的擾動(dòng)問題：在一個(gè)破壞掉系統(tǒng)的可積性的擾動(dòng)下，一個(gè)可積系統(tǒng)的環(huán)面會(huì)發(fā)生什么變化？對于小擾動(dòng)，橢圓型不動(dòng)點(diǎn)周圍的大部分環(huán)面發(fā)生輕微的畸變，但仍保持它們的拓?fù)浣Y(jié)構(gòu)（KAM定理）。但是相空間中鄰接區(qū)域則變成混沌的，給出的截面圖看起來像是隨機(jī)撒的一些點(diǎn)子。在這些混沌區(qū)域中又嵌入別的橢圓型不動(dòng)點(diǎn)和別的混沌區(qū)域，構(gòu)成層次結(jié)構(gòu)。第61頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四He

21、non-Heiles位勢對可積性的重大偏離必須用數(shù)值方法研究。進(jìn)行這一研究的一個(gè)方便的例子是這個(gè)位勢原來是Henon和Heiles在研究恒星穿過星系的軌道時(shí)引入的。恒星的運(yùn)動(dòng)被認(rèn)為限制在二維平面上。第62頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四第63頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四Hamiltion運(yùn)動(dòng)方程第64頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四有了運(yùn)動(dòng)方程，給定一組初始條件（x, y, px, py）就可以計(jì)算出一條軌道。這里使用4階Runge-Kutta方法。第65頁，共81頁，2022年，5月20日，6點(diǎn)35分，星期四選擇初條件考慮束縛軌道總能量 E1/6為了方便滿足這個(gè)條件，我們直接選定E。第66頁，共81頁