《數(shù)值積分與常微分方程的數(shù)值解法.ppt》由會員分享，可在線閱讀，更多相關(guān)《數(shù)值積分與常微分方程的數(shù)值解法.ppt（40頁珍藏版）》請在裝配圖網(wǎng)上搜索。

● Euler法及其改進(jìn) ● Runge-Kutta法,● 梯形法 ● Simpson法 ●離散點數(shù)據(jù)的求積,第三章 數(shù)值積分與常微分方程的數(shù)值解法,,數(shù)值積分,常微分方程的數(shù)值解法,1.f(x)函數(shù)形式已知，但其積分不能表示成初等函數(shù)的閉合形式 2. f(x)函數(shù)形式未知，但其離散數(shù)據(jù)表已給出,,3-1-1 梯形法——方法原理,基本思想： 復(fù)化求積，即從近似計算為出發(fā)點，用有限項 的求和計算來代替從而求出定積分的近似值。,定步長：,求f(x)在[a,b]上的定積分,,y=f(x),,,a,b,xk-1,xk,h,,Ik,h——步長,變步長：,N個區(qū)間，h，T1,2N個區(qū)間，h/2，T2,,,|T2-T1|EPS,（xk-1,xk）,,xk-1/2,其中：,3-1-1 梯形法——方法原理,,,,,,,,,,例：Debye-Einstein公式推導(dǎo)得到計算固體熱容的公式為,其中：,?D為Debye溫度，R為氣體常數(shù)8.314JK-1mol-1,已知固體的Debye溫度如下：,求在50，100，298.15，500，1500K時，各固體的熱容。,3-2-1-1 Simpson法——問題的提出,,,求積分,Simpson法是把積分區(qū)間分割成有限個小區(qū)間，在每個小區(qū)間上采用二次拋物線來近似被積函數(shù)f(x)的圖形，近似求出小區(qū)間的面積，然后再將有限個小區(qū)間相加得到被積函數(shù)的近似值。,,y=f(x),,,xi-1,xi+1,,xi,y=g(x),h,h,Si,定步長:,3-2-1-2 Simpson法——方法原理,,變步長：,其中：,3-2-1-2 Simpson法——方法原理,,,,,,,,,判據(jù)：,3-2-1-2 Simpson法——方法原理,3-2-1-3 Simpson法——程序框圖,Simp(A,B,EPS,S2,F),,N=1, H=B-A, S1=0, T1=H*(F(A)+F(B))/2,DO K=1，N,,S=0,S=S+F(A+(K-1/2)*H),T2=(T1+H*S)/2, S2=T2+(T2-T1)/3, D=|S2-S1|,|S2|1,,D=|(S2-S1)/S2|,DEPS,,,,No,,RETURN,Yes,,No,N=N*2 H=H/2 T1=T2 S1=S2,,Yes,,,,,,,,,3-2-1-4 Simpson法——應(yīng)用示例,開始,,輸入： Debye溫度T（5），精度EPS，溫度THETA,,輸出：固體的熱容Cv,,結(jié)束,調(diào)用Simpson積分法子程序計算式右方積分值S2,,,計算：XM=THETA/T(I) (I=1,N),輸入：積分上下限A=10-4，B=XM,B=0,,,,,Yes,,No,固體的熱容Cv=9R/XM**3*S2,,顯示程序 顯示輸出,,3-1-3 –1 離散點數(shù)據(jù)的求積——方法原理,實驗時，得不到變量間的關(guān)系式，只測量到（xi,yi） 的離散點數(shù)據(jù)。,,,a,b,,方法 ：,1.用插值程序求任意 點的函數(shù)值。,一元三點Lagrange插值：,2．用Simpson求積程 序計算[a,b]區(qū)間中離 散點下的面積。,Simp(M,A,B,X,Y,EPS,S2),,N=1, H=B-A, （1）； S1=0, T1=H*(F(A)+F(B))/2,DO K=1，N,,S=0,（2）； S=S+F(A+(K-1/2)*H),T2=(T1+H*S)/2, S2=T2+(T2-T1)/3, D=|S2-S1|,|S2|1,,D=|(S2-S1)/S2|,DEPS,,,,No,,RETURN,Yes,,No,N=N*2 H=H/2 T1=T2 S1=S2,,Yes,,,,,,,,,(1)調(diào)用Lagrange一元三點插值F(A),F(B),(2)調(diào)用Lagrange一元三點插值F(A+(K-1/2)*H),3-1-3 –2 離散點數(shù)據(jù)的求積——程序框圖,純氣體的逸度由定義：,（1）,代入（1）,積分，并取極低壓力下氣體視為理想氣體，得,逸度：,φ為逸度系數(shù),（2）,例1： 實際氣體逸度的計算,已知p~Vm數(shù)據(jù),3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,開始,,輸入：數(shù)據(jù)點數(shù)N，精度EPS，溫度T 壓力p和摩爾體積Vm的實驗數(shù)據(jù)X(I),Y(I) （I=1，N）,,輸出：B, FI, FF,,結(jié)束,調(diào)用離散點求積子程序計算（2）式右方積分值S,,,計算：Y(I)=1/X(I)-Y(I)/RT (I=1,N),輸入：要計算的壓力P，積分上下限A=0，B=P,B=0,,,,,Y,,N,逸度系數(shù)FI=EXP(S),逸度FF=FI*B,,例2：已知固體Pb的熱容Cp~溫度T數(shù)據(jù)，求從15K到550K的固體Pb的焓變。,例3：分子標(biāo)準(zhǔn)熵S及Cp~T數(shù)據(jù)，求500K時的熵S值。,T1：298.15K T2:500K,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,已知數(shù)據(jù),,,,例4：合成氨反應(yīng),焓變?H與溫度T數(shù)據(jù)，已知623K下Kp1，求773K下Kp2。,,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,開始,,輸入：焓變與溫度的實驗數(shù)據(jù)X(I),Y(I) （I=1，N）,,輸出：KP2,,結(jié)束,調(diào)用離散點求積子程序計算積分值S,,,計算：XI=X(I),Y(I)=Y(I)/(XI*XI) (I=1,N),輸入：積分上、下限T2,T1及KP1,,計算KP2=KP1*EXP(S/R),,顯示程序 顯示輸出,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,,,,例5：已知：由A、B組成的二元混合物經(jīng)色譜分析 得到兩個分開的峰，時間和峰高的數(shù)據(jù)如下：,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,求A、B兩種物質(zhì)相對含量之比。,在色譜圖上，色譜峰的面積與色譜分析中各物質(zhì)的含量成正比。,開始,,輸入：A、B兩物質(zhì)的時間X與濃度峰高Y的實驗數(shù)據(jù) X(I),Y(I) X1(I),Y1(I),,輸出：S/S1,,結(jié)束,兩次調(diào)用離散點積分法子程序計算A，B物質(zhì)的峰面積S，S1 （其中調(diào)用Lagrange插值法子程序計算任意點的函數(shù)值）,,,輸入：A、B兩物質(zhì)和積分上、下限A，B A1，B1,計算A，B物質(zhì)的相對含量之比S/S1,,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,,例6：在簡單蒸餾釜內(nèi)蒸餾1000 Kg含C2H5OH質(zhì)量分?jǐn)?shù)為60%和H2O質(zhì)量分?jǐn)?shù)為40%的混合液。蒸餾結(jié)束時，殘液中含C2H5OH質(zhì)量分?jǐn)?shù)為5%。試求殘液的質(zhì)量是多少千克？該體系的汽液平衡數(shù)據(jù)如下：其中 x為液相中C2H5OH的質(zhì)量分?jǐn)?shù)，y 為汽相中C2H5OH的質(zhì)量分?jǐn)?shù)。,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,,,簡單蒸餾的雷利公式為：,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,,,式中，F(xiàn)為原料液量，W為殘液量； xF為原料液組成，xW為殘液組成。,開始,,輸入：原料液量F, x, 1/(y-x),,輸出：W,,結(jié)束,兩次調(diào)用離散點積分法子程序計算右側(cè)積分S （其中調(diào)用Lagrange插值法子程序計算任意點的函數(shù)值）,,,輸入：積分上、下限XW，XF,計算殘液量W＝F/exp(S),,顯示程序 顯示輸出,3-1-3 –3 離散點數(shù)據(jù)的求積——應(yīng)用示例,3-2 常微分方程的數(shù)值解法——引言,一階常微分方程的初值問題：,一階常微分方程組的初值問題：,數(shù)值解法：尋求解y(x)在一系列離散點 上的近似值 , 使y與x的關(guān)系近似滿足y=F(x) 步長： 假定h為定值,（1）,數(shù)值解的特點：,找一個遞推公式,（1）式積分,（2）,3-2 常微分方程的數(shù)值解法——引言,,,數(shù)值積分,3-2-1-1 Euler法及其改進(jìn)——問題的提出,例：異丙苯氧化反應(yīng)的動力學(xué),已知：t=0時的[RH]、[ROOH]，120℃時的反應(yīng)速率常數(shù)k， 計算0-14h每隔2h的[ROOH]。,解：,設(shè)x=[ROOH]，則[RH]=c -x,,右端利用矩形求積：得,——Euler公式 顯式,幾何意義：,取切線的端點Pi+1 作為yi+1,,,xi,xi+1,,y=y(x),,pi,pi+1,h,,,3-2-1-2 Euler法及其改進(jìn)——方法原理,,,h,提高精度：,（2）式右端利用梯形法求積：,隱式,改進(jìn)：預(yù)報~校正,預(yù)報：,校正：,編程：,,3-2-1-2 Euler法及其改進(jìn)——方法原理,,,,EULER(F,X,Y,N,H),,X0=X(1),,DO I=1，N-1,X(I+1)=X0+I*H YP=Y(I)+H*F(X(I),Y(I)) YC= Y(I)+H*F(X(I+1),YP) Y(I+1)=(YP+YC)/2,,,,結(jié)束,,3-2-1-3 Euler法及其改進(jìn)——程序框圖,開始,,輸入：c,k,t0,X0,t(I),,調(diào)用Euler法子程序計算X(I),,輸出：X(I),t(I),,結(jié)束,例1：異丙苯氧化反應(yīng)的動力學(xué),3-2-1-4 Euler法及其改進(jìn)——應(yīng)用示例,顯示程序 顯示輸出,例2：氣相色譜儀及其實驗過程的仿真,根據(jù)物料平衡的原理，可以得出塔板j上氣相物質(zhì)的 物料平衡方程式：,(1≤j≤N，N為塔板總數(shù)）,氣相各組分i（包括載氣）在塔板j 上的物料平衡方程式：,3-2-1-4 Euler法及其改進(jìn)——應(yīng)用示例,初始化,求各板壓力 Pj和氣相體積流量Vj,求各板氣相摩爾流量Gj和固定相物質(zhì)的量Lj,各組分物料衡算求出各板中各組分的滯料量(Euler法),求各組分氣、液相的摩爾分?jǐn)?shù)yj、 xj,用文件記錄xj、yj的值,,,,,求各板壓力，給各流股賦值,,,,,,輸入到繪圖軟件形成圖形,,,各 板 循 環(huán),模型計算 動態(tài)流程圖,,,,3-2-2-1 Runge-Kutta法——問題的提出,例：,基元反應(yīng)：,已知：[A]0=[B]0=1molL-1, [C]0=[D]0=[E]0=0 1molL-1, k1=1.0min-1， k2=0.25min-1， k3=0.5min-1， 求各組分濃度隨時間的變化（0-20min）.,基本思想：,求平均斜率,差商：,由微分中值定理，存在,使,得,由,,(xi,xi+1)上平均斜率,3-2-2-2 Runge-Kutta法——方法原理,二階RK公式：（改進(jìn)Euler）,改進(jìn)Euler： 二點斜率,,通式：,滿足：,,3-2-2-2 Runge-Kutta法——方法原理,——變形Euler,三階RK公式：（三點斜率）,或,3-2-2-2 Runge-Kutta法——方法原理,四階RK公式：（四點斜率）,或,3-2-2-2 Runge-Kutta法——方法原理,常微分方程組四階RK公式：,變步長：判據(jù)：,（方程）,（方程組）,3-2-2-2 Runge-Kutta法——方法原理,RK4(F,X,Y,N,H),,H2=H/2,H6=H/6,,DO I=1，N-1,XI=X(I), YI=Y(I) X12=XI+H2 , X11=XI+H X(I+1)=X11 RK1=F(XI,YI) RK2=F(X12,YI+H2*RK1) RK3=F(X12, YI+H2*RK2) RK4=F(X11,YI+H*RK3) Y(I+1)=YI+H6*(RK1+2(RK2+RK3)+RK4),,,,RETURN,,3-2-2-3 Runge-Kutta法——程序框圖,定步長：,RK41(F,X0,Y0,Y,H,K),X=X0,Y=Y0 H2=H/2,H6=H/6,DO I=1，K,X12=XI+H2 , X11=XI+H X(I+1)=X11 RK1=F(XI,YI) RK2=F(X12,YI+H2*RK1) RK3=F(X12, YI+H2*RK2) RK4=F(X11,YI+H*RK3) Y=Y+H6*(RK1+2(RK2+RK3)+RK4) X=X+H,變步長：,RETURN,,,,,,,方程組：,DO I=1,N (N為方程個數(shù)) CALL RK4( ),3-2-2-3 Runge-Kutta法——程序框圖,開 始,,賦值：Y0,N=5,TO=0,T1=20,H0=1.0,EPS=1E-6,H=0,M=(T1-T0)/H0,,DO I=1，M,K=1 CALL RK41(T0,Y0,Y1,H0,N,K),,,,K=K+K,H=H0/K,CALL K41(T0,Y0,Y2,H0,N,K) ES=0,,DO J=1，N,,,,ES=ES+ABS(Y2(J)-Y1(J))/Y2(J),,,N0,J=1,N YO(J)=Y2(J),,輸出T0,H,Y2(J),結(jié)束,,,,yes,J=1,N Y1(J)=Y2(J),,,ESEPS,顯示程序 顯示輸出,3-2-2-4 Runge-Kutta法——應(yīng)用示例,綜合應(yīng)用示例：,● 化學(xué)反應(yīng)動力學(xué)的計算與計算機(jī)模擬 ● 二組元汽-液相平衡的計算機(jī)模擬 ● 在固體催化劑表面上對氣體吸附等溫 方程式的計算及過程的計算機(jī)模擬,