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

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

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

专业：**************班级：**********

班内序号：**

学号：********

姓名：******

任课教师：******

北京邮电大学

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

实验目标

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

y‘=y−2x

,

y0=1,

0≤x≤4

其中解析解为 y x=1+2x

实验原理

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

p n+1=u n−3+4

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

m n+1=p n+1+112

121

(c n−p n)

c n+1=1

9u n−u n−2+

3

h[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，u0

2）置h=(b-a)/N，t0=a，n=1

3）计算f n-1=f(t n-1,u n-1)

K1=hf n-1

K2=hf(t n-1+h/2, u n-1+K1/2)

K3=hf(t n-1+h/2, u n-1+K2/2)

K4=hf(t n-1+h, u n-1+K3)

u n= u n-1+1/6(K1 +2K2 +2K3 +K4)

t n=a+nh

4）输出(t n ,u n)

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

否则，置n+1→n，0→p0，0→c0，转6. 6）计算

f3=f(t3,u3)

t= t3+h

p=u0+4/3(2 f3–f2 +2f1)

m=p+112/121(c0-p0)

c=1/8(9u3- u1)+3/8h[f(t,m)+ 2 f3–f2]

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