数值积分与常微分方程的数值解法

1、Euler法及其改进 Runge-Kutta法,梯形法 Simpson法 离散点数据的求积,第三章 数值积分与常微分方程的数值解法,数值积分,常微分方程的数值解法,1.f(x)函数形式已知，但其积分不能表示成初等函数的闭合形式 2. f(x)函数形式未知，但其离散数据表已给出,3-1-1 梯形法方法原理,基本思想： 复化求积，即从近似计算为出发点，用有限项 的求和计算来代替从而求出定积分的近似值,定步长,求f(x)在a,b上的定积分,y=f(x,a,b,xk-1,xk,h,Ik,h步长,变步长,N个区间，h，T1,2N个区间，h/2，T2,T2-T1|EPS,xk-1,xk,xk-1/2,其中

2、,3-1-1 梯形法方法原理,例：Debye-Einstein公式推导得到计算固体热容的公式为,其中,D为Debye温度，R为气体常数8.314JK-1mol-1,已知固体的Debye温度如下,求在50，100，298.15，500，1500K时，各固体的热容,3-2-1-1 Simpson法问题的提出,求积分,Simpson法是把积分区间分割成有限个小区间，在每个小区间上采用二次抛物线来近似被积函数f(x)的图形，近似求出小区间的面积，然后再将有限个小区间相加得到被积函数的近似值,y=f(x,xi-1,xi+1,xi,y=g(x,h,h,Si,定步长,3-2-1-2 Simpson法方法原理

3、,变步长,其中,3-2-1-2 Simpson法方法原理,判据,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法应用示例,开始,输入： Deb

4、ye温度T（5），精度EPS，温度THETA,输出：固体的热容Cv,结束,调用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 离散点数据的求积方法原理,实验时，得不到变量间的关系式，只测量到（xi,yi） 的离散点数据,a,b,方法,1.用插值程序求任意 点的函数值,一元三点Lagrange插值,2用Simpson求积程 序计算a,b区间中离 散点下的面积,Simp(M,A,B,X,Y,EPS,S2,N=1, H

5、=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)调用Lagrange一元三点插值F(A),F(B,2)调用Lagrange一元三点插值F(A+(K-1/2)*H,3-1-3 2 离散点数据的求积程序框图,纯气体的逸度由定义,1,代入（1,积分，并取极低压力下气体视为理想气体，得,逸度

6、,为逸度系数,2,例1： 实际气体逸度的计算,已知pVm数据,3-1-3 3 离散点数据的求积应用示例,开始,输入：数据点数N，精度EPS，温度T 压力p和摩尔体积Vm的实验数据X(I),Y(I) （I=1，N,输出：B, FI, FF,结束,调用离散点求积子程序计算（2）式右方积分值S,计算：Y(I)=1/X(I)-Y(I)/RT (I=1,N,输入：要计算的压力P，积分上下限A=0，B=P,B=0,Y,N,逸度系数FI=EXP(S),逸度FF=FI*B,例2：已知固体Pb的热容Cp温度T数据，求从15K到550K的固体Pb的焓变,例3：分子标准熵S及CpT数据，求500K时的熵S值,T1：

7、298.15K T2:500K,3-1-3 3 离散点数据的求积应用示例,已知数据,例4：合成氨反应,焓变H与温度T数据，已知623K下Kp1，求773K下Kp2,3-1-3 3 离散点数据的求积应用示例,开始,输入：焓变与温度的实验数据X(I),Y(I) （I=1，N,输出：KP2,结束,调用离散点求积子程序计算积分值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 离散点数据的求积应用示例,例5：已知：由A、B组成的二元混合物经色谱分析 得到两个分开的峰，

8、时间和峰高的数据如下,3-1-3 3 离散点数据的求积应用示例,求A、B两种物质相对含量之比,在色谱图上，色谱峰的面积与色谱分析中各物质的含量成正比,开始,输入：A、B两物质的时间X与浓度峰高Y的实验数据 X(I),Y(I) X1(I),Y1(I,输出：S/S1,结束,两次调用离散点积分法子程序计算A，B物质的峰面积S，S1 （其中调用Lagrange插值法子程序计算任意点的函数值,输入：A、B两物质和积分上、下限A，B A1，B1,计算A，B物质的相对含量之比S/S1,3-1-3 3 离散点数据的求积应用示例,例6：在简单蒸馏釜内蒸馏1000 Kg含C2H5OH质量分数为60%和H2O质量分

9、数为40%的混合液。蒸馏结束时，残液中含C2H5OH质量分数为5%。试求残液的质量是多少千克？该体系的汽液平衡数据如下：其中 x为液相中C2H5OH的质量分数，y 为汽相中C2H5OH的质量分数,3-1-3 3 离散点数据的求积应用示例,简单蒸馏的雷利公式为,3-1-3 3 离散点数据的求积应用示例,式中，F为原料液量，W为残液量； xF为原料液组成，xW为残液组成,开始,输入：原料液量F, x, 1/(y-x,输出：W,结束,两次调用离散点积分法子程序计算右侧积分S （其中调用Lagrange插值法子程序计算任意点的函数值,输入：积分上、下限XW，XF,计算残液量WF/exp(S,显示程序

10、显示输出,3-1-3 3 离散点数据的求积应用示例,3-2 常微分方程的数值解法引言,一阶常微分方程的初值问题,一阶常微分方程组的初值问题,数值解法：寻求解y(x)在一系列离散点 上的近似值 , 使y与x的关系近似满足y=F(x) 步长： 假定h为定值,1,数值解的特点,找一个递推公式,1）式积分,2,3-2 常微分方程的数值解法引言,数值积分,3-2-1-1 Euler法及其改进问题的提出,例：异丙苯氧化反应的动力学,已知：t=0时的RH、ROOH，120时的反应速率常数k， 计算0-14h每隔2h的ROOH,解,设x=ROOH，则RH=c -x,右端利用矩形求积：得,Euler公式 显式,

11、几何意义,取切线的端点Pi+1 作为yi+1,xi,xi+1,y=y(x,pi,pi+1,h,3-2-1-2 Euler法及其改进方法原理,h,提高精度,2）式右端利用梯形法求积,隐式,改进：预报校正,预报,校正,编程,3-2-1-2 Euler法及其改进方法原理,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,结束,3-2-1-3 Euler法及其改进程序框图,开始,输入：c,k,t0,X0,t(I,调用Euler法子程

12、序计算X(I,输出：X(I),t(I,结束,例1：异丙苯氧化反应的动力学,3-2-1-4 Euler法及其改进应用示例,显示程序 显示输出,例2：气相色谱仪及其实验过程的仿真,根据物料平衡的原理，可以得出塔板j上气相物质的 物料平衡方程式,1jN，N为塔板总数,气相各组分i（包括载气）在塔板j 上的物料平衡方程式,3-2-1-4 Euler法及其改进应用示例,初始化,求各板压力 Pj和气相体积流量Vj,求各板气相摩尔流量Gj和固定相物质的量Lj,各组分物料衡算求出各板中各组分的滞料量(Euler法,求各组分气、液相的摩尔分数yj、 xj,用文件记录xj、yj的值,求各板压力，给各流股赋值,输入

13、到绘图软件形成图形,各 板 循 环,模型计算 动态流程图,3-2-2-1 Runge-Kutta法问题的提出,例,基元反应,已知：A0=B0=1molL-1, C0=D0=E0=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公式：（改进Euler,改进Euler： 二点斜率,通式,满足,3-2-2-2 Runge-Kutta法方法原理,变形Euler,

14、三阶RK公式：（三点斜率,或,3-2-2-2 Runge-Kutta法方法原理,四阶RK公式：（四点斜率,或,3-2-2-2 Runge-Kutta法方法原理,常微分方程组四阶RK公式,变步长：判据,方程,方程组,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+

15、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为方程个数) 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,