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

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,其中：,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-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法应用示例,开始,输入： Debye温度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=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）,积分，并取极低压力下气体视为理想气体，得,逸度：,为逸度系数,（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：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组成的二元混合物经色谱分析 得到两个分开的峰，时间和峰高的数据如下：,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质量分数为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),显示程序 显示输出,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公式 显式,几何意义：,取切线的端点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法子程序计算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的值,求各板压力，给各流股赋值,输入到绘图软件形成图形,各 板 循 环,模型计算 动态流程图,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,三阶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+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,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),结束,yes,J=1,N Y1(J)=Y2(J),ESEPS,显示程序 显示输出,3-2-2-4 Runge-Kutta法应用示例,综合应用示例：, 化学反应动力学的计算与计算机模拟 二组元汽-液相平衡的计算机模拟 在固体催化剂表面上对气体吸附等温 方程式的计算及过程的计算机模拟,