rk_4.m

求解时滞微分方程的龙格库塔方法！用matlab编写的。

function RK_4(h,N)
%%%%%%%%%%%%%% 经典的四阶 Runge-Kutta方法求时滞微分方程
%%%%%%%%%%%%%% 经典的四阶 Runge-Kutta方法的系数
A=[0 0 0 0;1/2 0 0 0;0 1/2 0 0;0 0 1 0];
B=[1/6 1/3 1/3 1/6];
c=[0 1/2 1/2 1];
e=[1,1,1,1]';
t0=0;
m=1/h;
%%%%%%%%%%%%% 初始区间的函数值
for n = 1:m+1
    k = n-(m+1);
    Y(:,n)=[exp(-3*(k*h+c(1)*h)^2),exp(-2*(k*h+c(1)*h)),exp(-3*(k*h+c(2)*h)^2),exp(-2*(k*h+c(2)*h)),exp(-3*(k*h+c(3)*h)^2),exp(-2*(k*h+c(3)*h)),exp(-3*(k*h+c(4)*h)^2),exp(-2*(k*h+c(4)*h))]';
    y(:,n)=[1,1]';
    D(:,n)=[exp(-3*(k*h)^2),exp(-2*(k*h))]';
end
%%%%%%%%%%%%%%% 四阶 Runge-Kutta方法
for n = m+1:N
    t(n)=t0+n*h;
    L=[-exp(-3)-6*t(n),exp(-3)+6*t(n)-2;0,(-1/2)*exp(2)-2];
    M=[exp(-6*t(n)),-exp(-6*t(n));0,1/2];
    Y(:,n)=(eye(8)-h*kron(A,L))^(-1)*(h*kron(e,L)*y(:,n)+h*kron(e,M)*y(:,n-m)+h*kron(A,M)*Y(:,n-m));
    y(:,n+1)=y(:,n)+kron(B,eye(2))* Y(:,n);
    D(:,n+1)=[exp(-3*(t(n))^2),exp(-2*t(n))]';
end
t=[0:N]*h;
subplot(2,2,1);
plot(t,y(1,:),'.k',t,D(1,:),'r');
xlabel('the first element of solution vector');
subplot(2,2,2);
plot(t,y(2,:),'.k',t,D(2,:),'r');
xlabel('the second element of solution vector');