SCPD.Lecture.2011.08.数值积分与常微分方程.pptVIP

数值积分与常微分方程的数值解法 Runge-Kutta法 常微分方程组四阶RK公式： 变步长：判据： （方程） （方程组） 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 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( ) 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 * 数值积分与常微分方程的数值解法 数值积分 常微分方程的数值解法 梯形法 Simpson法 离散点数据的求积 Euler法及其改进 Runge-Kutta法 1.f(x)函数形式已知，但其积分 不能表示成初等函数的闭合形式 2. f(x)函数形式未知，但其离散 数据表已给出 梯形法 方法原理： 基本思想：复化求积，即从近似计算为出发点，用有限项 的求和计算来代替从而求出定积分的近似值。 定步长： 求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 其中： Simpson法 问题的提出 例：Debye-Einstein公式推导得到计算固体热容的公式为 其中： ?D为Debye温度，R为气体常数8.314JK-1mol-1 已知固体的Debye温度如下： 1910 281 227 420 398 315 215 88 ?D/*K C NaCl KCl Fe Al Cu Ag Pb 求在50，100，298.15，500，1500K时，各固体的热容。 方法原理： Simpson法 求积分 Simpson法是把积分区间分割成有限个小区间，在每个小区间上采用二次抛物线来近似被积函数f(x)的图形，近似求出小区间的面积，然后再将有限个小区间相加得到被积函数的近似值。 x y y=f(x) xi-1 xi+1 xi y=g(x) h h Si 定步长: Simpson法 变步长： 其中： Simpson法 判据： 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 Simpson法 应用示例： 开始 输入： Debye温度T（5），精度EPS，温度THETA 输出：固体的热容Cv 结束 调用Simpson积分法子程序计算式右方积分值S2 计算：XM=THETA/T(I) (I=1,N) 输入：积分上下限A=10-4，B=XM B=0 Y N 固体的热容Cv=9R/XM**3*S2 离散点数据的求积 方法原理 实验时，得不到变量间的关系式，只测量到（xi,yi） 的离散点数据。 x y a b 方法 ： 1.用插值