multilinear4ad.m

来自「计算动力学系统的分岔图」· M 代码 · 共 84 行

function ytayl4 = multilinear4AD(func,q1,q2,q3,q4,x0,p1,n) 
%
%----------------------------------------------------
% This file computes  D^(n)*q1*q2*q3*q4 using automatic differentiation, where 
%D^(n) is the 4th order derivatives  of map f^(n). 
%----------------------------------------------------



taylorder=4; 
if (q1==q2)
    if (q1==q3)
        if (q1==q4)
            hc=q1;             
                y1 = Bvv(func,x0,q1,p1,taylorder,n);             
             
        else
           
                y14 = Bvv(func,x0,3.0*q1+q4,p1,taylorder,n);            
             
                y24 = Bvv(func,x0,3.0*q1-q4,p1,taylorder,n);
             
                y34 = Bvv(func,x0,q1+q4,p1,taylorder,n);
           
                y44 = Bvv(func,x0,q1-q4,p1,taylorder,n);
             
               y1=1/192.0*(y14-y24-3*y34+3*y44);
        end        
    else
        
           y61 = Bvv(func,x0,2*q1+q3+q4,p1,taylorder,n);
         
           y62 = Bvv(func,x0,2*q1+q3-q4,p1,taylorder,n);
         
           y63 = Bvv(func,x0,2*q1-q3+q4,p1,taylorder,n);
         
           y64 = Bvv(func,x0,2*q1-q3-q4,p1,taylorder,n);
         
           y65 = Bvv(func,x0,q3+q4,p1,taylorder,n);
         
           y66 = Bvv(func,x0,-q3+q4,p1,taylorder,n);
        
          y1=1/192.0*(y61-y62-y63+y64-2*y65+2*y66);
        
    end    
else
    
    y81 = Bvv(func,x0,q1+q2+q3+q4,p1,taylorder,n);
    
    y82 = Bvv(func,x0,q1+q2+q3-q4,p1,taylorder,n);
    
    y83 = Bvv(func,x0,q1+q2-q3+q4,p1,taylorder,n);
    
    y84 = Bvv(func,x0,q1+q2-q3-q4,p1,taylorder,n);    
   
    y85 = Bvv(func,x0,q1-q2+q3+q4,p1,taylorder,n);
    
    y86 = Bvv(func,x0,q1-q2+q3-q4,p1,taylorder,n);
   
    y87 = Bvv(func,x0,q1-q2-q3+q4,p1,taylorder,n);
   
    y88 = Bvv(func,x0,q1-q2-q3-q4,p1,taylorder,n);
    
    y1=1/192.0*(y81-y82-y83+y84-y85+y86+y87-y88); 
    
end      
ytayl4=24*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