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

第三章数值积分与常微分方程的数值解法第1页，课件共40页，创作于2023年2月3-1-1梯形法——方法原理基本思想：复化求积，即从近似计算为出发点，用有限项的求和计算来代替从而求出定积分的近似值。定步长：求f(x)在[a,b]上的定积分xyy=f(x)abxk-1xkhIkh——步长第2页，课件共40页，创作于2023年2月变步长：N个区间，h，T12N个区间，h/2，T2|T2-T1|<EPS（xk-1,xk）xk-1/2其中：3-1-1梯形法——方法原理第3页，课件共40页，创作于2023年2月例：Debye-Einstein公式推导得到计算固体热容的公式为其中：

D为Debye温度，R为气体常数8.314JK-1mol-1已知固体的Debye温度如下：PbAgCuAlFeKClNaClC

D/*K882153153984202272811910求在50，100，298.15，500，1500K时，各固体的热容。3-2-1-1Simpson法——问题的提出第4页，课件共40页，创作于2023年2月求积分

Simpson法是把积分区间分割成有限个小区间，在每个小区间上采用二次抛物线来近似被积函数f(x)的图形，近似求出小区间的面积，然后再将有限个小区间相加得到被积函数的近似值。xyy=f(x)xi-1xi+1xiy=g(x)hhSi定步长:

3-2-1-2Simpson法——方法原理第5页，课件共40页，创作于2023年2月变步长：

其中：3-2-1-2Simpson法——方法原理第6页，课件共40页，创作于2023年2月判据：3-2-1-2Simpson法——方法原理第7页，课件共40页，创作于2023年2月3-2-1-3Simpson法——程序框图Simp(A,B,EPS,S2,F)N=1,H=B-A,S1=0,T1=H*(F(A)+F(B))/2DOK=1，NS=0S=S+F(A+(K-1/2)*H)T2=(T1+H*S)/2,S2=T2+(T2-T1)/3,D=|S2-S1||S2|>1D=|(S2-S1)/S2|D<EPSNoRETURNYesNoN=N*2H=H/2T1=T2S1=S2Yes第8页，课件共40页，创作于2023年2月3-2-1-4Simpson法——应用示例开始输入：Debye温度T（5），精度EPS，温度THETA输出：固体的热容Cv结束调用Simpson积分法子程序计算式右方积分值S2计算：XM=THETA/T(I)(I=1,N)输入：积分上下限A=10-4，B=XMB=0YesNo固体的热容Cv=9R/XM**3*S2显示程序显示输出第9页，课件共40页，创作于2023年2月3-1-3–1离散点数据的求积——方法原理实验时，得不到变量间的关系式，只测量到（xi,yi）的离散点数据。

xyab方法：1.用插值程序求任意点的函数值。一元三点Lagrange插值：

2．用Simpson求积程序计算[a,b]区间中离散点下的面积。

第10页，课件共40页，创作于2023年2月Simp(M,A,B,X,Y,EPS,S2)N=1,H=B-A,（1）；S1=0,T1=H*(F(A)+F(B))/2DOK=1，NS=0（2）；S=S+F(A+(K-1/2)*H)T2=(T1+H*S)/2,S2=T2+(T2-T1)/3,D=|S2-S1||S2|>1D=|(S2-S1)/S2|D<EPSNoRETURNYesNoN=N*2H=H/2T1=T2S1=S2Yes(1)调用Lagrange一元三点插值F(A),F(B)(2)调用Lagrange一元三点插值F(A+(K-1/2)*H)3-1-3–2离散点数据的求积——程序框图第11页，课件共40页，创作于2023年2月纯气体的逸度由定义：（1）代入（1）积分，并取极低压力下气体视为理想气体，得逸度：φ为逸度系数（2）例1：实际气体逸度的计算已知p~Vm数据3-1-3–3离散点数据的求积——应用示例第12页，课件共40页，创作于2023年2月开始输入：数据点数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=PB=0YN逸度系数FI=EXP(S),逸度FF=FI*B第13页，课件共40页，创作于2023年2月例2：已知固体Pb的热容Cp~温度T数据，求从15K到550K的固体Pb的焓变。例3：分子标准熵S及Cp~T数据，求500K时的熵S值。

T1：298.15KT2:500K3-1-3–3离散点数据的求积——应用示例第14页，课件共40页，创作于2023年2月

已知数据

例4：合成氨反应焓变

H与温度T数据，已知623K下Kp1，求773K下Kp2。T/K623.0648.0637.0698.0723.0748.0773.0

H/kJ·mol-1-50.7879-51.1390-51.4738-51.7943-52.1006-52.3929-52.6715

3-1-3–3离散点数据的求积——应用示例第15页，课件共40页，创作于2023年2月开始输入：焓变与温度的实验数据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离散点数据的求积——应用示例第16页，课件共40页，创作于2023年2月

例5：已知：由A、B组成的二元混合物经色谱分析得到两个分开的峰，时间和峰高的数据如下：3-1-3–3离散点数据的求积——应用示例求A、B两种物质相对含量之比。

在色谱图上，色谱峰的面积与色谱分析中各物质的含量成正比。

时间t/sec峰高/mm时间t/sec峰高/mm38013044032500.5531725866440266147569274218690282139768290510843298011919第17页，课件共40页，创作于2023年2月开始输入：A、B两物质的时间X与浓度峰高Y的实验数据

X(I),Y(I)X1(I),Y1(I)输出：S/S1结束两次调用离散点积分法子程序计算A，B物质的峰面积S，S1（其中调用Lagrange插值法子程序计算任意点的函数值）

输入：A、B两物质和积分上、下限A，BA1，B1计算A，B物质的相对含量之比S/S1

3-1-3–3离散点数据的求积——应用示例第18页，课件共40页，创作于2023年2月

例6：在简单蒸馏釜内蒸馏1000Kg含C2H5OH质量分数为60%和H2O质量分数为40%的混合液。蒸馏结束时，残液中含C2H5OH质量分数为5%。试求残液的质量是多少千克？该体系的汽液平衡数据如下：其中x为液相中C2H5OH的质量分数，y为汽相中C2H5OH的质量分数。3-1-3–3离散点数据的求积——应用示例xx0.0250.050.100.150.200.250.35.003.222.402.222.202.272.440.350.400.450.500.550.600.652.642.903.293.74403850296.66第19页，课件共40页，创作于2023年2月简单蒸馏的雷利公式为：3-1-3–3离散点数据的求积——应用示例式中，F为原料液量，W为残液量；

xF为原料液组成，xW为残液组成。

第20页，课件共40页，创作于2023年2月开始输入：原料液量F,x,1/(y-x)输出：W结束两次调用离散点积分法子程序计算右侧积分S（其中调用Lagrange插值法子程序计算任意点的函数值）

输入：积分上、下限XW，XF计算残液量W＝F/exp(S)

显示程序显示输出3-1-3–3离散点数据的求积——应用示例第21页，课件共40页，创作于2023年2月3-2常微分方程的数值解法——引言一阶常微分方程的初值问题：

一阶常微分方程组的初值问题：数值解法：寻求解y(x)在一系列离散点上的近似值,使y与x的关系近似满足y=F(x)

步长：假定h为定值（1）第22页，课件共40页，创作于2023年2月数值解的特点：找一个递推公式（1）式积分

（2）3-2常微分方程的数值解法——引言数值积分第23页，课件共40页，创作于2023年2月3-2-1-1Euler法及其改进——问题的提出例：异丙苯氧化反应的动力学已知：t=0时的[RH]、[ROOH]，120℃时的反应速率常数k，计算0-14h每隔2h的[ROOH]。解：

设x=[ROOH]，则[RH]=c-x第24页，课件共40页，创作于2023年2月右端利用矩形求积：得

——Euler公式显式几何意义：

取切线的端点Pi+1作为yi+1xyxixi+1y=y(x)pipi+1h3-2-1-2Euler法及其改进——方法原理h第25页，课件共40页，创作于2023年2月提高精度：（2）式右端利用梯形法求积：隐式

改进：预报~校正预报：

校正：

编程：

3-2-1-2Euler法及其改进——方法原理第26页，课件共40页，创作于2023年2月EULER(F,X,Y,N,H)X0=X(1)DOI=1，N-1X(I+1)=X0+I*HYP=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-3Euler法及其改进——程序框图第27页，课件共40页，创作于2023年2月开始输入：c,k,t0,X0,t(I)调用Euler法子程序计算X(I)输出：X(I),t(I)结束例1：异丙苯氧化反应的动力学3-2-1-4Euler法及其改进——应用示例显示程序显示输出第28页，课件共40页，创作于2023年2月例2：气相色谱仪及其实验过程的仿真根据物料平衡的原理，可以得出塔板j上气相物质的物料平衡方程式：

(1≤j≤N，N为塔板总数）

塔板j

示意图固定相流动相Gj-1Gjyi,jLjxi,j气相各组分i（包括载气）在塔板j上的物料平衡方程式：

3-2-1-4Euler法及其改进——应用示例第29页，课件共40页，创作于2023年2月初始化求各板压力Pj和气相体积流量Vj求各板气相摩尔流量Gj和固定相物质的量Lj各组分物料衡算求出各板中各组分的滞料量(Euler法)求各组分气、液相的摩尔分数yj、xj用文件记录xj、yj的值求各板压力，给各流股赋值输入到绘图软件形成图形各板循环模型计算动态流程图第30页，课件共40页，创作于2023年2月3-2-2-1Runge-Kutta法——问题的提出例：

基元反应：

k2已知：[A]0=[B]0=1mol·L-1,[C]0=[D]0=[E]0=01mol·L-1,k1=1.0min-1，k2=0.25min-1，k3=0.5min-1，求各组分浓度随时间的变化（0-20min）.第31页，课件共40页，创作于2023年2月基本思想：

求平均斜率

差商：

由微分中值定理，存在

使

得由(xi,xi+1)上平均斜率3-2-2-2Runge-Kutta法——方法原理第32页，课件共40页，创作于2023年2月二阶RK公式：（改进Euler）改进Euler：二点斜率通式：满足：3-2-2-2Runge-Kutta法——方法原理第33页，课件共40页，创作于2023年2月——变形Euler三阶RK公式：（三点斜率）或3-2-2-2Runge-Kutta法——方法原理第34页，课件共40页，创作于2023年2月四阶RK公式：（四点斜率）或3-2-2-2Runge-Kutta法——方法原理第35页，课件共40页，创作于2023年2月常微分方程组四阶RK公式：

变步长：判据：

（方程）（方程组）3-2-2-2Runge-Kutta法——方法原理第36页，课件共40页，创作于2023年2月RK4(F,X,Y,N,H)H2=H/2,H6=H/6DOI=1，N-1XI=X(I),YI=Y(I)X12=XI+H2,X11=XI+HX(I+1)=X11RK1=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)RETURN3-2-2-3Runge-Kutta法——程序框图定步长：第37页，课件共40页，创作于2023年2月R