計算方法第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法2 27.1引入—三論PageRank算法設(shè)網(wǎng)頁i地重要性為Pri,鏈出數(shù)為Li。如果網(wǎng)頁i存在一個指向網(wǎng)頁A地鏈接,則表明i地所有者認(rèn)為A比較重要,從而把i地一部分重要性得分賦予A。一個頁面地PageRank是由所有鏈向它地頁面（鏈入頁面）地重要性經(jīng)過遞歸算法得到。假設(shè)世界上只有四張網(wǎng)頁:A,B,C,D3 37.1引入—三論PageRank算法網(wǎng)頁間地鏈接矩陣M為(7.1.1)其中M地元素mij為0時表示第i個網(wǎng)頁沒有到第j個網(wǎng)頁地鏈接,否則mij為1/Li4 47.1引入—三論PageRank算法將網(wǎng)頁地重要性Pr看作是不斷變化地變量,則每次地變化為M·Pr-Pr=(M-I)·Pr。用建立下面地常微分方程地方法進(jìn)行求解。(7.1.2)5 57.1引入—三論PageRank算法分別取Pr地初值為[0.25,0.25,0.25,0.25]T [1,0,0,0]T6 67.1引入常微分方程少數(shù)簡單類型地常微分方程能求得精確解析有些常微分方程求解過程極為復(fù)雜多數(shù)情況只能使用近似解法求得其近似解常微分方程數(shù)值解法很有必要7 7第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法8 87.2常微分方程初值問題一階常微分方程地定解問題(7.2.1)其中:x∈[a,b]為自變量=(y1,y2,…,yd)∈Rd,y=y(x)為向量函數(shù),(x,y):RRd→Rd稱為右端向量場y0∈Rd稱為初值。當(dāng)給定向量場f(x,y)與初值y0,在區(qū)間[a,b]上求函數(shù)值y(x),使其滿足方程組(7.2.1)地問題稱為初值問題。在d=1地情況,方程組(7.2.1)為簡單地常微分方程初值問題。9 9常微分方程初值問題地解地存在性與唯一性定理定理7.1假設(shè)f(x,y)在區(qū)域G={(x,y)|a≤x≤b,|y|<∞內(nèi)連續(xù),并對y滿足,利普希茨(Lipschitz)條件,即存在常數(shù)L>0,使得(7.2.2)對所有x∈[a,b]與任意y1,y2∈R成立,則初值問題(7.2.1)在[a,b]上存在唯一解y(x),若f∈Ck(G),k≥0,那么其解y(x)∈Ck+1([a,b])107.2常微分方程初值問題例7.1馬爾薩斯模型:馬爾薩斯在分析人口出生與死亡情況地資料后發(fā)現(xiàn),人口凈增長率r基本上是一個常數(shù),提出了著名地人口指數(shù)增長模型。(7.2.3)其中,t表示時間,N表示人口數(shù)117.2常微分方程初值問題例7.2一重物垂直作用于彈簧所引起地振蕩,當(dāng)運動阻力與速度平方成正比時,可用一個二階常微分方程進(jìn)行描述:(7.2.4)12例7.2若令則上述二階常微分方程可以化為等價地一階常微分方程組(7.2.5)13第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法14 147.3歐拉方法及其改進(jìn)對于初值問題(7.2.1),在求解區(qū)間[a,b]上取等間距節(jié)點a=x0<x1<x2<…<xN-1<xN=b(7.3.1)其中?xn=xn+1-xn=h稱為積分網(wǎng)格地步長常微分方程(7.2.1)初值問題地數(shù)值解采用數(shù)值算法計算出初值問題精確解y(x)在節(jié)點x1,x2,…,xN上地函數(shù)值地近似值y1,y2,…,yN。常用地數(shù)值算法:幾何方法,數(shù)值微分,數(shù)值積分與Taylor展開15 157.3歐拉方法及其改進(jìn)如果在計算yn時該迭代公式只用到已經(jīng)求出地yn-1,而不使用y1,y2,…,yn-2中地任何一個則稱此算法為單步方法否則稱之為多步方法167.3.1歐拉方法17 177.3.1歐拉方法取等距節(jié)點組,此時有xn=a+nh,n=0,1,…,N,?xn≡h=(b-a)/N若以向前差商近似代替(7.2.1)中地導(dǎo)數(shù)y‘(xn),得(7.3.2)187.3.1歐拉方法將其看作等式,并記y(xn)地近似值為yn,則有(7.3.3)19歐拉方法Octave程序functiony=Euler(f,y0,a,b,N)%EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=b%withinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a,y0);fori=2:Ny(i)=y(i-1)+h*f(a+(i-1)*h,y(i-1));endendfunction207.3.2歐拉方法地幾何解釋21 217.3.2歐拉方法地幾何解釋前向歐拉方法地求解過程可以看作:從初值點(x0,y0)出發(fā),以斜率f(x0,y0)做切線在經(jīng)過步長h后與另一條積分曲線相交于(x1,y1)依次得到y(tǒng)1,y2,…,yN依次連接(x0,y0),(x1,y1),…,(xN,yN)地折線去近似方程y’(x)=f(x,y(x))過點(x0,y0)地積分曲線該方法又稱為折線法227.3.2歐拉方法地幾何解釋237.3.3歐拉方法地截斷誤差24 247.3.3歐拉方法地截斷誤差假設(shè)y(x)是y’(x)=f(x,y(x))任意一個二次連續(xù)可微地解,則將y(x)在點xn處進(jìn)行Taylor展開,可得(7.3.4)若取x=xn+1,則有(7.3.5)257.3.3歐拉方法地截斷誤差余項(7.3.6)為向前歐拉公式(7.3.3)地誤差稱為向前歐拉公式地截斷誤差267.3.4向后歐拉方法27 277.3.4向后歐拉方法用去近似導(dǎo)數(shù)有(7.3.7)記y(xn)地近似值為yn,則有(7.3.8)28"預(yù)報-校正"法先以向前歐拉公式給出yn+1地初值(7.3.9)然后再用迭代公式(7.3.10)計算yn+1地近似值。只用(7.3.10)式計算一次,而不進(jìn)行迭代29"預(yù)報-校正"法Octave代碼functiony=Euler(f,y0,a,b,N)EulermethodtosolveOrdinarydifferentialequationy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a,y0);fori=2:Ny(i)=y(i-1)+h*f(a+(i-1)*h,y(i-1));endendfunction30第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法31 317.4.1梯形方法常微分方程地定解問題(7.2.1)式中地將其在[xn,xn+1]區(qū)間上進(jìn)行積分,可得(7.4.1)如果等式右端地積分可以求出,就可以得到y(tǒng)(xn+1)可以采用不同地數(shù)值積分公式進(jìn)行近似求解327.4.1梯形方法用梯形求積公式,可得(7.4.2)其中(7.4.3)為梯形求積公式余項337.4.1梯形方法記y(xn)地近似值為yn,則有(7.4.4)稱此遞推公式為梯形公式347.4.2改進(jìn)歐拉格式35 357.4.2改進(jìn)歐拉格式例7.3分別用向前歐拉法,基于"預(yù)報-校正"技術(shù)地向后歐拉法與改進(jìn)歐拉格式求解初值問題36 36例7.3解:取步長h=0.1,則xn=0.1n(n=0,1,…,10),并由y0=1,分別采用三種方法進(jìn)行求解。向前歐拉法地迭代公式為:37 37例7.3解:基于"預(yù)報-校正"技術(shù)地向后歐拉法地迭代公式為:38 38例7.3解:改進(jìn)歐拉格式地迭代公式為:39 39例7.3地比較結(jié)果xn向前Euler法向后Euler法改進(jìn)Euler解析解11.10001.09181.09591.095421.19181.17631.18411.183231.27741.25461.26621.264941.35821.32781.34341.341651.43511.39641.41641.414261.50901.46091.48601.483271.58031.52161.55251.549281.64981.57861.61651.612591.71781.63211.67821.6733101.78481.68191.73791.732140 40例7.3地比較結(jié)果41 41第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法42 427.5.1龍格-庫塔方法地基本思想43 437.5.1龍格-庫塔方法地基本思想設(shè)y(x)是初值問題(7.2.1)地解,將y(xn+1)在xn處進(jìn)行泰勒展開,只保留常數(shù)項,有(7.5.1)其中f(ξ,y(ξ))為y(x)在區(qū)間[xn,xn+1]上地平均斜率用區(qū)間[xn,xn+1]內(nèi)更多節(jié)點地斜率地線性組合得到平均斜率地更高逼近44 447.5.2龍格-庫塔方法45 457.5.2龍格-庫塔方法設(shè)y(x)是初值問題(7.2.1)地解,并充分光滑,則將y(x)在xn處進(jìn)行泰勒展開,取x=xn+1,有(7.5.2)其中(7.5.3)為余項467.5.2龍格-庫塔方法將該余項略去,并記y(xn)地近似值為yn,則有(7.5.4)由微分地鏈?zhǔn)椒▌t可知(7.5.5)477.5.2龍格-庫塔方法將該余項略去,并記y(xn)地近似值為yn,則有(7.5.4)由微分地鏈?zhǔn)椒▌t可知(7.5.5)487.5.2龍格-庫塔方法用f(x,y)地函數(shù)值地線性組合近似(7.5.4)式右側(cè)y(x)地各階導(dǎo)數(shù)之與,即(7.5.6)其中(7.5.7)497.5.3二級龍格-庫塔格式50 507.5.3二級龍格-庫塔格式利用(7.5.5)式將(7.5.2)式改寫為(7.5.8)其中f,f’x,f’y分別代表f(x,y(x))及其對x與y地一階偏導(dǎo)數(shù)在(xn,yn)上取值。517.5.3二級龍格-庫塔格式下面分析近似公式(7.5.6)與(7.5.7),級數(shù)M與參數(shù)αm,λm,μm地取值都不是唯一地不同地取法可以得到不同地迭代格式設(shè)M=2,由二元函數(shù)地泰勒展開式可以得到(7.5.9)527.5.3二級龍格-庫塔格式將其代入(7.5.7)地第二式,進(jìn)而將(7.5.7)式代入(7.5.6)式,可得(7.5.10)比較(7.5.8)與(7.5.10)兩式中h與h2地系數(shù),并令對應(yīng)地各項系數(shù)相等,有(7.5.11)537.5.3二級龍格-庫塔格式令則得二級龍格-庫塔格式為(7.5.12)547.5.3二級龍格-庫塔格式令則得二級Runge-Kutta格式為(7.5.13)55二級龍格-庫塔格式Octave代碼functiony=midpoint_RuggeKutta(f,y0,a,b,N)midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N,1);y(1)=y0+h*f(a+h/2,y0+h*f(a,y0)/2);fori=2:Ny(i)=y(i-1)+h*(f(a+(i-1)*h+h/2,y(i-1))+h*f(a+(i-1)*h,y(i-1))/2);endendfunction567.5.4四級龍格-庫塔格式57 577.5.4四級龍格-庫塔格式當(dāng)p=4時,取M=4(7.5.16)587.5.4四級龍格-庫塔格式進(jìn)一步利用Taylor展開式與(7.5.5)式可以推出(7.5.17)其中,是把(7.5.16)式ki中地yn與ki-1用y(xn)與代替后得到地表達(dá)式。并有(7.5.18)四級龍格-庫塔格式地截斷誤差為h地5階無窮小量59四級龍格-庫塔格式Octave程序functiony=Fourth_RuggeKutta(f,y0,a,b,N)midpoint_RuggeKuttamethodtosolveODEy'=f(x,y),a<x<=bwithinitialconditiony(a)=y0.[a,b]willbedividedbyNsegmentsevenly.h=(b-a)/N;y=zeros(N+1,1);y(1)=y0;fori=2:N+1k1=f(a+(i-2)*h,y(i-1));k2=f(a+(i-2)*h+h/2,y(i-1)+k1*h/2);k3=f(a+(i-2)*h+h/2,y(i-1)+k2*h/2);k4=f(a+(i-1)*h,y(i-1)+k3*h);y(i)=y(i-1)+h*(k1+2*k2+2*k3+k4)/6;endendfunction607.5.4四級龍格-庫塔格式例7.4用四階龍格-庫塔格式求解下面初值問題并與例7.3結(jié)果進(jìn)行比較61例7.4解:取步長h=0.1,四階龍格-庫塔格式為62例7.4地結(jié)果比較xn向前Euler法向后Euler法改進(jìn)Euler四階Runge-解析解Kutta11.10001.09181.09591.09541.095421.19181.17631.18411.18321.183231.27741.25461.26621.26491.264941.35821.32781.34341.34161.341651.43511.39641.41641.41421.414261.50901.46091.48601.48321.483271.58031.52161.55251.54921.549281.64981.57861.61651.61251.612591.71781.63211.67821.67331.6733101.78481.68191.73791.73211.732163第七章常微分方程初值問題地數(shù)值解7.1引入7.2常微分方程初值問題7.3歐拉方法及其改進(jìn)7.4梯形方法7.5龍格-庫塔方法7.6常微分方程組地數(shù)值解法64 647.6常微分方程組地數(shù)值解法65 657.6常微分方程組地數(shù)值解法一階常微分方程組(7.6.1)其向前歐拉格式為(7.6.2)667.6常微分方程組地數(shù)值解法改進(jìn)地歐拉格式為(7.6.3)四階龍格-庫塔格式為(7.6.4)677.6常微分方程組地數(shù)值解法例7.5采用向前歐拉格式,