數(shù)值積分與常微分方程的數(shù)值解法.ppt

Euler法及其改進(jìn) Runge-Kutta法, 梯形法 Simpson法 離散點(diǎn)數(shù)據(jù)的求積,第三章 數(shù)值積分與常微分方程的數(shù)值解法,數(shù)值積分,常微分方程的數(shù)值解法,1.f(x)函數(shù)形式已知，但其積分不能表示成初等函數(shù)的閉合形式 2. f(x)函數(shù)形式未知，但其離散數(shù)據(jù)表已給出,3-1-1 梯形法方法原理,基本思想： 復(fù)化求積，即從近似計(jì)算為出發(fā)點(diǎn)，用有限項(xiàng) 的求和計(jì)算來(lái)代替從而求出定積分的近似值。,定步長(zhǎng)：,求f(x)在a,b上的定積分,y=f(x),a,b,xk-1,xk,h,Ik,h步長(zhǎng),變步長(zhǎng)：,N個(gè)區(qū)間，h，T1,2N個(gè)區(qū)間，h/2，T2,|T2-T1|EPS,（xk-1,xk）,xk-1/2,其中：,3-1-1 梯形法方法原理,例：Debye-Einstein公式推導(dǎo)得到計(jì)算固體熱容的公式為,其中：,D為Debye溫度，R為氣體常數(shù)8.314JK-1mol-1,已知固體的Debye溫度如下：,求在50，100，298.15，500，1500K時(shí)，各固體的熱容。,3-2-1-1 Simpson法問(wèn)題的提出,求積分,Simpson法是把積分區(qū)間分割成有限個(gè)小區(qū)間，在每個(gè)小區(qū)間上采用二次拋物線來(lái)近似被積函數(shù)f(x)的圖形，近似求出小區(qū)間的面積，然后再將有限個(gè)小區(qū)間相加得到被積函數(shù)的近似值。,y=f(x),xi-1,xi+1,xi,y=g(x),h,h,Si,定步長(zhǎng):,3-2-1-2 Simpson法方法原理,變步長(zhǎng)：,其中：,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)用示例,開(kāi)始,輸入： Debye溫度T（5），精度EPS，溫度THETA,輸出：固體的熱容Cv,結(jié)束,調(diào)用Simpson積分法子程序計(jì)算式右方積分值S2,計(jì)算：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 離散點(diǎn)數(shù)據(jù)的求積方法原理,實(shí)驗(yàn)時(shí)，得不到變量間的關(guān)系式，只測(cè)量到（xi,yi） 的離散點(diǎn)數(shù)據(jù)。,a,b,方法 ：,1.用插值程序求任意 點(diǎn)的函數(shù)值。,一元三點(diǎn)Lagrange插值：,2用Simpson求積程 序計(jì)算a,b區(qū)間中離 散點(diǎn)下的面積。,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一元三點(diǎn)插值F(A),F(B),(2)調(diào)用Lagrange一元三點(diǎn)插值F(A+(K-1/2)*H),3-1-3 2 離散點(diǎn)數(shù)據(jù)的求積程序框圖,純氣體的逸度由定義：,（1）,代入（1）,積分，并取極低壓力下氣體視為理想氣體，得,逸度：,為逸度系數(shù),（2）,例1： 實(shí)際氣體逸度的計(jì)算,已知pVm數(shù)據(jù),3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,開(kāi)始,輸入：數(shù)據(jù)點(diǎn)數(shù)N，精度EPS，溫度T 壓力p和摩爾體積Vm的實(shí)驗(yàn)數(shù)據(jù)X(I),Y(I) （I=1，N）,輸出：B, FI, FF,結(jié)束,調(diào)用離散點(diǎn)求積子程序計(jì)算（2）式右方積分值S,計(jì)算：Y(I)=1/X(I)-Y(I)/RT (I=1,N),輸入：要計(jì)算的壓力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及CpT數(shù)據(jù)，求500K時(shí)的熵S值。,T1：298.15K T2:500K,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,已知數(shù)據(jù),例4：合成氨反應(yīng),焓變H與溫度T數(shù)據(jù)，已知623K下Kp1，求773K下Kp2。,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,開(kāi)始,輸入：焓變與溫度的實(shí)驗(yàn)數(shù)據(jù)X(I),Y(I) （I=1，N）,輸出：KP2,結(jié)束,調(diào)用離散點(diǎn)求積子程序計(jì)算積分值S,計(jì)算：XI=X(I),Y(I)=Y(I)/(XI*XI) (I=1,N),輸入：積分上、下限T2,T1及KP1,計(jì)算KP2=KP1*EXP(S/R),顯示程序 顯示輸出,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,例5：已知：由A、B組成的二元混合物經(jīng)色譜分析 得到兩個(gè)分開(kāi)的峰，時(shí)間和峰高的數(shù)據(jù)如下：,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,求A、B兩種物質(zhì)相對(duì)含量之比。,在色譜圖上，色譜峰的面積與色譜分析中各物質(zhì)的含量成正比。,開(kāi)始,輸入：A、B兩物質(zhì)的時(shí)間X與濃度峰高Y的實(shí)驗(yàn)數(shù)據(jù) X(I),Y(I) X1(I),Y1(I),輸出：S/S1,結(jié)束,兩次調(diào)用離散點(diǎn)積分法子程序計(jì)算A，B物質(zhì)的峰面積S，S1 （其中調(diào)用Lagrange插值法子程序計(jì)算任意點(diǎn)的函數(shù)值）,輸入：A、B兩物質(zhì)和積分上、下限A，B A1，B1,計(jì)算A，B物質(zhì)的相對(duì)含量之比S/S1,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,例6：在簡(jiǎn)單蒸餾釜內(nèi)蒸餾1000 Kg含C2H5OH質(zhì)量分?jǐn)?shù)為60%和H2O質(zhì)量分?jǐn)?shù)為40%的混合液。蒸餾結(jié)束時(shí)，殘液中含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 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,簡(jiǎn)單蒸餾的雷利公式為：,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,式中，F(xiàn)為原料液量，W為殘液量； xF為原料液組成，xW為殘液組成。,開(kāi)始,輸入：原料液量F, x, 1/(y-x),輸出：W,結(jié)束,兩次調(diào)用離散點(diǎn)積分法子程序計(jì)算右側(cè)積分S （其中調(diào)用Lagrange插值法子程序計(jì)算任意點(diǎn)的函數(shù)值）,輸入：積分上、下限XW，XF,計(jì)算殘液量WF/exp(S),顯示程序 顯示輸出,3-1-3 3 離散點(diǎn)數(shù)據(jù)的求積應(yīng)用示例,3-2 常微分方程的數(shù)值解法引言,一階常微分方程的初值問(wèn)題：,一階常微分方程組的初值問(wèn)題：,數(shù)值解法：尋求解y(x)在一系列離散點(diǎn) 上的近似值 , 使y與x的關(guān)系近似滿足y=F(x) 步長(zhǎng)： 假定h為定值,（1）,數(shù)值解的特點(diǎn)：,找一個(gè)遞推公式,（1）式積分,（2）,3-2 常微分方程的數(shù)值解法引言,數(shù)值積分,3-2-1-1 Euler法及其改進(jìn)問(wèn)題的提出,例：異丙苯氧化反應(yīng)的動(dòng)力學(xué),已知：t=0時(shí)的RH、ROOH，120時(shí)的反應(yīng)速率常數(shù)k， 計(jì)算0-14h每隔2h的ROOH。,解：,設(shè)x=ROOH，則RH=c -x,右端利用矩形求積：得,Euler公式 顯式,幾何意義：,取切線的端點(diǎn)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ù)報(bào)校正,預(yù)報(bào)：,校正：,編程：,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)程序框圖,開(kāi)始,輸入：c,k,t0,X0,t(I),調(diào)用Euler法子程序計(jì)算X(I),輸出：X(I),t(I),結(jié)束,例1：異丙苯氧化反應(yīng)的動(dòng)力學(xué),3-2-1-4 Euler法及其改進(jìn)應(yīng)用示例,顯示程序 顯示輸出,例2：氣相色譜儀及其實(shí)驗(yàn)過(guò)程的仿真,根據(jù)物料平衡的原理，可以得出塔板j上氣相物質(zhì)的 物料平衡方程式：,(1jN，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),模型計(jì)算 動(dòng)態(tài)流程圖,3-2-2-1 Runge-Kutta法問(wèn)題的提出,例：,基元反應(yīng)：,已知：A0=B0=1molL-1, C0=D0=E0=0 1molL-1, k1=1.0min-1， k2=0.25min-1， k3=0.5min-1， 求各組分濃度隨時(shí)間的變化（0-20min）.,基本思想：,求平均斜率,差商：,由微分中值定理，存在,使,得,由,(xi,xi+1)上平均斜率,3-2-2-2 Runge-Kutta法方法原理,二階RK公式：（改進(jìn)Euler）,改進(jìn)Euler： 二點(diǎn)斜率,通式：,滿足：,3-2-2-2 Runge-Kutta法方法原理,變形Euler,三階RK公式：（三點(diǎn)斜率）,或,3-2-2-2 Runge-Kutta法方法原理,四階RK公式：（四點(diǎn)斜率）,或,3-2-2-2 Runge-Kutta法方法原理,常微分方程組四階RK公式：,變步長(zhǎng)：判據(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法程序框圖,定步長(zhǎng)：,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,變步長(zhǎng)：,RETURN,方程組：,DO I=1,N (N為方程個(gè)數(shù)) CALL RK4( ),3-2-2-3 Runge-Kutta法程序框圖,開(kāi) 始,賦值：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)動(dòng)力學(xué)的計(jì)算與計(jì)算機(jī)模擬 二組元汽-液相平衡的計(jì)算機(jī)模擬 在固體催化劑表面上對(duì)氣體吸附等溫 方程式的計(jì)算及過(guò)程的計(jì)算機(jī)模擬,