第一章 基本運算(計算物理matlab)

1、第一章 基本數(shù)學運算本章要學會計算的物理問題：雙原子分子的振動能級本章內(nèi)容插值14數(shù)值積分3數(shù)值微分2擬合45方程求根分子振動的半經(jīng)典量子化461.1 插值x x0 x1 x2 xny y0 y1 y2 yn給定一組離散的數(shù)據(jù)尋找一個解析形式的函數(shù) (x)， 滿足 (xi) = yi, i=1,2, n 問題的提出最常用的函數(shù)形式是多項式，稱為多項式插值。Y(x0,y0)(x1,y1)(xn,yn)x0 x1xnXx x0 x1 x2y y0 y1 y222102)(xaxaax222221012121100202010yxaxaayxaxaayxaxaa最直觀的方法以二階為例設(shè)插值函數(shù)為帶入

2、數(shù)據(jù)，得離散數(shù)據(jù)為解這個方程，即可得到相應系數(shù) a0, a1, a2。N階呢nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaa22101121211000202010 nnnxaxaxaax2210)(設(shè)插值函數(shù)為n階多項式x x0 x1 xny y0 y1 yn離散數(shù)據(jù)為帶入數(shù)據(jù)，得當n很大時，求解這個方程計算量太大，需另尋它法拉格朗日插值從一階說起x x0 x1y y0 y1插值方程為一直線方程，可表示為對離散數(shù)據(jù)為了推廣到高階，將其寫成更對稱的形式x0 x1A0 (x)10A1 (x)01函數(shù)值函數(shù)值函數(shù)函數(shù)節(jié)點節(jié)點1010 xxxxA0101xxxxA滿足其中0122

3、021()()()()xxxxAxxxx221100)()()()(yxAyxAyxAxyx0 x1x2A0 (x)100A1 (x)010A2 (x)001函數(shù)值函數(shù)值函數(shù)函數(shù)節(jié)點節(jié)點更進一步二階插值y1y2x0 x1XYOy=f(x)x2y01200102()()()()xxxxAxxxx2011210()()()()xxxxAxxxxjnjjyxAxy)()(0njiiijijxxxxxA0)(這稱之為拉格朗日多項式插值。一般的N階多項式插值x x0 x1 xny y0 y1 yn已知離散數(shù)據(jù)插值多項式為其中滿足過高階的插值可能導致嚴重的振蕩行為，即Runge現(xiàn)象。是否階數(shù)越高，效果越好

4、？怎樣改進？)(xf)(10 xL-5 5例：連續(xù)函數(shù)215)(xxf可以看出，L10(x)的誤差在區(qū)間兩端非常大在區(qū)間-5,5上取等距插值節(jié)點討論用分段低次插值，最簡單的就是分段線性插值不光滑！解決的辦法分段插值區(qū)間a,b 有離散點：a= x0 x1 三點兩點 注意誤差隨步長 h 的減小先減小再增大注意舍入誤差微分公式涉及兩個很接近的 f 值相減。當步長過小時，計算機的舍入誤差會使導數(shù)計算不準確。高階導數(shù)二階導數(shù)更高階的階導數(shù)以此類推這也容易從二階導數(shù)的定義直接得出1.4 數(shù)值積分牛頓-萊布尼茲公式被積函數(shù) f(x) 有解析表達式f(x) 的原函數(shù) F(x) 為初等函數(shù))()()(aFbFd

5、xxfIab但是這要求為什么要數(shù)值積分？1) f(x) 沒有解析表達式2) f(x)有解析表達式，但原函數(shù)不是初等函數(shù) ，例如我們面臨的的問題x0.10.20.30.40.5f(x)44.5688.5它們的原函數(shù)都不是初等函數(shù)數(shù)值積分的出發(fā)點積分轉(zhuǎn)化為求和將區(qū)間 a,b 分割為 n 等份，每個小區(qū)間的寬度為 h=(b-a)/n在每個區(qū)間 xi, xi+1 進行線性插值線性近似梯形法則二階多項式近似辛普生法則在區(qū)間 xi-1, xi+1 上對 f(x) 進行二階插值。不同積分方法的結(jié)果比較注意，此時誤差隨h減小而減小，舍入誤差并不重要，這是因為積分公式中，所有 f 的值的符號都一樣。選取更多的點

