portraitzzy.m

这是本人收集的一些程序源代码

function [DX1,DX2,DP]=portraitzzy(x1,x2,h)
% PORTRAITZZY		采用"一步仿真"计算状态变量斜率和状态导数的二次方根
% x1,x2				分别给定"状态平面"上的格点坐标
% h						给定积分计算采用的时间步长
% DX1,DX2				轨线斜率在状态坐标轴上的投影长度
% DP					状态导数向量的长度
opts=simset('solver','ode5','fixedstep',h);		%采用ode5定步长积分算法	<7>
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)]);%设置状态初值		<12>
      [t,x,y]=sim('exm09634_1',h,opts);	%步长为 h 的"一步仿真"		<13>
      dx1=x(2,1)-x1(ii);					%计算x1的变化率
      dx2=x(2,2)-x2(jj);					%计算x2的变化率
      L=sqrt(dx1^2+dx2^2);					%计算状态轨线长度的变化率。
      Z(jj,ii)=L;
      if L>1.e-10								%若状态轨线变化率大于"零"阈值
         DX1(jj,ii)=dx1/L;DX2(jj,ii)=dx2/L;	%计算各状态变量的近似斜率	<19>
         										%注意下标次序。这是绘图指令格式要求		
      end
   end
end
DP=Z/h;										%状态导数向量的长度