multilinear3ad.m

计算动力学系统的分岔图

function ytayl3 = multilinear3AD(func,q1,q2,q3,x0,p1,n) 
%
%----------------------------------------------------
% This file computes  C^(n)*q1*q2*q3 using automatic differentiation, where 
%C^(n) is the 3rd order derivatives  of map f^(n). 
%----------------------------------------------------

taylorder=3;
if q1==q2
    if q1==q3        
       y1=Bvv(func,x0,q1,p1,taylorder,n);
    else        
       y1p3=Bvv(func,x0,q1+q3,p1,taylorder,n);
       
       y1m3=Bvv(func,x0,q1-q3,p1,taylorder,n);
       
       y3=Bvv(func,x0,q3,p1,taylorder,n);
       
       y1=1/6.0*(y1p3-y1m3)-1/3.0*y3;
    end
else
      
       y123= Bvv(func,x0,q1+q2+q3,p1,taylorder,n);
       
       y12m3= Bvv(func,x0,q1+q2-q3,p1,taylorder,n);
       
       y1m23= Bvv(func,x0,q1-q2+q3,p1,taylorder,n);  
       
       y1m2m3= Bvv(func,x0,q1-q2-q3,p1,taylorder,n);
       
       y1=1/24.0*(y123-y12m3-y1m23+y1m2m3);
end 

ytayl3=6*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