6、，進行更高階的插值可以得到更高階的算法，如 Bode規(guī)則（四階插值） 高階的算法simpson 3/8 算法（三階插值）過高階的插值可能導致嚴重的振蕩行為，從而給出被積函數(shù)不準確的插值。所以為了得到更高精度，往往用低階插值，同時減小 h。反常積分分為兩類： 積分區(qū)間有限,在積分區(qū)間內(nèi)被積函數(shù)有奇點 積分區(qū)間為無限反常積分的處理通過積分變量的變換，將反常積分變換為普通積分策略（1） 積分區(qū)間內(nèi)含有奇點的積分積分區(qū)間內(nèi)含有奇點的積分在 x=1 處有一個奇點，假設(shè)函數(shù) g 在這點的值有限，則積分為有限值。做變換 t=(1-x)1/2，積分變?yōu)椋?）無限區(qū)間的積分）無限區(qū)間的積分g(x) 在 x 很大

7、時趨于常數(shù)，積分為有限值。做變換 t=x-1，積分變?yōu)?quad 用自適應辛普森算法。根據(jù)積分精度的需要，自動調(diào)節(jié)積分取點的數(shù)目。 調(diào)用格式為 quad(函數(shù)句柄, 積分下限，積分上限)quad(x)sin(x), 0.5, 0.6)Matlab自帶的積分指令1.5 數(shù)值求根高于四次方程的一般代數(shù)方程沒有一般形式的代數(shù)解更不用說更為復雜的方程1110( )(0)nnnnnf xa xaxa xaa( )( ) 10cos x cosh x 阿貝爾阿貝爾(18021829)求 f (x) = 0 的根原理：若 f Ca, b，且 f (a) f (b) 0，則 f (x)=0 在 (a, b)

8、上必有一根。yxbaf (x)x*1.5.1 二分法（搜索法）給定有根區(qū)間 a, b ( f(a) f(b) 0) 和 精度 1. 令 x = (a+b)/22. 如果 b a ， 停機，輸出 x3. 如果 f (a) f (x) 0 , 則令 b = x，否則令 a = x, 返回第1步二分法的算法實現(xiàn)abx1x2abx* 簡單易用 穩(wěn)妥，對 f (x) 要求不高，只要連續(xù)即可收斂 收斂速度慢 二分法的優(yōu)缺點 例1: 用二分法求方程 在區(qū)間 (1,2)內(nèi) 的實根, 要求誤差限為 。01523 xx21011.11.21.31.41.51.61.71.81.92-4-3-2-10123453(

9、 )251f xxx二分法例題 解：令 f (1)0 記 I0=1,2 , x0 =(1+2)/2=1.5 因為 f (x0) f (1)0 得 I1=1.5, 2 , x1 =(1.5+2)/2=1.75 f (x1) f (1.5)0 得 I2=1.5, 1.75 , x2 =(1.5+1.75)/2=1.625 . I6=1.681875, 1.6875, I7=1.671875, 1.679688 b7 - a7=0.781310-2 10-2 x*x7 =1.6758152)(3xxxf二分法例題 例1: 用二分法求方程 在區(qū)間 (1,2)內(nèi) 的實根, 要求誤差限為 。01523 x

10、x210 xyx*x0)()(1kkkkxfxfxx 0100()()f xxxfx 幾何意義1.5.2 牛頓法迭代形式為牛頓迭代法1: 初始化 x0 , ,置 i:=02: 如果| f(xi ) | ，則停止. 3: 計算 xi+1:=xi - f (xi) / f (xi)4: 如果 | f (xi+1) | ，則停止.5: i:=i+1, 轉(zhuǎn)至3.1()()iiiif xxxfx 牛頓法的算法構(gòu)造例1： 利用牛頓迭代法求解 f(x)=ex-1.5-tan-1x 的零點。初始點 x0=-7.0 解： f (x0)=-0.70210-1，f (x)=ex-(1+x2)-1 計算迭代格式: 計

11、算結(jié)果如下表：(取|f(x) |=10-10)k x f (x) 0 -7.0000 -0.07018881 -10.6771 -0.02256662 -13.2792 -0.004366023 -14.0537 -0.000239024 -14.1011 -7.99585e-0075 -14.1013 -9.00833e-0121()()kkkkf xxxfx 例1： 利用牛頓迭代法求解 f(x)=ex-1.5-tan-1x 的零點。初始點 x0=-7.0 Newtons Method 收斂性依賴于x0 的選取。x*x0 x0 x0算法說明x0 x1切線 割線 切線斜率割線斜率11)()()

