multilinear5ad.m

计算动力学系统的分岔图

function ytayl5 = multilinear5AD(func,q1,q2,q3,q4,q5,x0,p1,n) 
%
%----------------------------------------------------
% This file computes  E^(n)*q1*q2*q3*q4*q5 using automatic differentiation, where 
%C^(n) is the 5th order derivatives  of map f^(n). 
%----------------------------------------------------


taylorder=5;

if (q1==q2 & q1 ==q3)
    if (q1==q4 & q1==q5)

     y1 = Bvv(func,x0,q1,p1,taylorder,n);
 
    else
                
        y16 = Bvv(func,x0,3.0*q1+2.0*q4,p1,taylorder,n);                
              
        y26 =  Bvv(func,x0,3.0*q1-2.0*q4,p1,taylorder,n);
        
        y36 =  Bvv(func,x0,3.0*q1,p1,taylorder,n);
        
        y46 = Bvv(func,x0,q1,p1,taylorder,n);
        
        y56 = Bvv(func,x0,q1+2.0*q4,p1,taylorder,n);
        
        y66 = Bvv(func,x0,q1-2.0*q4,p1,taylorder,n);
        
        y1 = (y16 + y26 - 2.0*y36 + 6.0*y46 - 3.0*y56 - 3.0*y66)/1920.0;
    end
    else        
                
         y18 =Bvv(func,x0,q1+q2+q3+q4+q5,p1,taylorder,n);    
        
         y28 = Bvv(func,x0,q1+q2+q3+q4-q5,p1,taylorder,n);
        
         y38 = Bvv(func,x0,q1+q2+q3-q4-q5,p1,taylorder,n);
        
         y48 = Bvv(func,x0,q1+q2+q3-q4+q5,p1,taylorder,n);      
        
         y58 = Bvv(func,x0,q1+q2-q3+q4+q5,p1,taylorder,n);
        
         y68 = Bvv(func,x0,q1+q2-q3+q4-q5,p1,taylorder,n);
        
         y78 = Bvv(func,x0,q1+q2-q3-q4-q5,p1,taylorder,n);
        
         y88 = Bvv(func,x0,q1+q2-q3-q4+q5,p1,taylorder,n);
        
         y1= (y18 - y28 + y38 - y48 - y58 + y68 - y78 + y88)/1920.0;
        
         y18 =  Bvv(func,x0,q1-q2+q3+q4+q5,p1,taylorder,n);
        
         y28 =  Bvv(func,x0,q1-q2+q3+q4-q5,p1,taylorder,n);
        
         y38 =  Bvv(func,x0,q1-q2+q3-q4-q5,p1,taylorder,n);
        
         y48 =  Bvv(func,x0,q1-q2+q3-q4+q5,p1,taylorder,n);
        
         y58 =  Bvv(func,x0,q1-q2-q3+q4+q5,p1,taylorder,n);
        
         y68 = Bvv(func,x0,q1-q2-q3+q4-q5,p1,taylorder,n);
        
         y78 =  Bvv(func,x0,q1-q2-q3-q4-q5,p1,taylorder,n);
        
         y88 = Bvv(func,x0,q1-q2-q3-q4+q5,p1,taylorder,n);
      
         y1 = y1 - (y18 - y28 + y38 - y48 - y58 + y68 - y78 + y88)/1920.0;
        
end
ytayl5=120*tcs(y1);
%--------------------------------------------------------------------------
function y1 = Bvv(mapsf,x0,hc,p1,taylorder,n)
   
   % Convert to "active" variables:
   s = adtayl(0,taylorder); %Base point & Taylor order
   y1= x0 + s.*hc;
   for i=1:n     
     y1 = mapsf(0, y1,p1{:});
   end