基本運算計算物理

基本運算計算物理第1頁，共73頁，2023年，2月20日，星期四本章內(nèi)容插值14數(shù)值積分3數(shù)值微分2擬合45方程求根分子振動的半經(jīng)典量子化46第2頁，共73頁，2023年，2月20日，星期四1.1插值xx0

x1

x2……

xnyy0

y1

y2

……

yn給定一組離散的數(shù)據(jù)尋找一個解析形式的函數(shù)φ(x)，滿足φ(xi)=yi,i=1,2,…,n問題的提出最常用的函數(shù)形式是多項式，稱為多項式插值。Y(x0,y0)(x1,y1)(xn,yn)x0x1xnX第3頁，共73頁，2023年，2月20日，星期四xx0

x1

x2yy0

y1

y2最直觀的方法——以二階為例設(shè)插值函數(shù)為帶入數(shù)據(jù)，得離散數(shù)據(jù)為解這個方程，即可得到相應(yīng)系數(shù)a0,a1,a2。第4頁，共73頁，2023年，2月20日，星期四N階呢設(shè)插值函數(shù)為n階多項式xx0

x1

……

xnyy0

y1

……

yn離散數(shù)據(jù)為帶入數(shù)據(jù)，得當(dāng)n很大時，求解這個方程計算量太大，需另尋它法第5頁，共73頁，2023年，2月20日，星期四拉格朗日插值——從一階說起xx0

x1yy0

y1插值方程為一直線方程，可表示為對離散數(shù)據(jù)第6頁，共73頁，2023年，2月20日，星期四為了推廣到高階，將其寫成更對稱的形式x0x1A0

(x)10A1

(x)01函數(shù)值函數(shù)節(jié)點滿足其中第7頁，共73頁，2023年，2月20日，星期四x0x1x2A0

(x)100A1

(x)010A2

(x)001函數(shù)值函數(shù)節(jié)點更進(jìn)一步——二階插值y1y2x0x1XYOy=f(x)x2y0第8頁，共73頁，2023年，2月20日，星期四這稱之為拉格朗日多項式插值。一般的——N階多項式插值xx0

x1

……

xnyy0

y1

……

yn已知離散數(shù)據(jù)插值多項式為其中滿足第9頁，共73頁，2023年，2月20日，星期四過高階的插值可能導(dǎo)致嚴(yán)重的振蕩行為，即Runge現(xiàn)象。是否階數(shù)越高，效果越好？怎樣改進(jìn)？-55例：連續(xù)函數(shù)可以看出，L10(x)的誤差在區(qū)間兩端非常大在區(qū)間[-5,5]上取等距插值節(jié)點討論第10頁，共73頁，2023年，2月20日，星期四用分段低次插值，最簡單的就是分段線性插值不光滑！解決的辦法——分段插值第11頁，共73頁，2023年，2月20日，星期四區(qū)間[a,b]有離散點：a=x0<x1<…<xn=b，及其對應(yīng)的函數(shù)值yi(i=0,1,…,n)，如果函數(shù)S(x)滿足條件：(2)在每個子區(qū)間[xi

,xi+1](i=0,1,…,n-1)上是三次多項式則稱S(x)是y=f(x)的三次樣條插值函數(shù)確定S(x)需要4n個條件,而我們只給出了4n-2個條件端點函數(shù)值：2個內(nèi)節(jié)點函數(shù)值及連續(xù)條件：2(n-1)

個內(nèi)節(jié)點一、二階導(dǎo)數(shù)連續(xù)條件：2(n-1)

個解決之道——在每個區(qū)間上，不是用線性函數(shù)而是用三階多項式進(jìn)行插值。(1)在[a,b]上具有二階連續(xù)導(dǎo)數(shù)，即(3)樣條插值第12頁，共73頁，2023年，2月20日，星期四需補充2個邊界條件，根據(jù)實際情況有不同的取法例如，可取周期邊界條件：第13頁，共73頁，2023年，2月20日，星期四樣條插值的例子第14頁，共73頁，2023年，2月20日，星期四Matlab指令interp1(X,Y,Xi，method)X,Y 數(shù)據(jù)點及其函數(shù)值Xi 待求的數(shù)據(jù)點method 插值方法，可取’linear’,’spline’等返回值

