常微分方程的数值解法2010.ppt

对一阶常微分方程的初值问题，其一般形式是

(1)
在下面的讨论中，假定f(x,y)连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使得

则初值问题(1)的解必定存在且唯一。常微分方程的数值解法
所谓数值解法，就是要求问题(1)的解在若干点：

处的近似值yi(i=0,1,2…n)的方法，yi称为问题(1)的数值解。相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定hn为定长，称为定步长，这时节点可表示为

数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。一、欧拉（Euler）法
1.1Euler公式
欧拉（Euler）方法是解初值问题的最简单的数值方法。采用差分方法解初值问题

在的近似解
则Euler法的计算格式DIMENSIONX(0:10),Y(0:10)
DOUBLEPRECISIONX,Y,H,X0,Y0,N1
F(X,Y)=X-Y
open(2,file="aout.txt")
H=1.0/10
N=10
X(0)=0.0
Y(0)=1.0
DOI=0,N-1
Y(I+1)=Y(I)+F(X(I),Y(I))*H
X(I+1)=X(I)+H
ENDDO
10 FORMAT(1X,2D15.6/)
WRITE(*,10)(X(I),Y(I),I=1,10)
WRITE(2,10)(X(I),Y(I),I=1,10)
END1.2梯形公式
为了提高精度,对方程的两端在区间上
积分得，

改用梯形方法计算其积分项，即1.3改进的欧拉公式
显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。
先用欧拉公式(2)求出一个初步的近似值
,称为预测值,它的精度不高,再用梯形公式(5)对它校正一次,即迭代一次,求得yi+1,称为校正值,这种预测—校正方法称为改进的欧拉公式：DIMENSIONX(0:10),Y(0:10)
DOUBLEPRECISIONX,Y,H,X0,Y0,N1
F(X,Y)=Y-2*X/Y
H=1.0/10
N=10
X(0)=0.0
Y(0)=1.0
DOI=0,N-1
X(I+1)=(I+1)*H
Y(I+1)=Y(I)+H*F(X(I),Y(I))
Y(I+1)=Y(I)+(F(X(I),Y(I))+F(X(I+1),Y(I+1)))*H/2
ENDDO
10 FORMAT(1X,2D15.6/)
WRITE(*,10)(X(I),Y(I),I=1,10)
END设为节点上的近似解，
则有改进的Euler格式为二、龙格-库塔（Runge-Kutta）法
2.1龙格-库塔(Runge-Kutta)法的基本思想
Euler公式可改写成2.1高阶龙格-库塔(Runge-Kutta)法的构造yn+1在(xn，yn)处的Taylor展开式为：要使近似公式(8)的局部截断误差为O(h3)，则应要求(10)和(11)式前三项相同：DIMENSIONX(0:10),Y(0:10)
DOUBLEPRECISIONX,Y,H,X0,Y0,N1,k1,k2
F(X,Y)=Y-2*X/Y
H=1.0/10
N=10
X(0)=0.0
Y(0)=1.0
DOI=0,N-1
X(I+1)=(I+1)*H
K1=F(X(I),Y(I))
K2=F(X(I+1),Y(I)+H*K1)
Y(I+1)=Y(I)+H*(K1+K2)/2.0
ENDDO
10 FORMAT(1X,2D15.6/)
WRITE(*,10)(X(I),Y(I),I=1,10)
END2.2四阶龙格—库塔法DIMENSIONX(0:10),Y(0:10)
DOUBLEPRECISIONX,Y,H,X0,Y0,N1,K1,K2,K3,K4
F(X,Y)=X-Y
H=0.2
N=5
X(0)=0.0
Y(0)=1.0
DOI=0,N-1
X(I+1)=(I+1)*H
K1=F(X(I),Y(I))
K2=F(X(I)+0.5*H,Y(I)+0.5*H*K1)
K3=F(X(I)+0.5*H,Y(I)+0.5*H*K2)
K4=F(X(I)+H,Y(I)+H*K3)
Y(I+1)=Y(I)+H*(K1+2*K2+2*K3+K4)/6
ENDDO
10 FORMAT(1X,2D15.6/)
WRITE(*,10)(X(I),Y(I),I=1,10)
END三、一阶常微分方程组
基于一阶常微分方程初值问题的各种数值解法可类似得到一阶常微分方程组的各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。
3.1一阶常微分方程组的初值问题SUBROUTINEGELR1(T,Y,M,H,N,Z,F,D)
DIMENSIONY(