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

解微分方程的欧拉法，龙格-库塔法及其ＭATLＡB简单实例

欧拉方法（Ｅｕler method)用以对给定初值的常微分方程(即初值问题)求解ﻫ

分为前进EULER法、后退EULER法、改进的EULEＲ法。ﻫ缺点:ﻫ欧拉法简单地取

切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。ﻫ改进欧拉格式：ﻫ为提高精度,需要在

欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。ﻫ算法为：

微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值。

对于常微分方程:

x∈[a,b］

ｙ(a) = y0

可以将区间[ａ,ｂ]分成n段，那么方程在第xi点有y＇（xi) = f（ｘｉ,y （xi))，再用向前差商近似代替导数则为：

在这里，ｈ是步长，即相邻两个结点间的距离。因此可以根据ｘi点和yｉ点的数值计算出ｙi+1来：

i=０,1,2,L

这就是向前欧拉格式。

改进的欧拉公式：

将向前欧拉公式中的导数f（ｘi,yi）改为微元两端导数的平均，即

上式便是梯形的欧拉公式。

可见,上式是隐式格式,需要迭代求解。为了便于求解，使用改进的欧拉公式：

数值分析中，龙格-库塔法（Rｕnge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为ｆ（ｘn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想:

在区间[xn，ｘn+1]内多取几个点,将他们的斜率加权平均,作为导数的近似。

龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“ＲK4”或者就是“龙格库塔法”。ﻫ

令初值问题表述如下。

则，对于该问题的RK4由如下方程给出：

其中

这样，下一个值(yｎ+1)由现在的值（y n)加上时间间隔（h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均:

k１是时间段开始时的斜率;ﻫ k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn ＋ h/2的值；ﻫ k3也是中点的斜率,但是

这次采用斜率k２决定y值；ﻫ k４是时间段终点的斜率,其y值用k3决定。当四个斜率取平均时,中点的斜率有更大的权值：

RK4法是四阶方法，也就是说每步的误差是h5阶,而总积累误差为h4阶。

注意上述公式对于标量或者向量函数（y可以是向量）都适用。

例子：

下面给出了数值求解该微分方程的简单程序。

其中y1,y2，y３,ｙ4分别为向前欧拉公式,改进的欧拉公式，4级4阶龙格－库塔公式及精确解。

h＝0．１;ﻫｘ=0:ｈ:１；ﻫy1=zｅrｏs(size(x))；

y1(1)=1;ﻫy2=zeros(size（x)）;ﻫy2(1)＝１；ﻫy3＝zｅｒoｓ(sｉze(ｘ))；

ﻫy3(1)=1;

ﻫｆor ｉ1=2:lengtｈ(x）ﻫｙ1(i1)=y1(i1-1)+ｈ＊(y1(i１-1）－２＊ｘ(ｉ1－1）／y1(ｉ1-１）);

ｋ1＝y２(i1-1）-2*x(i1-1）／y2(i1－1）;

k2=y2(i1-1）+ｈ*k1－2＊x（ｉ１)/(y2(i1-1)+h＊k1）;ﻫy2(i１)＝y2(ｉ1－1)+ｈ*(k1+k2)/２;ﻫ

ｋ1=y2(i１－1)-２*x(i1-1)/y2(ｉ１-1）；

k2＝y2（i1-１）+h＊k１/2-2＊（x(ｉ1－1)＋h/2）/（ｙ2(ｉ1－

１)+h*k1/2);ﻫk3＝ｙ2(i1-１)＋h*k2/2-２＊（x(ｉ1-１)＋h/2）／(y2(i１-1）+h*ｋ２/2）;

k4=ｙ2(i1-1)+h*ｋ3-2*(x(i1-１）+ｈ)／(ｙ2（i1-１)+h*k3);

y3(i1)=y3(i１－１)+（k1＋2*k2+2*k3＋k４)*h/6;

eｎdﻫy4=sqrt（１+2*x);

%plｏｔ(ｘ，y1,x，y2,x,y３,x,y4)ﻫ%ｌｅgeｎd(＇y1','y２','y３','ｙ４')

plot(ｘ，ｙ4-y1,ｘ,y4-y2,x,y4-y3)

ｌegeｎd('y１','ｙ2','ｙ3')