數(shù)據(jù)點Xi的函數(shù)值第15頁，共73頁，2023年，2月20日，星期四1.2擬合給定一組離散數(shù)據(jù)：要求一個函數(shù)p(x)，按照最小二乘原理，使得達(dá)到最小值，這一方法被稱為擬合，其中p(x)稱為擬合曲線。xx1

……

xnyy1

……

yn第16頁，共73頁，2023年，2月20日，星期四則問題變?yōu)榍骯k(k=0,1,2,…,m),使得一種常見的擬合曲線為多項式擬合，即選取擬合曲線p(x)為m次多項式極小。第17頁，共73頁，2023年，2月20日，星期四由多元函數(shù)取極值條件，得這就是正則方程。解這個方程即可得到擬合函數(shù)的形式。第18頁，共73頁，2023年，2月20日，星期四當(dāng)最高冪次為1時，即為最小二乘法的擬合直線，正則方程簡化為第19頁，共73頁，2023年，2月20日，星期四Matlab指令polyfit(X,Y,N)X,Y 數(shù)據(jù)點及其函數(shù)值N 擬合多項式的最高冪次返回值

為按降冪排列的多項次系數(shù)第20頁，共73頁，2023年，2月20日，星期四例子i123456789xi-1-0.75-0.5-0.2500.250.50.751yi-0.22090.32950.88261.43922.00032.56453.13343.76014.2836設(shè)所求的最小二乘二次擬合多項式為相應(yīng)的正則方程為第21頁，共73頁，2023年，2月20日，星期四其解為a0=2:0034;a1=2:2625;a2=0:0378，所求多項式為第22頁，共73頁，2023年，2月20日，星期四實例——密立根油滴實驗第23頁，共73頁，2023年，2月20日，星期四指數(shù)擬合如果數(shù)據(jù)點的分布近似于指數(shù)曲線，則可考慮用指數(shù)函數(shù)去擬合數(shù)據(jù)。同時對數(shù)據(jù)yi

