这个脉冲微分方程的数值模拟程序有什么问题.
function DrawTraj6()PeriodT=1; p=0.55;
NumOfStep=50;NumOfPer=5;
Init_x1=0.3; Init_x2=0.5; Init_t=0;
NStep1=NumOfStep;
x1=zeros((NumOfStep+2)*NumOfPer,1);
x2=zeros((NumOfStep+2)*NumOfPer,1);
t=zeros((NumOfStep+2)*NumOfPer,1);
ppp=odeset('RelTol',1e-8,'AbsTol',1e-8);
for i=1:NumOfPer
if i==1
x1(1)=Init_x1;x2(1)=Init_x2;
else
Init_x1=x1((i-1)*(NumOfStep+2));
Init_x2=(1-p)*x2((i-1)*(NumOfStep+2));
x1((i-1)*(NumOfStep+2)+1)=Init_x1;
x2((i-1)*(NumOfStep+2)+1)=Init_x2;
end
t((i-1)*(NumOfStep+2)+1)=(i-1)*PeriodT+Init_t;
sol=ode45(@OdePP,[(i-1)*PeriodT+Init_t,i*PeriodT+Init_t],[Init_x1,Init_x2],ppp);
j=linspace(PeriodT*(i-1)+Init_t,i*PeriodT+Init_t,NStep1);
x1(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,1);
x2(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,2);
t(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=j;
end
save TrajData
figure(1);
plot(x1,x2,'k');
box on
xlabel('x_1');
ylabel('x_2');
hgsave('Traj');
figure(2)
plot(t,x2,'k');
xlabel('t');
ylabel('x_2');
hgsave('tseries_x1');
function dxdt=OdePP(t,x)
r1=8; b1=0.1; a1=2.5; r2=2; a2=0.58; k1=5; k2=5;
dxdt=[x(1)*(r1-b1*x(1)-a1*x(2)/(x(1)+k1))
x(2)*(r2-a2*x(2)/(x(1)+k2))
];