龙格库塔方法的Miline-Hamming预测-校正算法实验报告知识讲解

龙格库塔方法的

M i l i n e-H a m m i n g预测-校正算法实验报告

2011-2012学年第2学期实验报告

实验名称：微分方程数值解实验学院：******

专业：**************

班级：**********

班内序号：**

学号：********

姓名：******

任课教师：******

北京邮电大学

时间：****年**月**日

实验目标

用多环节Miline-Hamming 预测-校正算法求下列方程的解

{y‘=y ?2x

y ,y (0)=1, 0≤x ≤4 其中解析解为 y (x )=√1+2x

实验原理

计算龙格库塔显示公式计算预测值，然后用隐式公式进行校正。Miline-Hamming 预测-校正公式为

{

p n+1=u n?3+4

3

h(2f n ?f n?1+f n?2)

m n+1=p n+1+112

121(c n

?p n )

c n+1=18

(9u n ?u n?2)+38h[f (t n+1,m n+1)+2f n ?f n?1]

u n+1=c n+1?9(c n+1

?p n+1)

其对应的算法流程为

1） 输入a ，b ，f(t,u)，N ，u 0 2） 置h=(b-a)/N ，t 0=a ，n=1 3） 计算f n-1=f(t n-1,u n-1)

K 1=hf n-1

K 2=hf(t n-1+h/2, u n-1+K 1/2) K 3=hf(t n-1+h/2, u n-1+K 2/2) K 4=hf(t n-1+h, u n-1+K 3)

u n = u n-1+1/6(K 1 +2K 2 +2K 3 +K 4) t n =a+nh

4） 输出(t n ,u n )

5） 若n<3，置n+1→n ，返回3；

否则，置n+1→n ，0→p 0，0→c 0，转6.

6） 计算

f 3=f(t 3,u 3) t= t 3+h

p=u 0+4/3(2 f 3 –f 2 +2f 1) m=p+112/121(c 0-p 0)

c=1/8(9u 3- u 1)+3/8h[f(t,m)+ 2 f 3 –f 2] u=c-9/121(c-p)

输出(t,u)

7）如n

否则停止。

实验过程

我们不妨设步长h=0.2，编程实现如下：

clear

clf

clc

%直接求解微分方程

y=dsolve('Dy=y-2*t/y','y(0)=1','t');

%Miline-Hamming预测-校正法

h=0.2;

t=0:h:4;

n=length(t);

u=zeros(1,n);

u(1)=1;

zbu(1,1)=t(1);

zbu(2,1)=u(1);

f=zeros(1,n);

p=zeros(1,n);

c=zeros(1,n);

m=zeros(1,n);

for i=2:n

if i-1<=3

f(i-1)=u(i-1)-2*t(i-1)/u(i-1);

k1=h*f(i-1);

k2=h*((u(i-1)+k1/2)-2*(t(i-1)+h/2)/(u(i-1)+k1/2)); k3=h*((u(i-1)+k2/2)-2*(t(i-1)+h/2)/(u(i-1)+k2/2)); k4=h*((u(i-1)+k3)-2*(t(i-1)+h)/(u(i-1)+k3));

u(i)=u(i-1)+1/6*(k1+2*k2+2*k3+k4);

zbu(1,i)=t(i);

zbu(2,i)=u(i);

else

f(i-1)=u(i-1)-2*t(i-1)/u(i-1);

p(i)=u(i-4)+4/3*h*(2*f(i-1)-f(i-2)+2*f(i-3));

m(i)=p(i)+112/121*(c(i-1)-p(i-1));

c(i)=1/8*(9*u(i-1)-u(i-3))+3/8*h*(m(i)-2*t(i)/m(i)+2*f(i-1)-f(i-2));

u(i)=c(i)-9/121*(c(i)-p(i));

zbu(1,i)=t(i);

zbu(2,i)=u(i);

end

end

zbu

%作图

plot(t,u,'r*','markersize',10)

hold on,

ezplot(y,[0,4])

hold on,

title('Miline-Hamming?¤2a-D￡?y??·¨')

grid on

legend('Miline-Hamming?¤2a-D￡?y??·¨','?a???a')

%解真值

h=0.2;

t=0:h:4;

n=length(t);

for i=1:n

y(i)=(1+2*t(i))^(1/2);

zby(1,i)=t(i);

zby(2,i)=y(i);

end

zby

我们可以得到计算后的结果图像如图1所示

图1 Miline-Hamming预测-校正法与解析解比较（h=0.2）同时我们可以得到Miline-Hamming预测-校正法和解析解在各点处的数值分别如下表1所示：

表1 Miline-Hamming预测-校正法与解析解在各点数值比较（h=0.1）

为了评判Miline-Hamming预测-校正法的算法精度，在这里我们利用相对误差的概念进行评判。对于Miline-Hamming预测-校正法的每个的估计值有：

相对误差=|估计值－真值|

真值

从而我们可以通过计算得到如下的相对误差表：

t坐00.20000.40000.60000.8000 1.0000 1.2000 1.4000 1.6000

t坐 1.8000 2.0000 2.2000 2.4000 2.6000 2.8000 3.0000 3.2000 3.4000

很明显，当我们对各点处的相对误差取平均后，该平均值小于0.01。因

此，我们可以认为Miline-Hamming预测-校正法的在h=0.2时的算法精度相对较高，所得到的结果与真值较为接近。

接下来我们在对h=0.5和h=0.1的情况进行计算，可以得到结果如下表3和表4所示

00.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 t坐

标

H=.5 1.0000 1.4155 1.7355 2.0084 2.2517 2.4886 2.7454 3.0750 3.5964 真 1.0000 1.4142 1.7321 2.0000 2.2361 2.4495 2.6458 2.8284 3.0000

t坐00.10000.20000.30000.40000.50000.60000.70000.8000

值

表4 Miline-Hamming预测-校正法在各点值及相对误差比较（h=0.1）

其中表3中相对误差的平均值为0.0393。而表4中的误差值小于是10的-3次方，在此不列出。

可以的到结果图像如下图2和图3所示：

图1 Miline-Hamming预测-校正法与解析解比较（h=0.5）

图1 Miline-Hamming预测-校正法与解析解比较（h=0.1）

实验结果

通过以上计算，我们可以得到如下的结论：

1．Miline-Hamming预测-校正法的计算精度相对较高，，当步长h=0.2时平均相对误差已小于0.01，因此可以认为这种方法可以得到和解析解较为接近的数值解。

2．伴随着步长的增加，我们可以发现相对误差的平均值随之减小。因此我们认为当步长h越小时计算精度越高。因此在计算能力允许的范围内，选取步长越小可以得到更加精确的结果。

3．在利用Miline-Hamming预测-校正法的过程中，前4次迭代的结果会对第五轮求得的数值产生影响。因此，一旦前四轮轮迭代所得的结果有偏差，下一轮结果的偏差将大于之前的偏差。因此会导致伴随迭代次数的增加而产生更大的偏差。