也取對數(shù)得Inyi,利用數(shù)據(jù)組(xi,Inyi)求出最小二乘擬合直線。再取指數(shù)即得到原數(shù)據(jù)組的最小擬合指數(shù)對上式取對數(shù)第24頁，共73頁，2023年，2月20日，星期四例子ti0.20.30.40.50.60.70.8Ii3.162.381.751.341.000.740.56已知一發(fā)射源的發(fā)射強度具有指數(shù)形式I=I0e-αt,現(xiàn)有一組觀測數(shù)據(jù)如下第25頁，共73頁，2023年，2月20日，星期四WithfourparametersIcanfitanelephant,andwithfiveIcanmakehimwigglehistrunk.關(guān)于擬合，要提醒的第26頁，共73頁，2023年，2月20日，星期四1.3數(shù)值微分為什么要學(xué)習(xí)數(shù)值微分？我們碰到的函數(shù)往往沒有解析形式，例如可能是如下數(shù)表形式x0.10.20.30.40.5f(x)44.5688.5這就需要借助于數(shù)值微分更重要的，數(shù)值微分是其它很多數(shù)值方法的基礎(chǔ)第27頁，共73頁，2023年，2月20日，星期四求X=0處的導(dǎo)數(shù)x…-2h-h0h2h…f(x)…f-2f-1f0f1f2…已知第28頁，共73頁，2023年，2月20日，星期四方法：泰勒展開x=h處的函數(shù)值x=-h處的函數(shù)值第29頁，共73頁，2023年，2月20日，星期四兩點公式向前差分向后差分第30頁，共73頁，2023年，2月20日，星期四三點公式（中心差分公式）第31頁，共73頁，2023年，2月20日，星期四進(jìn)一步的改進(jìn)——5點公式雖然精度更大，但是計算量也更大！通常的計算三點公式已足夠。第32頁，共73頁，2023年，2月20日，星期四例子：不同微分方法的精度比較不同方法計算dsinx/dx|x=1=0.540302的誤差第33頁，共73頁，2023年，2月20日，星期四不同公式的精度：五點>三點>兩點注意誤差隨步長h的減小先減小再增大注意舍入誤差微分公式涉及兩個很接近的f值相減。當(dāng)步長過小時，計算機的舍入誤差會使導(dǎo)數(shù)計算不準(zhǔn)確。第34頁，共73頁，2023年，2月20日，星期四高階導(dǎo)數(shù)二階導(dǎo)數(shù)更高階的階導(dǎo)數(shù)以此類推這也容易從二階導(dǎo)數(shù)的定義直接得出第35頁，共73頁，2023年，2月20日，星期四1.4數(shù)值積分牛頓-萊布尼茲公式被積函數(shù)f(x)有解析表達(dá)式f(x)的原函數(shù)F(x)為初等函數(shù)但是這要求為什么要數(shù)值積分？第36頁，共73頁，2023年，2月20日，星期四1)f(x)沒有解析表達(dá)式2)f(x)有解析表達(dá)式，但原函數(shù)不是初等函數(shù)，例如我們面臨的的問題x0.10.20.30.40.5f(x)44.5688.5它們的原函數(shù)都不是初等函數(shù)第37頁，共73頁，2023年，2月20日，星期四數(shù)值積分的出發(fā)點積分轉(zhuǎn)化為求和將區(qū)間[a,b]分割為n等份，每個小區(qū)間的寬度為h=(b-a)/n第38頁，共73頁，2023年，2月20日，星期四在每個區(qū)間[xi,xi+1]進(jìn)行線性插值線性近似——梯形法則第39頁，共73頁，2023年，2月20日，星期四二階多項式近似——辛普生法則在區(qū)間[xi-1,xi+1]上對f(x)進(jìn)行二階插值。第40頁，共73頁，2023年，2月20日，星期四不同積分方法的結(jié)果比較注意，此時誤差隨h減小而減小，舍入誤差并不重要，這是因為積分公式中，所有f的值的符號都一樣。第41頁，共73頁，2023年，2月20日，星期四選取更多的點，進(jìn)行更高階的插值可以得到更高階的算法，如Bode規(guī)則（四階插值）高階的算法simpson3/8算法（三階插值）第42頁，共73頁，2023年，2月20日，星期四過高階的插值可能導(dǎo)致嚴(yán)重的振蕩行為，從而給出被積函數(shù)不準(zhǔn)確的插值。所以為了得到更高精度，往往用低階插值，同時減小h。第43頁，共73頁，2023年，2月20日，星期四反常積分分為兩類：積分區(qū)間有限,在積分區(qū)間內(nèi)被積函數(shù)有奇點積分區(qū)間為無限反常積分的處理通過積分變量的變換，將反常積分變換為普通積分策略第44頁，共73頁，2023年，2月20日，星期四（1）

積分區(qū)間內(nèi)含有奇點的積分在x=1

處有一個奇點，假設(shè)函數(shù)g

在這點的值有限，則積分為有限值。做變換t=(1-x)1/2，積分變?yōu)榈?5頁，共73頁，2023年，2月20日，星期四（2）無限區(qū)間的積分g(x)在x很大時趨于常數(shù)，積分為有限值。做變換t=x-1，積分變?yōu)榈?6頁，共73頁，2023年，2月20日，星期四quad

用自適應(yīng)辛普森算法。根據(jù)積分精度的需要，自動調(diào)節(jié)積分取點的數(shù)目。調(diào)用格式為quad(函數(shù)句柄,積分下限，積分上限)quad(@(x)sin(x),0.5,0.6)Matlab自帶的積分指令第47頁，共73頁，2023年，2月20日，星期四1.5數(shù)值求根高于四次方程的一般代數(shù)方程沒有一般形式的代數(shù)解更不用說更為復(fù)雜的方程阿貝爾(1802—1829)第48頁，共73頁，2023年，2月20日，星期四求f(x)=0的根原理：若f

