解微分方程欧拉法-R-K法及其MATLAB实例

1、解微分方程的欧拉法，龙格-库塔法及其MATLAB简 单实例欧拉方法(E u ler method)用以对给定初值的常微分方程(即初值问题)求解山 分为前进EULER法、后退EULER法、改进的EULER法。金缺点3欧拉法简单地取 切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越 大。因此欧拉格式一般不用于实际计算。亠改进欧拉格式：金为提高精度，需要在 欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。 改进欧拉法的精度为二阶。A算法为：微分方程的本质特征是方程中含有导数项，数值解法的笫一步就是设法消除其导 数值。对于常微分方程：xG a, by (a

2、) = yO可以将区间a, b分成n段，那么方程在第xi点有* (xi) = f ( x i ,y (xi),再用向前差商近似代替导数则为：(畑+1畑)=畑畑)在这里，h是步长，即相邻两个结点间的距离。因此可以根据xi点和yi点的 数值计算出yi+1来：如=y-i +hx yi)i 二 0,1,2, L这就是向前欧拉格式。改进的欧拉公式：将向前欧拉公式中的导数f( xi,yi)改为微元两端导数的平均，即上式便是梯形的欧拉公式。可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式:数值分析中，龙格-库塔法（Ru nge-Kutta）是用于模拟常微分方程的解的重要的 一类隐式或显式

3、迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为f （x n, yn）,而改进的欧拉公式将导数项取为两端导数的平均。龙格-库塔方法的基本思想：在区间xn, xn+1内多取儿个点，将他们的斜率加权平均，作为导数的近似。龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“ RK4”或者就是 “龙格库塔法”。令初值问题表述如下。则，对于该问题的RK4 111如下方程给出：如+1 = yn + 6（应1 + 22 + 2 应3 + 魅）其中ki = f (tn, yn)fcz = f (tn + 另如 + fcl = f (如+,?k +用2也=f (tn + h,y

4、n + hk3)这样，下一个值（几G由现在的值（）加上时间间隔（力）和一个估算的斜率的乘 积决定。该斜率是以下斜率的加权平均：k 1是时间段开始时的斜率严k2是时间段中点的斜率，通过欧拉法采用斜率kl来决定y在点tn + h/2的值Mk3也是中点的斜率，但是这次采用斜率k 2决定y值沖k 4是时间段终点的斜率，其y值用k3决定。当四个斜率取平均时，中点的斜率有更大的权值：ki + 22 + 2 应3 + 际slope =.bRK4法是四阶方法，也就是说每步的误差是片阶，而总积累误差为月阶。注意上述公式对于标量或者向量函数（y可以是向量）都适用。例子:下面给出了数值求解该微分方程的简单程序。其中

5、yl,y2, y3, y 4分别为向前欧拉公式，改进的欧拉公式，4级4阶龙格一库 塔公式及精确解。h二0. 1严 x=0:h:l ; Ayl=zeros(size(x);yl(l) = l;y2=zeros(size ( x) py2(lj = 1 ; y3 = ze r o s (s i ze(x); W 二 1;a for i l=2:lengt h (x )金 y l(il)=yl(il<)+ h * (yl(i 1-1) - 2 * x(i i-D /yi(ii-i);kl = y2(il-l) -2*x(il-l) /y2(il -1);k2=y2(il-l) + h*kl-2

6、*x( i l)/(y2(il-l)+h * kl)? y2(il) = y2(i l-l)+h*(kl+k2)/2?kl=y2(il -l)-2*x(il-l)/y2(i 1-1);k2 = y2 (il-1 ) +h*k 1/2-2* ( x(i 1 - l) + h/2 ) / ( y2(il- l)+h*kl/2)?k3= y2(il-l) + h*k2/2-2 * ( x(i l-l) + h/2 ) / (y2(i 1-1) +h*k2/2);k4= y 2(il-l)+h* k3-2*(x(il-1 ) + h) / (y 2 (il-l)+h*k3); y3(il)-y3(il - 1)+ ( kl + 2*k2+2*k3 + k4)*h/6;e n cMy4二sqrt ( 1 +2*x);%plo t (x , yl,x f y2,x,y3,