// ***  duffing.sci  ***
// Determining equation for Duffing equation
// by Shin'ichi OISHI
// 1999 12.19
function r=duf(v)
   k=0.1;
   B=0.1;

   ss=v('dim');
   s=ss(1)-1;
   x=fourier(s);
   x('cos')=mat(s,1);
   x('sin')=mat(s,1);

   s1=ss(1)/2;
   for i=1:s1,
      x('cos')(i*2+1)=v(i*2+1);
      x('sin')(i*2+1)=v(i*2+2);
   end;

   y=fourier_dif(fourier_dif(x)) + k * fourier_dif(x) +x*x*x;
   y('cos')(3)=y('cos')(3)-B;

   r=mat(ss(1),1);

   for i=1:s1,
      r(i*2+1)=y('cos')(i*2+1);
      r(i*2+2)=y('sin')(i*2+1);
   end;
endfunction
     
function [x,y]=duf2(v)
   k=0.1;
   B=5;

   ss=v('dim');
   s=ss(1)-1;
   x=fourier(s);
   x('cos')=mat(s,1);
   x('sin')=mat(s,1);

   s1=ss(1)/2;
   for i=1:s1,
      x('cos')(i*2+1)=v(i*2+1);
      x('sin')(i*2+1)=v(i*2+2);
   end;

   y=fourier_dif(fourier_dif(x)) + k * fourier_dif(x) +x*x*x;
   y('cos')(3)=y('cos')(3)-B;
endfunction

function [acp,p,c,C1,r_n]=duffing_est(v)
   k=0.1;
   B=5;
   down();
   ppp = 2* (%pi - 0.000000000000001);
   qq = tanh(ppp) - 1e-14;
   up();
   bb=sqrt(2 *(1+k*k)); 
   pp = 2* (%pi + 0.000000000000001);
   dd = sqrt(pp/qq) + 1e-14;

   ss=v('dim');
   s1=ss(1);
   eint=i(-1e-15,1e-15);
   for j=1:s1, v(j+2)=v(j+2)+eint; end; 
   y=init_dif(v);
   zz=duf(y);
   [jl,ju]=int_jacobi(zz);

   R=inv(ju);

   nor1 = norm(R,'fro');
   jc=(jl+ju)/2;
   jr=jc-jl;
   Gu=R*jc-eye(s1,s1)+abs(R)*jr;
   down();
   Gd=R*jc-eye(s1,s1)+(-abs(R))*jr;
   up();
   Gmax=max(abs(Gu),abs(Gd));
   nG = norm(Gmax, 'fro');
   down();
   dnor = 1-nG; 
   up();
   nor = nor1/dnor;
   [x,y] = duf2(v);
   r_n = int_fourier_norm(y);
   k1 = 3 * int_fourier_maxnorm(x*x);
   down();
   cd = sqrt((s1+1)^4 + k^2 * (s1+1)^2);
   up();
   c = 1/cd;
   down();
   C1 = 1+(-c)*k1*(1+ nor * k1);
   up();   
   C=( (1+c) * (1+ nor * k1)+nor) /C1;
   p= 2 * C * r_n;
   a = 6 * bb * dd * (int_fourier_norm(x) + p)/(sqrt(2)-1e-14);
   acp = a * C *p;

endfunction