lesliegower.m

计算动力学系统的分岔图

% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,b1,b2,ma,mj,my,mz,cjj,cja,cyy,cyz,cjy,cjz,cyj,cya)
jacp=[[kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4)),0,0,0,0,0,-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(1),-b1*kmrgd(2)^2/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2,0,0,-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(3),-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(4),0,0];[0,0,-kmrgd(2),-kmrgd(1),0,0,0,0,0,0,0,0,0,0];[0,kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4)),0,0,0,0,0,0,-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(3),-b2*kmrgd(4)^2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2,0,0,-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(1),-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(2)];[0,0,0,0,-kmrgd(3),-kmrgd(4),0,0,0,0,0,0,0,0]];
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,b1,b2,ma,mj,my,mz,cjj,cja,cyy,cyz,cjy,cjz,cyj,cya)
hess1=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjj^2,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjj+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja*cjj,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*cjj,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cjj];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyj^2,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya*cyj,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*cyj,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyj+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cyj];[0,0,0,0]];
hess2=[[-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjj+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja*cjj,-2*b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cja+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja^2,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjy+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*cja,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjz+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cja];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya*cyj,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya^2,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*cya,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cya+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cya];[0,0,0,0]];
hess3=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*cjj,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjy+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*cja,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy^2,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cjy];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*cyj,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*cya,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy^2,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyy+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cyy];[0,0,0,0]];
hess4=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cjj,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjz+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cja,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*cjy,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz^2];[0,0,0,0];[-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyj+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cyj,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cya+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cya,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyy+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*cyy,-2*b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyz+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz^2];[0,0,0,0]];
hess(:,:,1) =hess1;
hess(:,:,2) =hess2;
hess(:,:,3) =hess3;
hess(:,:,4) =hess4;
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,b1,b2,ma,mj,my,mz,cjj,cja,cyy,cyz,cjy,cjz,cyj,cya)
hessp1=[[-kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjj,1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))-kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cja,-kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjy,-kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*cjz];[0,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp2=[[0,0,0,0];[0,0,0,0];[-kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyj,-kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cya,-kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyy,1/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))-kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*cyz];[0,0,0,0]];
hessp3=[[0,0,0,0];[0,-1,0,0];[0,0,0,0];[0,0,0,0]];
hessp4=[[0,0,0,0];[-1,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp5=[[0,0,0,0];[0,0,0,0];[0,0,0,0];[0,0,-1,0]];
hessp6=[[0,0,0,0];[0,0,0,0];[0,0,0,0];[0,0,0,-1]];
hessp7=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjj*kmrgd(1)-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2,-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(1)+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja*kmrgd(1),2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*kmrgd(1),2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*kmrgd(1)];[0,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp8=[[2*b1*kmrgd(2)^2/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjj,-2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2+2*b1*kmrgd(2)^2/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja,2*b1*kmrgd(2)^2/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy,2*b1*kmrgd(2)^2/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz];[0,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp9=[[0,0,0,0];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyj*kmrgd(3),2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya*kmrgd(3),2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*kmrgd(3)-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2,-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(3)+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*kmrgd(3)];[0,0,0,0]];
hessp10=[[0,0,0,0];[0,0,0,0];[2*b2*kmrgd(4)^2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyj,2*b2*kmrgd(4)^2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya,2*b2*kmrgd(4)^2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy,-2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2+2*b2*kmrgd(4)^2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz];[0,0,0,0]];
hessp11=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjj*kmrgd(3),-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(3)+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja*kmrgd(3),2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*kmrgd(3)-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2,2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*kmrgd(3)];[0,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp12=[[2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjj*kmrgd(4),-b1/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2*kmrgd(4)+2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cja*kmrgd(4),2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjy*kmrgd(4),2*b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^3*cjz*kmrgd(4)-b1*kmrgd(2)/(1+cjj*kmrgd(1)+cja*kmrgd(2)+cjy*kmrgd(3)+cjz*kmrgd(4))^2];[0,0,0,0];[0,0,0,0];[0,0,0,0]];
hessp13=[[0,0,0,0];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyj*kmrgd(1)-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya*kmrgd(1),2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*kmrgd(1),-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(1)+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*kmrgd(1)];[0,0,0,0]];
hessp14=[[0,0,0,0];[0,0,0,0];[2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyj*kmrgd(2),2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cya*kmrgd(2)-b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2,2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyy*kmrgd(2),-b2/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^2*kmrgd(2)+2*b2*kmrgd(4)/(1+cyj*kmrgd(1)+cya*kmrgd(2)+cyy*kmrgd(3)+cyz*kmrgd(4))^3*cyz*kmrgd(2)];[0,0,0,0]];