脉冲时滞微分方程matlab方程

function maichon_shizhi1()

clear;clc

PeriodT=10;q=0.5;

%NumOfStep 是每半个周期取多少点画图，NumOfPer 是画多少个周期的相图

NumOfStep=50;NumOfPer=10;

%以下是初值和初始时刻

Inx1=1;Inx2=1;Inx3=1;Inxt=0;

for j=1:20

x1=zeros((NumOfStep+1)*NumOfPer,1);

x2=zeros((NumOfStep+1)*NumOfPer,1);

x3=zeros((NumOfStep+1)*NumOfPer,1);

t=zeros((NumOfStep+1)*NumOfPer,1);

temp1=0;temp2=0;temp3=0;

%延迟时间值

lags=[0.5];

%下面是控制误差的参数，越小则运行速度越慢

ppp=odeset('RelTol',1e-8,'AbsTol',1e-8);

for i=1:NumOfPer  %循环一次画一个周期的轨线

    

    if i==1

       x111= Inx1;

       x211= Inx2;

       x311= Inx3;       

    else

  %以下三个语句确定脉冲后的值         

       x111=temp1+2.1;

       x211=temp2+0.5;

       x311=temp3;

       

       Inx1=x111;

       Inx2=x211;

       Inx3=x311;

    end

    sol1 = dde23(@myddefun,lags,[Inx1,Inx2,Inx3],[(i-1)*PeriodT+Inxt,(i-0.5)*PeriodT+Inxt]);

%     

%     sol1.x;

%     sol1.y(1,:);

%     sol1.y(2,:);

%     sol1.y(3,:);   

 

    temp1=sol1.y(1,end);

    temp2=sol1.y(2,end);

    temp3=sol1.y(3,end);

     %以下三个语句确定脉冲后的值    

        Inx1 = 0.5 * temp1;

        Inx2 = temp2 + 5*temp1;

        Inx3 = temp3;

    

       x1((i-1)*(NumOfStep+1)+1)=Inx1;

       x2((i-1)*(NumOfStep+1)+1)=Inx2;

       x3((i-1)*(NumOfStep+1)+1)=Inx3;

       

    t((i-1)*(NumOfStep+1)+1)=(i-0.5)*PeriodT+Inxt;

    sol=ode45(@OdePP,[(i-0.5)*PeriodT+Inxt,i*PeriodT+Inxt],[Inx1,Inx2,Inx3],ppp);

    j=linspace(PeriodT*(i-0.5+1/NumOfStep)+Inxt,i*PeriodT+Inxt,NumOfStep);

    x1(((i-1)*(NumOfStep+1)+2):(i*(NumOfStep+1)))=deval(sol,j,1);

    x2(((i-1)*(NumOfStep+1)+2):(i*(NumOfStep+1)))=deval(sol,j,2);

    x3(((i-1)*(NumOfStep+1)+2):(i*(NumOfStep+1)))=deval(sol,j,3);

    t(((i-1)*(NumOfStep+1)+2):(i*(NumOfStep+1)))=j;

    

    t1=t(((i-1)*(NumOfStep+1)+1):(i*(NumOfStep+1)));

    x11=x1(((i-1)*(NumOfStep+1)+1):(i*(NumOfStep+1)));

    x21=x2(((i-1)*(NumOfStep+1)+1):(i*(NumOfStep+1)));

    x31=x3(((i-1)*(NumOfStep+1)+1):(i*(NumOfStep+1)));

    

    temp1=x11(end,1);

    temp2=x21(end,1);

    temp3=x31(end,1);

    

    t1_zhuanzhi= (t1).';

    x11_zhuanzhi=(x11).';

    x21_zhuanzhi=(x21).';

    x31_zhuanzhi=(x31).';

    

   

%     figure(1)      

%

%     plot(sol1.x,sol1.y(1,:),(t1).',x11_zhuanzhi);     

%     ylabel('S1(t)');

%     hold on

%     axis([(i-1)*PeriodT,i*PeriodT,0,2]);

%

%     figure(2)    

%     plot(sol1.x,sol1.y(2,:),(t1).',x21_zhuanzhi);

%     ylabel('S2(t)');

%     hold on

%     figure(3)    

%     plot(sol1.x,sol1.y(3,:),(t1).',x31_zhuanzhi);

%     ylabel('x(t)');

%     hold on

    figure(4)

    plot3(horzcat(sol1.x,t1_zhuanzhi),horzcat(sol1.y(1,:),x11_zhuanzhi),horzcat(sol1.x,t1_zhuanzhi),horzcat(sol1.y(3,:),x31_zhuanzhi),horzcat(sol1.x,t1_zhuanzhi),horzcat(sol1.y(2,:),x21_zhuanzhi));

    xlabel('S1(t)');

    ylabel('x(t)');

    zlabel('S2(t)');

    hold on

    

%     figure(5)

%     plot3(horzcat(sol1.y(1,:),x11_zhuanzhi),horzcat(sol1.y(3,:),x31_zhuanzhi),horzcat(sol1.y(2,:),x21_zhuanzhi));

% %     xlabel('S1(t)');

% %     ylabel('x(t)');

% %     zlabel('S2(t)');

%     hold on

end

Inx1=Inx1+0.5;

Inx2=Inx2+0.3;

Inx3=Inx3+0.2;

end

function dx= myddefun(t,x,Z)

dx=[

    2.5-x(1);

    -x(2)-(x(2)*x(3))/(0.5*(1+x(2)));

    -x(3)+exp(-0.5)*Z(2,1)*Z(3,1)/(1+Z(2,1))];

function dxdt=OdePP(t,x)

%以下参数是系统参数

dxdt=[ 2.5-x(1)

       -x(2)

       -0.5*x(3) ];