solutionofillsystem.m

病态系统求解方法比较：龙格库塔4固定步长、变步长、吉尔法

function [d,y21,y31,y41]=allnew(intialH,N)
%-----------------------------------------------------------------------
%  输入变量为初始步长和采样点数
%  输出为固定步长法、变步长法(ode45）、Gear法的误差
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
%                        求解:
%                        DY=AY,
%                  A=[-21  19  -20;
%                      19 -21   20;
%                      40 -40  -40];
%                   Y0=[1 0 -1]
%-----------------------------------------------------------------------
%        求A的特征值
%-----------------------------------------------------------------------
A=[-21  19  -20;    
    19 -21   20;
    40 -40  -40];
y0=[1 0 -1];
d=eig(A);
%-----------------------------------------------------------------------
%      RK4固定步长方法(编制RK4公式)
%      当步长大于0.05时,开始发散
%-----------------------------------------------------------------------
tic
h=intialH;  
y2(1,:)=y0;
t(1)=0;
for i=1:N                 

    K1=f(t(i),y2(i,:));

    K2=f(t(i)+h/2,y2(i,:)+h*K1/2);

    K3=f(t(i)+h/2,y2(i,:)+h*K2/2);

    K4=f(t(i)+h,y2(i,:)+h*K3);

    y2(i+1,:)=y2(i,:)+h*(K1+2*K2+2*K3+K4)/6;

    t(i+1)=t(i)+h;
end
toc
%-----------------------------------------------------------------------
%      RK4变步长方法(调用ode45函数)
%-----------------------------------------------------------------------
tic
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,y3] = ode45(@rigid,[0:intialH:intialH*N],[1 0 -1],options);  
toc
%-----------------------------------------------------------------------
%      RKM变步长法
%-----------------------------------------------------------------------
%tic
%h=0.01;  %当步长大于0.05时,开始发散
%N=500;
%y3(1,:)=y0;
%t(1)=0;
%for i=1:N                 

%    K1=f(t(i),y3(i,:));
 
%    K2=f(t(i)+h/3,y3(i,:)+h*K1/3);

 %   K3=f(t(i)+h/2,y3(i,:)+h*(K1+K2)/6);
 
 %   K4=f(t(i)+h/2,y3(i,:)+h*(K1+3*K3)/8);
 
 %   K5=f(t(i)+h,y3(i,:)+h*(K1-3*K3+4*K4)/2);

%    y3(i+1,:)=y3(i,:)+h*(K1+4*K4+K5)/6;

 %   t(i+1)=t(i)+h;
%end
%toc
%-----------------------------------------------------------------------
%     病态系统求解方法--ode15s法(改进的吉尔法)
%-----------------------------------------------------------------------
tic
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,y4] = ode15s(@rigid,[0:intialH:intialH*N],[1 0 -1],options);  
toc
%-----------------------------------------------------------------------
%      解析解法
%-----------------------------------------------------------------------
y11=1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)+1/2*exp(-40*t).*sin(40*t) ;
y12=-1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)-1/2*exp(-40*t).*sin(40*t);
y13=exp(-40*t).*sin(40*t)-exp(-40*t).*cos(40*t);
%-----------------------------------------------------------------------
%     显示图像
%-----------------------------------------------------------------------
subplot(3,1,1)
plot(t,y11,'-g',t,y2(:,1),'-.r',t,y3(:,1),'--b',t,y4(:,1),'.m');
xlabel('t/sec'),ylabel('y1')         

subplot(3,1,2)
plot(t,y13,'-g',t,y2(:,2),'-.r',t,y3(:,2),'--b',t,y4(:,2),'.m');
xlabel('t/sec'),ylabel('y2')

subplot(3,1,3)
plot(t,y13,'-g',t,y2(:,3),'-.r',t,y3(:,3),'--b',t,y4(:,3),'.m');
xlabel('t/sec'),ylabel('y3')

%-----------------------------------------------------------------------
%        微分方程
%-----------------------------------------------------------------------
function dy =f(t,y)    
dy = zeros(1,3);    
dy=(A*y')';
end
function dy = rigid(t,y)
dy = zeros(3,1);   
dy=A*y;
end
%-----------------------------------------------------------------------
%     计算精度
%-----------------------------------------------------------------------
y21=zeros(1,4);
y31=zeros(1,4);
y41=zeros(1,4);
for i=1:1:N+1
    y21(1)=(y2(i,1)-y11(i))^2+y21(1);
    y21(2)=(y2(i,2)-y12(i))^2+y21(2);
    y21(3)=(y2(i,3)-y13(i))^2+y21(3);
end
y21(4)=y21(1)+y21(2)+y21(3);
for i=1:1:N+1
    y31(1)=(y3(i,1)-y11(i))^2+y31(1);
    y31(2)=(y3(i,2)-y12(i))^2+y31(2);
    y31(3)=(y3(i,3)-y13(i))^2+y31(3);
end
y31(4)=y31(1)+y31(2)+y31(3);
for i=1:1:N+1
    y41(1)=(y4(i,1)-y11(i))^2+y41(1);
    y41(2)=(y4(i,2)-y12(i))^2+y41(2);
    y41(3)=(y4(i,3)-y13(i))^2+y41(3);
end
y41(4)=y41(1)+y41(2)+y41(3);
y21,y31,y41
%-----------------------------------------------------------------------
end