eighta.m

判断微分方程是否为病态方程

function []=rk4()
%-----------------------------------------------------------------------
%求解:
%      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];
K=1;
d=eig(A)
%-----------------------------------------------------------------------
%      解析解法
%-----------------------------------------------------------------------
tic
t0=0:0.01:5;
s=sym('s');
b=[s,0,0;
    0,s,0;
    0,0,s];
c=b-A;
d=inv(c);
ee=ilaplace(d);
%y1=[0;0;0];
y1= ee * [1;0;-1];   %解析解
subplot(4,3,1)
ezplot(y1(1),t0); hold on
subplot(4,3,2)
ezplot(y1(2),t0); hold on
subplot(4,3,3)
ezplot(y1(3),t0);hold on
toc
%-----------------------------------------------------------------------
%      RK4固定步长方法(编制RK4公式)
%-----------------------------------------------------------------------
tic
h=0.02;  %当步长大于0.05时,开始发散
N=5/h;
y2(:,1)=y0;
t(1)=0;
for i=1:N-1                 

    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
subplot(4,3,4)
plot(t,y2(1,:),'-r');
hold on
subplot(4,3,5)
plot(t,y2(2,:),'-.r');
hold on
subplot(4,3,6)
plot(t,y2(3,:),'--r');
hold on
toc
%-----------------------------------------------------------------------
%      RK4变步长方法(调用ode45函数)
%-----------------------------------------------------------------------
tic
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,y3] = ode45(@f,[0:0.01:2],[1 0 -1],options);  
subplot(4,3,7)
plot(t,y3(:,1),'-b');
hold on
subplot(4,3,8)
plot(t,y3(:,2),'-.b');
hold on
subplot(4,3,9)
plot(t,y3(:,3),'--b');
hold on
toc
%-----------------------------------------------------------------------
%     病态系统求解方法--ode15s法
%-----------------------------------------------------------------------
tic
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[t,y4] = ode15s(@f,[0:0.01:2],[1 0 -1],options);  
subplot(4,3,10)
plot(t,y4(:,1),'-m');
hold on
subplot(4,3,11)
plot(t,y4(:,2),'-.m');
hold on
subplot(4,3,12)
plot(t,y4(:,3),'--m');
hold on
toc
%-----------------------------------------------------------------------
%        微分方程
%-----------------------------------------------------------------------
function dy =f(t,y)    
dy = zeros(3,1);    
dy(1) =-21*y(1)+19*y(2)-20*y(3);
dy(2) =19*y(1)-21*y(2)+20*y(3);
dy(3) =40*y(1)-40*y(2)-40*y(3);
end
%-----------------------------------------------------------------------
%     计算精度
%-----------------------------------------------------------------------
for i=1:1:200
    t=0.01*i;
    y1(i)=1/2*exp(-40*t).*cos(40*t)+1/2*exp(-2*t)+1/2*exp(-40*t).*sin(40*t);
end
%-----------------------------------------------------------------------
%     加图像注释
%-----------------------------------------------------------------------
%title('病态系统求解实例');
%xlabel('t/sec');
%ylabel('y');
%text(3,0.9,'绿线--解析解','color','g');
%text(3,0.8,'蓝线--RK4变步长的ode45法','color','b');
%text(3,0.7,'红线--RK4固定步长法','color','r');
%text(3,0.6,'紫线--适合刚性变阶的ode15s法','color','m');
%text(2.75,0.5,'y(1)--实线  y(2)--点划线  y(3)--虚线');
%-----------------------------------------------------------------------
size(y1)
size(y2)
size(y3)
size(y4)
end