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

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法是把积分区间分割成有限个小区