mainnark4v.m

自己编写的变步长龙格库塔方法解强非线性微分方程

clear;
clc
% % % global e kesi w 
% % % e=0;
% % % kesi=0.01;
% % % 
% % % y0=[0;0];
% % % t0=1;
% % % t1=500;
% % % w=0.56; 

global Tm Ta i12 J1 J2 wa wn kesi aomiga b r1 r2 rn km  wen

k=3.12e7; c=0.1; rn=0.034; r1=0.0768;r2=0.06;b=0.00005;
Tm=4;
Ta=11.45;   % 14.3965  %33.04899989;  % 0.67  0.07  33.4-33.5值得关注    20.5  19.8-19.9
i12=1.28;
J1=0.005;
J2=0.0036;
wa=20000;
n=600;
we=n*pi;  %大轮齿数z＝30  f=2*n*pi/60 / 2*pi/z   w=2*pi*f

km=k;
c12=c;
me=J1*J2/(r2*r2*J1+r1*r1*J2);
wn=sqrt(km/me);
wen=we/wn;
kesi=0.04;

% 
%  for nn=1:50
%      
%  Ta=15+0.1*nn;
%we=(200+10*nn)/60*2*pi*pi/0.1;
% wen=we/wn;
%wen=1.059;
wen=0.7;


t0=0;
t1=1000;
y0=[1;0];

format long;
[t,y]=nark4v([t0 t1],y0,1e-8);%时间历程图
% % % [t,y]=rkf45v([t0 t1],y0,1e-12);%时间历程图

plot(t,y(:,1));
figure
plot(y(3000:end,1),y(3000:end,2));%相图


maxk=round(max(t/2/pi*wen));%poincare 截面
 k=1;
    while (k-1)*2*pi/wen+200<t1
    d=t-(k-1)*2*pi/wen-200;
    [z,n]=sort(abs(d));
    tl=t(n(1));
    tr=t(n(2));
    y1l=y(n(1),1);
    y1r=y(n(2),1);
    y2l=y(n(1),2);
    y2r=y(n(2),2);
   if abs(tl)+abs(tr)<=1e-15
       yy1(k)=y1l;
       yy2(k)=y2l;
   else
       yy=rk4(tl,[y1l;y2l],(k-1)*2*pi/wen+200-tl);
       yy1(k)=yy(1);
       yy2(k)=yy(2);
       
%       Q=polyfit([tl,tr],[y1l,y1r],1); 
%       yy1(k)=polyval(Q,(k-1)*2*pi/w+200); 
%       Q=polyfit([tl,tr],[y2l,y2r],1); 
%       yy2(k)=polyval(Q,(k-1)*2*pi/w+200); 
  end 
  k=k+1;
end
  
figure
axis([-4,4,-4,6]);
plot(yy1,yy2,'r.');
    
format short;