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

●Euler法及其改进●Runge-Kutta法●梯形法●Simpson法●离散点数据的求积第三章数值积分与常微分方程的数值解法数值积分常微分方程的数值解法1.f(x)函数形式已知，但其积分不能表示成初等函数的闭合形式2.f(x)函数形式未知，但其离散数据表已给出3-1-1梯形法——方法原理基本思想：复化求积，即从近似计算为出发点，用有限项的求和计算来代替从而求出定积分的近似值。定步长：),,2,1(,/)(,nknabhkhaxk求f(x)在[a,b]上的定积分xyy=f(x)abxk-1xkhIk11111)()]()([2)]()([21nkknkkknkknxfhbfafhxfxfhITbadxxf)(h——步长变步长：N个区间，h，T12N个区间，h/2，T2|T2-T1|EPS（xk-1,xk）xk-1/22/)()(221)(2)]()([4)]}()([221)]()([221{12/1112/1112/12/112SHTxfhTxfhxfxfhxfxfhxfxfhTnnkknnknkkkknkkkkkn其中：hHxfSxfhbfafhTnkknkkn,)(,)()]()([2112/113-1-1梯形法——方法原理例：Debye-Einstein公式推导得到计算固体热容的公式为m0243md)1/(/9xxxVxexeXRC其中：TX/DmD为Debye温度，R为气体常数8.314JK-1mol-1已知固体的Debye温度如下：PbAgCuAlFeKClNaClCD/*K882153153984202272811910求在50，100，298.15，500，1500K时，各固体的热容。3-2-1-1Simpson法——问题的提出求积分badxxfS)(Simpson法是把积分区间分割成有限个小区间，在每个小区间上采用二次抛物线来近似被积函数f(x)的图形，近似求出小区间的面积，然后再将有限个小区间相加得到被积函数的近似值。xyy=f(x)xi-1xi+1xiy=g(x)hhSi定步长:nabh2/)()]2()(4)([3d)(22hxfhxfxfhxrqxpxSiiihxxiii3-2-1-2Simpson法——方法原理})]2(2))12((4[)]()({[312ninihafhiafbfafhS})]())2/1((2[2/)]()({[31nknkhafhkafbfafhS变步长：3/)(22nnnnTTTShHxfSxfhbfafhTnkknkkn,)(,)()]()([212/111其中：2/)(2SHTTnn3-2-1-2Simpson法——方法原理判据：EPSSSSSEPSSSSnnnnnnn22222/)(,1,1时当时当3-2-1-2Simpson法——方法原理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|DEPSNoRETURNYesNoN=N*2H=H/2T1=T2S1=S2Yes3-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显示程序显示输出3-1-3–1离散点数据的求积——方法原理实验时，得不到变量间的关系式，只测量到（xi,yi）的离散点数据。xyab方法：1.用插值程序求任意点的函数值。一元三点Lagrange插值：22])([)()(iijiikikjkyxxxxxpxf2．用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))/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|DEPSNoRETURNYesNoN=N*2H=H/2T1=T2S1=S2Yes(1)调用Lagrange一元三点插值F(A),F(B)(2)调用Lagrange一元三点插值F(A+(K-1/2)*H)3-1-3–2离散点数据的求积——程序框图纯气体的逸度由定义：pVfRTdlndd实m,pRTVVpRTVV实m,实m,实m,m,-理（1）代入（1）ppRTfRTd)(lnd积分，并取极低压力下气体视为理想气体，得pppRTpRTpf00dd1lnln逸度：pfφ为逸度系数（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=PB=0YN逸度系数FI=EXP(S),逸度FF=FI*B例2：已知固体Pb的热容Cp~温度T数据，求从15K到550K的固体Pb的焓变。TCHTTpd21例3：分子标准熵S及Cp~T数据，求500K时的熵S值。TTCSSTTpd)/(2112T1：298.15KT2:500K3-1-3–3离散点数据的求积——应用示例)exp(12SKKppSTRTHKKTTpp2112dln2已知数据TRTH~2TH~例4：合成氨反应322NHH23N21焓变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.67153-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两种物质相对含量之比。在色谱图上，色谱峰的面积与色谱分析中各物质的含量成正比。时间t/sec峰高/mm时间t/sec峰高/mm38013044032500.5531725866440266147569274218690282139768290510843298011919开始输入：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/S13-1-3–3离散点数据的求积——应用示例例6：在简单蒸馏釜内蒸馏1000Kg含C2H5OH质量分数为60%和H2O质量分数为40%的混合液。蒸馏结束时，残液中含C2H5OH质量分数为5%。试求残液的质量是多少千克？该体系的汽液平衡数据如下：其中x为液相中C2H5OH的质量分数，y为汽相中C2H5OH的质量分数。3-1-3–3离散点数据的求积——应用示例xy1xy1xx0.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简单蒸馏的雷利公式为：3-1-3–3离散点数据的求积——应用示例式中，F为原料液量，W为残液量；xF为原料液组成，xW为残液组成。FWxxxydxWFln开始输入：原料液量F,x,1/(y-x)输出：W结束两次调用离散点积分法子程序计算右侧积分S（其中调用Lagrange插值法子程序计算任意点的函数值）输入：积分上、下限XW，XF计算残液量W＝F/exp(S)显示程序显示输出3-1-3–3离散点数据的求积——应用示例3-2常微分方程的数值解法——引言一阶常微分方程的初值问题：00)(),(yxyyxfdxdyy一阶常微分方程组的初值问题：niyxyyyyxfyinii,,2,1)(),,,,(0021数值解法：寻求解y(x)在一系列离散点上的近似值,使y与x的关系近似满足y=F(x)步长：假定h为定值nxxx10nyyy,,,101iixxh（1）数值解的特点：找一个递推公式ihxxniyxyxiii0,0,0,,2,1)()(（1）式积分11d),(diiiixxxxxyxfy1d),()()(1iixxiixyxfxyxy（2）3-2常微分方程的数值解法——引言数值积分3-2-1-1Euler法及其改进——问题的提出例：异丙苯氧化反应的动力学ROOH)(RHHCHC]o[7356已知：t=0时的[RH]、[ROOH]，120℃时的反应速率常数k，计算0-14h每隔2h的[ROOH]。解：2/1]RH][ROOH[d]d[ROOHkt设x=[ROOH]，则[RH]=c-x)()(dd002/1txxxxcktx1d),()()(1iixxiixyxfxyxy右端利用矩形求积：得),(1iiiiyxhfyy——Euler公式显式几何意义：取切线的端点Pi+1作为yi+1xyxixi+1y=y(x)pipi+1h3-2-1-2Euler法及其改进——方法原理11(,)iiiiiiyyfxyxxh提高精度：（2）式右端利用梯形法求积：)],(),([2111iiiiiiyxfyxfhyy隐式改进：预报~校正预报：),(~1iiiiyxhfyy校正：)]~,(),([2111iiiiiiyxfyxfhyy编程：)(21),(),(111cpiiiiciiipyyyyxhfyyyxhfyy3-2-1-2Euler法及其改进——方法原理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法及其改进——程序框图开始输入：c,k,t0,X0,t(I)调用Euler法子程序计算X(I)输出：X(I),t(I)结束例1：异丙苯氧化反应的动力学3-2-1-4Euler法及其改进——应用示例显示程序显示输出例2：气相色谱仪及其实验过程的仿真根据物料平衡的原理，可以得出塔板j上气相物质的物料平衡方程式：jijcijjjxLdtGdtGdU,11(1≤j≤N，N为塔板总