C[a,b]，且f(a)·f(b)<0，則f(x)=0在(a,b)上必有一根。yxbaf(x)x*1.5.1二分法（搜索法）第49頁，共73頁，2023年，2月20日，星期四給定有根區(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*第50頁，共73頁，2023年，2月20日，星期四

簡單易用

穩(wěn)妥，對f(x)要求不高，只要連續(xù)即可收斂

收斂速度慢二分法的優(yōu)缺點第51頁，共73頁，2023年，2月20日，星期四

例1:用二分法求方程在區(qū)間(1,2)內(nèi)的實根,要求誤差限為。二分法——例題第52頁，共73頁，2023年，2月20日，星期四

解：令

f(1)<0,f(2)>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.6758二分法——例題

例1:用二分法求方程在區(qū)間(1,2)內(nèi)的實根,要求誤差限為。第53頁，共73頁，2023年，2月20日，星期四xyx*x0幾何意義1.5.2牛頓法迭代形式為第54頁，共73頁，2023年，2月20日，星期四[牛頓迭代法]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.牛頓法的算法構(gòu)造第55頁，共73頁，2023年，2月20日，星期四例1：利用牛頓迭代法求解f(x)=ex-1.5-tan-1x的零點。初始點x0=-7.0

第56頁，共73頁，2023年，2月20日，星期四解：f(x0)=-0.702×10-1，f'(x)=ex-(1+x2)-1

計算迭代格式:

計算結(jié)果如下表：(取|f(x)|<=10-10)kxf(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-012例1：利用牛頓迭代法求解f(x)=ex-1.5-tan-1x的零點。初始點x0=-7.0

第57頁，共73頁，2023年，2月20日，星期四注：Newton’sMethod收斂性依賴于x0

的選取。x*x0x0x0算法說明第58頁，共73頁，2023年，2月20日，星期四x0x1切線

割線

切線斜率割線斜率任意2個初值x0

和x1可以啟動這個遞推關(guān)系。1.5.3弦割法第59頁，共73頁，2023年，2月20日，星期四三種求根算法的比較第60頁，共73頁，2023年，2月20日，星期四二分法最穩(wěn)妥，但是效率最低。牛頓法效率最高，但是要求函數(shù)解析，容易計算導(dǎo)數(shù)。弦割法是前面兩種方法的折衷，既有較高的效率，又不必像牛頓法那樣必須計算函數(shù)f的導(dǎo)數(shù)。

如果初始猜測比較接近待求解，則其收斂速度幾乎與牛頓算法一樣快。三種求根算法的比較如果待求解附近，函數(shù)的行為不好，則自動的牛頓法和弦割法都可能無法收斂或收斂到錯誤的結(jié)果。保險的做法是先用二分法初步的定出解的位置，再用兩個自動方法中的一個定出解的精確位置。第61頁，共73頁，2023年，2月20日，星期四fzero(@(x)sin(x),[3,4]）Matlab自帶的求根指令X=fzero(函數(shù)句柄，猜測的初始值或搜尋的區(qū)間)fzero：求單變量函數(shù)的零點使用zeroin算法（結(jié)合了二分法、弦割法以及其它方法的一種綜合方法）。fzero(@(x)sin(x),3）第62頁，共73頁，2023年，2月20日，星期四1.6分子振動的半經(jīng)典量子化原子的相互作用勢為Lennard-Jones或6-12形式勢能最低處為rmin=21/6a，深度為V0第63頁，共73頁，2023年，2月20日，星期四能量為En的相對運動的振動態(tài)可以用一維薛定諤方程的束縛態(tài)解Ψn(r)來描述目標(biāo)：對給定的勢，求得能量En標(biāo)準(zhǔn)辦法：數(shù)值求解常微分方程本征問題我們這里采用的方法：半經(jīng)典量子化第64頁，共73頁，2023年，2月20日，星期四通過考慮原子核在勢場