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

《数值积分与常微分方程的数值解法.ppt》由会员分享，可在线阅读，更多相关《数值积分与常微分方程的数值解法.ppt（40页珍藏版）》请在装配图网上搜索。

● 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： 实际气体逸度的计算,已知p~Vm数据,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及Cp~T数据，求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,计算残液量W＝F/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上气相物质的 物料平衡方程式：,(1≤j≤N，N为塔板总数）,气相各组分i（包括载气）在塔板j 上的物料平衡方程式：,3-2-1-4 Euler法及其改进——应用示例,初始化,求各板压力 Pj和气相体积流量Vj,求各板气相摩尔流量Gj和固定相物质的量Lj,各组分物料衡算求出各板中各组分的滞料量(Euler法),求各组分气、液相的摩尔分数yj、 xj,用文件记录xj、yj的值,,,,,求各板压力，给各流股赋值,,,,,,输入到绘图软件形成图形,,,各 板 循 环,模型计算 动态流程图,,,,3-2-2-1 Runge-Kutta法——问题的提出,例：,基元反应：,已知：[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公式：（改进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法——应用示例,综合应用示例：,● 化学反应动力学的计算与计算机模拟 ● 二组元汽-液相平衡的计算机模拟 ● 在固体催化剂表面上对气体吸附等温 方程式的计算及过程的计算机模拟,