portraitzzy.m

来自「同上」· M 代码 · 共 25 行

function [DX1,DX2,DP]=portraitzzy(x1,x2,h)
% PORTRAITZZY		采用“一步仿真”计算状态变量斜率和状态导数的二次方根
% x1,x2				分别给定“状态平面”上的格点坐标
% h					给定积分计算采用的时间步长
% DX1,DX2			轨线斜率在状态坐标轴上的投影长度
% DP				状态导数向量的长度

opts=simset('solver','ode5','fixedstep',h);
n=length(x1);			
X1=zeros(n,n);X2=X1;
for ii=1:n;
   for jj=1:n;
      opts=simset(opts,'initialstate',[x1(ii),x2(jj)]);
      [t,x,y]=sim('exm080634_1',h,opts);	
      dx1=x(2,1)-x1(ii);
      dx2=x(2,2)-x2(jj);
      L=sqrt(dx1^2+dx2^2);
      Z(jj,ii)=L;
      if L>1.e-10								
         DX1(jj,ii)=dx1;DX2(jj,ii)=dx2;	         												
      end
   end
end
DP=Z/h;