大家好，想请教个问题，最近在做复杂系统同步时要用ddesd解变时滞的微分方程

复杂系统同步时要用ddesd解变时滞的微分方程




function ddeex(t)
t0 = 0;
tfinal = 5;
tspan = [t0, tfinal];
sol = ddesd(@ddefun,@delay,@hist,tspan);%主函数和绘图函数
plot(sol.x,sol.y,'.-')
xlabel('time')
ylabel('value')
legend('x','y')


function d=delay(t,y)  %变时滞函数

d=t-1/(1+exp(-t));


function x=hist(t)  %微分方程注组初值
 s=[0.4;0.6];
n=2;N=7;
x=rand(n*N,1);
x=[s;x];



function [dydt] = ddefun(t,y,Z,p)  %微分方程的函数
n=2; 
a13=1;a21=1;a32=1;a37=-1;a34=-1;a45=1;a47=1;a65=1;a76=1;a71=-1;
A=[2 -0.1;-5 1.5];C=[1 0;0 1];B=[-1.5 -0.1;-0.2 -1];
X=hist(t);s=X(1:2);x=X(3:16);
x1=x(1:n);x2=x(n+1:2*n);x3=x(2*n+1:3*n);x4=x(3*n+1:4*n);x5=x(4*n+1:5*n);x6=x(5*n+1:6*n);x7=x(6*n+1:7*n);  %因为微分方程中含有矩阵，所以通过此语句将微分方程组初值分别导入对应的微分方程中

u3=-20*25*(x3-s);     %控制器
dsdt=-C*s+A*f(s)+B*f(Z(1:2,1));   %领导节点状态函数的微分方程
dxdt=[
      
    -C*x1+A*f(x1)+B*f(Z(1:2,1))-25*(abs(a21)*(x1-sign(a21)*x2))-25*(abs(a71)*(x1-sign(a71)*x7)); %1到7跟随节点的状态函数的微分方程
    -C*x2+A*f(x2)+B*f(Z(5:6,1))-25*(abs(a32)*(x2-sign(a32)*x3));
    -C*x3+A*f(x3)+B*f(Z(1:2,1))-25*(abs(a13)*(x3-sign(a13)*x1))+u3;                                                 
     -C*x4+A*f(x4)+B*f(Z(9:10,1))-25*(abs(a34)*(x4-sign(a34)*x3));
     -C*x5+A*f(x5)+B*f(Z(11:12,1))-25*(abs(a45)*(x5-sign(a45)*x4))-25*(abs(a65)*(x5-sign(a65)*x6));
     -C*x6+A*f(x6)+B*f(Z(13:14,1))-25*(abs(a76)*(x6-sign(a76)*x7));
     -C*x7+A*f(x7)+B*f(Z(15:16,1))-25*(abs(a37)*(x7-sign(a37)*x3))-25*(abs(a47)*(x7-sign(a47)*x4));
    
];


dydt=[dsdt;dxdt];
function G=f(x)    %f函数
G=tanh(x);