magnus4.m

用magnus展开求解延迟微分方程！这是一种保结构算法。通过两次转换得到的！

function Magnus4(h,N)
%%%%%%%%%%%%%%%%%%%%求解延迟微分方程的四阶Magnus方法程序

%h=0.01;N=200;
%h=0.1;N=80;
c1=(1/2-sqrt(3)/6)*h;
c2=(1/2+sqrt(3)/6)*h;
t0=0;
Y(:,1)=[0 0 0 1 0 0 2 1]'; %%%%%数值解初始向量值
for n=1:N
    t(n)=t0+n*h;
    k=n;
    A1=A(t(n)+c1);
    A2=A(t(n)+c2);
    U1=A1+A2;
    U2=A1-A2;
    U3=A1*A2-A2*A1;
    B=1/2*h*U1-sqrt(3)/12*h^2*U2+1/80*h^3*(U2*U3-U3*U2);
    Y(:,n+1)=expm(B)*Y(:,n);
   M(:,n+1)=[exp(-3*t(n)^2)+exp(-2*t(n));exp(-2*t(n))];
end
V1=Y(7,:);%%%%%%%%%%%数值解的第一个分量
V2=Y(8,:);%%%%%%%%%%%数值解的第二个分量
for n=1:N+1
    t(n)=t0+(n-1)*h;
   M1(n)=exp(-3*t(n)^2+exp(-2*t(n));%%%%%%%%%%%精确解的第一个分量
   M2(n)=exp(-2*t(n));%%%%%%%%%%%%%%%%%%%%%%%%%精确解的第二个分量
end
t=[0:N]*h;
subplot(2,2,1);
plot(t,Y(7,:),'.k',t,M1,'r');
xlabel('the first element of solution vector');
subplot(2,2,2);
plot(t,Y(8,:),'.k',t,M2,'r');
xlabel('the second element of solution vector');
%%%%%%%%%%%%%%%%%%%%%%%%%求解的常微分方程（经离散得到的）
function A=A(t)
C1=[-1 0 1 0;0 -1 0 1;0 0 -1 0;0 0 0 -1];
C2=[0 0 0 0;0 0 0 0;1 0 0 0;0 1 0 0];
C3=[0 0 0 0;0 0 0 0;exp(-6*t) -exp(-6*t) 0 0;0 1/2 0 0];
C4=[-1 0 1 0;0 -1 0 1;0 0 -exp(-3)-6*t exp(-3)+6*t-2; 0 0 0 -(1/2)*exp(2)-2];
A=[C1,C2;C3,C4];