脉冲微分方程的数值模拟
function [y1,y2,y3]=OdePP() %画函数图像PeriodT=0.5;p=2; %脉冲周期、脉冲量
NumOfStep=50;NumOfPer=4; %NumOfStep是每半个周期取多少点画图,NumOfPer 是画多少个周期的相图
Init_y1=1*10^6;Init_y2=1*10^3;Init_y3=0;Init_t=0; %初值和初始时刻
NStep1=fix(NumOfStep*1);
NStep2=NumOfStep-NStep1;
y1=zeros((NumOfStep+2)*NumOfPer,1);%存放第一维数据的数组
y2=zeros((NumOfStep+2)*NumOfPer,1);%存放第二维数据的数组
y3=zeros((NumOfStep+2)*NumOfPer,1);%存放第三维数据的数组
t=zeros((NumOfStep+2)*NumOfPer,1);%存放时间值数据的数组
ppp=odeset('RelTol',1e-1,'AbsTol',1e-1); %控制误差的参数
for i=1:NumOfPer %循环一次画一个周期的轨线
if i==1
y1(1)=Init-y1;y2(1)=Init-y2;y3(1)=Init-y3 ;
else
%下面3个语句确定在T时刻的脉冲值
Init-y1=y1*((i-1)*(NumOfStep+2));
Init-y2=y2*((i-1)*(NumOfStep+2));
Init-y3=y3*((i-1)*(NumOfStep+2))*(1-p);
y1*((i-1)*(NumOfStep+2)+1)=Init-y1;
y2*((i-1)*(NumOfStep+2)+1)=Init-y2;
y3*((i-1)*(NumOfStep+2)+1)=Init-y3;
end
t((i-1)*(NumOfStep+2)+1)=(i-1)*PeriodT+Init-t;
sol=ode45(@OdePP,[(i-1)*PeriodT+Init-t,(i-1+1)*PeriodT+Init-t],[Init-y1,Init-y2,Init-y3],ppp);
j=linspace(PeriodT*(i-1)+Init-t,(i-1+1)*PeriodT+Init-t,NStep1);
y1(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,1);
y2(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,2);
y3(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=deval(sol,j,3);
t(((i-1)*(NumOfStep+2)+2):((i-1)*(NumOfStep+2)+1+NStep1))=j;
%以上第一次脉冲完成
end
save TrajData %此语句用于保存所有内存数据
figure(1); %此语句绘出相图
plot3(y1,y2,y3,'k');
box on xlabel('y-1');ylabel('y-2');zlabel('y-3');
hgsave('Traj');%此语句用于保存绘出的图
figure(2);
plot(t,y1,'k');
xlabel('t');ylabel('y-1') ;
hgsave('tseries-y1');%绘制并保存第一维的时间序列图
function dydt=OdePP(t,y) %决定方程的函数
r=1; k=3;a=3;b=4; d=0.2;p=0.4 ; q=0.45; e=0.2; %系统参数
dy = zeros(3,1);
dy(1) =r*y(1)*(1-y(1)/k)-a*y(1)*y(2);
dy(2) =y(2)*((b-1)*a*y(1)-d-p*y(3));
dy(3) =y(3)*(q*y(2)-e);
请问谁做过脉冲微分方程的数值模拟?以上是我根据网上的资源编的脉冲程序,不过有问题运行不出来,请有经验的高手帮我看看,这里先谢谢啦!如果成功的话,有重酬!
[此贴子已经被作者于2018-7-14 17:41编辑过]