龙格库塔方法的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Ô¤²â-УÕýËã·¨')

grid on

legend('Miline-HammingÔ¤²â-УÕýËã·¨','½âÎö½â')

%解真值

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所示