数值积分算法误差分析

实验名称：数值积分算法误差分析

1.实验原理

1)欧拉法原理

在数学和计算机科学中，欧拉方法（Euler method）命名自它的发明者莱昂哈德·欧拉，是一种一阶数值方法，用以对给定初值的常微分方程(即初值问题)求解。它是一种解决常微分方程数值积分的最基本的一类显型方法(Explicit method)。

微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较低，一般不在工程中单独进行运算。所谓数值求解，就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。对于常微分方程：

0 )

(

]

,

[ ),

,

(

y a

y

b

a

x

y

x

f dx

dy

=∈

=

可以将区间[a,b]分成n段，那么方程在第x i点有y'(x i) =

f(x i,y(x i))，再用向前差商近似代替导数则为

：

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

1

，i=0,1,2,L

这就是欧拉格式，若初值y i+ 1是已知的，则可依据上式逐步算出数值解y1,y2,L。

1)龙哥库塔法原理

数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。

龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。令初值问题表述如下。

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

其中

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

∙

k 1是时间段开始时的斜率； ∙ k 2是时间段中点的斜率，通过欧拉法采用斜率k 1来决定y 在点

t n + h /2的值；

∙

k 3也是中点的斜率，但是这次采用斜率k 2决定y 值； ∙ k 4是时间段终点的斜率，其y 值用k 3决定。

当四个斜率取平均时，中点的斜率有更大的权值：

RK4法是四阶方法，也就是说每步的误差是h 5阶，而总积累误差为h 4阶。

2 实验方案

1)目标方程为

对该方程进行分析可知，该方程一阶微分方程，为y 对t 的微分。

2）子函数程序：

subf.m 文件

function dy = subf(t,y)

dy = (1-t)*y;

3）目标方程解析解

>> H = dsolve('Dy=(1-t)*y','y(0)=1','t') ⎪⎩⎪⎨⎧=-=1

)0()1(y y t dt dy

得到：H = exp(1)^(1/2)/exp((t - 1)^2/2)

3.实验程序

1)欧拉法仿真脚本程序

clear %欧拉法

h = 0.1;

tN = 10;

t = 0:h:tN;

N = length(t);

y(1) = 1;

j = 1;

for j = 1:N

dataY(j) = y(1);

yn = y+h*subf(t(j),y);

y = yn;

end

%解析解法

ezplot(' exp(1)^(1/2)/exp((t - 1)^2/2)',[0 10]) hold on

plot(t,dataY,'r+');

legend('解析解法','欧拉法')

2)龙哥库塔法仿真脚本程序

clear %龙哥库塔法

h = 0.1;

tN = 10;

t = 0:h:tN;

N = length(t);

y(1) = 1;

j = 1;

for j = 1:N

dataY(j) = y(1);

K1 = subf(t(j),y);

K2 = subf(t(j)+h/2,y+1/2*h*K1);

K3 = subf(t(j)+h/2,y+1/2*h*K2);

K4 = subf(t(j)+h,y+h*K3);

yn = y+h/6*(K1 + 2*K2 + 2*K3 +K4);

y = yn;

end

%解析解法

ezplot(' exp(1)^(1/2)/exp((t - 1)^2/2)',[0 10]) hold on

plot(t,dataY,'ro');

legend('解析解法','四阶龙哥库塔法')

5.仿真结果

（1）欧拉法步长h = 0.1

（2）欧拉法步长h = 0.01

（3）四阶龙哥库塔法步长h = 0.1

（4）四阶龙哥库塔法步长h = 0.01

（5）四阶龙哥库塔法步长h = 0.5

5.实验结果分析

由仿真结果可知,欧拉法取不同步长，其误差有较大差异，步长越小，误差越少，但计算量也将越大。因此使用时应选取合适的步长。四阶龙哥库塔法取不同步长，仿真曲线不变误差不变，由此可知，四阶龙哥库塔法精度与步长无关，但由步长为h = 0.5 在后面的数据有很大的偏离，因此步长不宜去太大。欧拉法与四阶龙格库塔法相比较，前者精度受步长影响较大，后者较小。