pengtaiyu2.m

四阶龙格库塔法,大学MALAB实践课程 需要的请下载

clear
close all
clc
select=input('1.4阶段龙格库塔法 2.阿达姆斯法 3.不同步长四阶龙格库塔法 4.病态方程h=0.01 5.病态方程h=0.04:');
if select==1
a=[-0.23 -0.56;0.56 0];
b=[1;0];
c=[1.12 5.67];
d=[0];u=10;
 x=[1;0];y=c*x+d*u;t=0;Tf=80;h=0.1;n=Tf/h;ti=x;
 for i=1:n
    k1=a*ti+b*u;
    k2=a*(ti+h/2*k1)+b*u;
    k3=a*(ti+h/2*k2)+b*u;
    k4=a*(ti+h*k3)+b*u;
    x1=ti+h*(k1+2*k2+2*k3+k4)/6;
    x=[x,x1];
    ti=x1;
    y=[y,c*ti];
    t=[t,t(end)+h];
end  
disp('四阶龙格-库塔法的前10步的迭代结果:')
    for j=1:10
     fprintf('y(%d)=%d\n',j,y(j))
 end
plot(t,y,'r')
figure(2);
plot(t,x(1,:),t,x(2,:));
grid
hold on
end
if select==2
X=[0;0];y=0;R=1;t=0;
Tf=80;
h=0.10;
n=Tf/h;
A=[-0.5 -0.6;0.6 0];
B=[1;0];
C=[1 5.6];
D=[0];
E=[1 0];
F=[0 1];
X10=0;X20=0;
   for i=1:n
    K1=A*X+B*R;
    K2=A*(X+h/2*K1)+B*R;
    K3=A*(X+h/2*K2)+B*R;
    K4=A*(X+h*K3)+B*R;
    X=X+h*(K1+2*K2+2*K3+K4)/6;
    X10=[X10,E*X];
    X20=[X20,F*X];
    y=[y,C*X];
    t=[t,t(i)+h];
end
for i=1:4
Adams(i)=y(i);
end
for j=4:n
  Adams(j+1)=y(j)+(h/24)*(55*y(j)-59*y(j-1)+37*y(j-2)-9*y(j-3));  
end
plot(t, Adams,'r')
grid on
end
if select==3
   a=[-0.5 -0.6;0.6 0];
b=[1;0];
c=[1 5.6];
d=[0];u=1;
 x=[1;0];y=c*x+d*u;t=0;Tf=80;h=0.8;n=Tf/h;ti=x;
for i=1:n
    k1=a*ti+b*u;
    k2=a*(ti+h/2*k1)+b*u;
    k3=a*(ti+h/2*k2)+b*u;
    k4=a*(ti+h*k3)+b*u;
    x1=ti+h*(k1+2*k2+2*k3+k4)/6;
    x=[x,x1];
    ti=x1;
    y=[y,c*ti];
    t=[t,t(end)+h];
end  
disp('四阶龙格-库塔法的前10步的迭代结果:')
    for j=1:10
     fprintf('y(%d)=%d\n',j,y(j))
 end
plot(t,y,'r')
figure(2);
plot(t,x(1,:),t,x(2,:));
grid
hold on
end
if select==4
a=[-21 19 -20;19 -21 20;40 -40 -40];
b=[1;0;0];
c=[-21 30 10];
d=[0];u=1;
 x=[1;0;-1];y=c*x+d*u;t=0;Tf=1;h=0.01;n=Tf/h;ti=x;
for i=1:n
    k1=a*ti+b*u;
    k2=a*(ti+h/2*k1)+b*u;
    k3=a*(ti+h/2*k2)+b*u;
    k4=a*(ti+h*k3)+b*u;
    x1=ti+h*(k1+2*k2+2*k3+k4)/6;
    x=[x,x1];
    ti=x1;
    y=[y,c*ti];
    t=[t,t(end)+h];
end  
disp('四阶龙格-库塔法的前10步的迭代结果:')
    for j=1:10
     fprintf('y(%d)=%d\n',j,y(j))
 end
plot(t,y,'r')
figure(2);
plot(t,x(1,:),t,x(2,:));
grid
hold on
end
if select==5
a=[-21 19 -20;19 -21 20;40 -40 -40];
b=[1;0;0];
c=[-21 30 10];
d=[0];u=1;
 x=[1;0;-1];y=c*x+d*u;t=0;Tf=1;h=0.04;n=Tf/h;ti=x;
for i=1:n
    k1=a*ti+b*u;
    k2=a*(ti+h/2*k1)+b*u;
    k3=a*(ti+h/2*k2)+b*u;
    k4=a*(ti+h*k3)+b*u;
    x1=ti+h*(k1+2*k2+2*k3+k4)/6;
    x=[x,x1];
    ti=x1;
    y=[y,c*ti];
    t=[t,t(end)+h];
end  
disp('四阶龙格-库塔法的前10步的迭代结果:')
    for j=1:10
     fprintf('y(%d)=%d\n',j,y(j))
 end
plot(t,y,'r')
figure(2);
plot(t,x(1,:),t,x(2,:));
grid
hold on
end