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

3-2-2-1 Runge-Kutta法——问题的提出 例： 基元反应： k2 已知：[A]0=[B]0=1mol·L-1, [C]0=[D]0=[E]0=0 1mol·L-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法——应用示例 综合应用示例： ● 化学反应动力学的计算与计算机模拟 ● 二组元汽-液相平衡的计算机模拟 ● 在固体催化剂表面上对气体吸附等温 方程式的计算及过程的计算机模拟 * ● Euler法及其改进 ● Runge-Kutta法 ● 梯形法 ● Simpson法 ●离散点数据的求积 第三章 数值积分与常微分方程的数值解法 数值积分 常微分方程的数值解法 1.f(x)函数形式已知，但其积分不能表示成初等函数的闭合形式 2. f(x)函数形式未知，但其离散数据表已给出 3-1-1 梯形法——方法原理 基本思想： 复化求积，即从近似计算为出发点，用有限项 的求和计算来代替从而求出定积分的近似值。 定步长： 求f(x)在[a,b]上的定积分 x y 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温度如下： Pb Ag Cu Al Fe KCl NaCl C ?D/*K 88 215 315 398 420 227 281 1910 求在50，100，298.15，500，1500K时，各固体的热容。 3-2-1-1 Simpson法——问题的提出 求积分 Simpson法是把积分区间分割成有限个小区