erk4.m

用指数龙格库塔方法求解时滞（延迟）微分方程！

function ERK4(h,tao,N)
%%%%%%%%%%%延迟微分方程的指数龙格库塔方法程序

c(1)=0;c(2)=1/2;c(3)=1;
m=ceil(tao/h);
Is=eye(3);
a=-10;
a0=-5;

syms u;
for i=1:3
    p=1;
    for j=1:3
        if j~=i
       p=p*((u/h)-c(j))/(c(i)-c(j));
    end
   end
   l(i)=p;
end

q=exp((h-u)*a);
for i=1:3
    B(i)=(1/h)*int(q*l(i),u,0,h);
end
for i=1:3
    for j=1:3
       A(i,j)=(1/h)*int(exp((c(i)*h-u)*a)*l(j),u,0,c(i)*h);
    end
end
A=numeric(A);
B=numeric(B);

for i=1:3
C(i)=exp(c(i)*h*a);
end
for n=1:m+1
    k=n-(m+1);
    for i=1:3
    t=k*h+c(i)*h;
    Y(n,i)=sin(t);
end
end
%%%%%%%%%%%延迟微分方程的指数龙格库塔方法
y(m+1)=0;
for n=m+1:N
    for i=1:3
  Y(n+1,i)=exp(h*a)*y(n)+h*A(i,1)*a0*Y(n-m,1)+h*A(i,2)*a0*Y(n-m,2)+h*A(i,3)*a0*Y(n-m,3);
  y(n+1)=exp(h*a)*y(n)+h*B(1)*a0*Y(n-m,1)+h*B(2)*a0*Y(n-m,2)+h*B(3)*a0*Y(n-m,3);
end
 end 
 R=exp(h*a)+h*B*a0*(Is-h*a0*A)^(-1)*C'
t=[1:N+1]*h;
plot(t,y,'k');
title('a= -10,a0= -5,初值为y(t)= sin(t)')