12、( kkkkkxxxfxfxf)()()(111 kkkkkkkxfxfxxxfxx任意2個初值 x0 和 x1可以啟動這個遞推關(guān)系。1.5.3 弦割法三種求根算法的比較 二分法最穩(wěn)妥，但是效率最低。 牛頓法效率最高，但是要求函數(shù)解析，容易計算導數(shù)。 弦割法是前面兩種方法的折衷，既有較高的效率，又不必像牛頓法那樣必須計算函數(shù) f 的導數(shù)。 如果初始猜測比較接近待求解，則其收斂速度幾乎與牛頓算法一樣快。 三種求根算法的比較如果待求解附近，函數(shù)的行為不好，則自動的牛頓法和弦割法都可能無法收斂或收斂到錯誤的結(jié)果。保險的做法是先用二分法初步的定出解的位置，再用兩個自動方法中的一個定出解的精確位置。fz

13、ero(x)sin(x),3,4）Matlab自帶的求根指令X=fzero(函數(shù)句柄，猜測的初始值或搜尋的區(qū)間)fzero：求單變量函數(shù)的零點使用zeroin算法（結(jié)合了二分法、弦割法以及其它方法的一種綜合方法）。fzero(x)sin(x),3）1.6 分子振動的半經(jīng)典量子化原子的相互作用勢為 Lennard-Jones 或6-12形式勢能最低處為rmin=21/6a，深度為V0能量為En的相對運動的振動態(tài)可以用一維薛定諤方程的束縛態(tài)解n(r)來描述目標：對給定的勢，求得能量En標準辦法：數(shù)值求解常微分方程本征問題我們這里采用的方法：半經(jīng)典量子化通過考慮原子核在勢場中做經(jīng)典運動，然后應用量子

14、化規(guī)則，可以定出其振動能 En。這就是所謂的Bohr-Sommerfeld-Wilson量子化法則。半經(jīng)典量子化Bohr 根據(jù)對應原理的思想得出了一個角動量子化的條件，即電子運動的角動量 J 只能是 的整數(shù)倍Bohr的量子化法則Somemerfeld 等為處理多自由度體系的周期運動的能量量子化，給出了推廣的量子化條件其中其中 qk, pk 代表一對共軛的正則坐標與動量， 代表對周期運動積分一個周期。推廣的量子化法則在相空間軌跡為勢場 V 中，原子核間距的經(jīng)典的運動，可以在能量 -V0E0上發(fā)生。原子核間距在 rin 和 rout 之間周期性的振動，對應相空間中的一條封閉軌道，軌道方程為其中n為非負整數(shù)量子化條件為即其中為了對 Lennard-Jones勢確定量子化條件，定義幾個無量綱量量子化條件變?yōu)槭菢硕然蟮奈粍?。無量綱化分子分子H2HDO221.724.8150若知道分子的轉(zhuǎn)動慣量（可以從分子轉(zhuǎn)動的能量得知）和離解能（把分子分解為組成它的兩個原子所需的能量），就能從實驗觀測測定參量 a 和 V0，從而定出 。問題1假設(shè)原子間相互作用勢為諧振子勢，利用半經(jīng)典量子化條件，解析求系統(tǒng)能級。其中問題2證明：L-J勢下，小幅振動的頻率之預期值為問題3對氫分子更適用的勢函數(shù)為調(diào)