常微分方程数值解上机

来源:互联网 发布:淘宝棉麻衬衣 编辑:程序博客网 时间:2024/04/27 17:20

   二步显式Adams法和Gear法求解 y0=1,步长分别为h=0.1h=0.05

1.       程序文本

      二步显式Adams

 

clc;

y(1)=1;

h=0.1;

y(2)=y(1)-2*h*y(1)+3*h;

n=1/h;

for i=2:n

 y(i+1)=y(i)-3*h*y(i)+h*y(i-1)+3*h;

end

t=0:h:1;

f=-0.5*exp(-2*t)+1.5;

r=plot(t,y,'*');

hold on,

g=plot(t,f);

set(g,'LineWidth',1,'color','blue');hold off;

w=abs(f(n+1)-y(n+1));

 

      二步显式Gear

clc;

y(1)=1;

h=0.05;

y(2)=y(1)+h*y(1);

n=1/h;

for i=2:n

y(i+1)=(4/3*y(i)-1/3*y(i-1)+2*h)/(1+4/3*h);

end

t=0:h:1;

f=-0.5*exp(-2*t)+1.5;

r=plot(t,y,'*');

hold on,

g=plot(t,f);

set(g,'LineWidth',1,'color','blue'),hold off;

w=abs(f(n+1)-y(n+1));

 

 应用显式中点法和44Gill方法求解初值问题 ,列出当步长h=0.2h=0.1时计算到终点t=1的数值解及其误差,并附算法的matlab程序。

1.       程序文本

    显式中点法

 clc;

y(1)=3;h=0.2;

t=0;

n=1/h;

for i=1:n

    t=t+h;

    Y=(1-0.5*h)*y(i)-0.5*h*t;

    y(i+1)=y(i)-h*t-h*Y-0.5*h^2;

end

t=0:h:1;

f=2*exp(-t)-t+1;

r=plot(t,y,'r*');

hold on,

g=plot(t,f);

set(g,'LineWidth',1,'color','blue'),hold off;

w=abs(f(n+1)-y(n+1))

y

 

 

44Gill方法

clc;

y(1)=3;h=0.05;

t=0;b=0.5^0.5

n=1/h;

for i=1:n

    t=t+h;

    Y1=y(i);

    Y2=y(i)-0.5*h*t-0.5*h*Y1;

    Y3=y(i)-0.5*h*t-(b-0.5)*h*Y1-(1-b)*h*Y2-0.5*(1-b)*h^2;

    Y4=y(i)-h*t+b*h*Y2-(1+b)*h*Y3-0.5*h^2;

    y(i+1)=y(i)-h*t-1/6*h*Y1-1/3*(1-b)*h*Y2-1/3*(1+b)*h*Y3-1/6*h*Y4-0.5*h^2;

end

t=0:h:1;

f=2*exp(-t)-t+1;

r=plot(t,y,'r*');

hold on,

g=plot(t,f);

set(g,'LineWidth',1,'color','blue'),hold off;

w=abs(f(n+1)-y(n+1))

y

 

原创粉丝点击