Return to taka_pfp.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Fixed a bug for a complex step size.

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_pfp.rr,v 1.11 2003/03/11 02:00:53 takayama Exp $ */

/* Table. It should be moved to ??? in a fugure. Experimental */
#define Xm_eval_in_quote  0
#define Xm_eval_in_series 0x100
#define Xm_eval_in_num    0x200
/* minor numbers */
#define Xm_eval_mask       0xff
#define Xm_by_series       0
#define Xm_by_int          1
#define Xm_by_runge_kutta  2
#define Xm_by_conn_formula 3

                        /* _need space */
Xm_eval = Xm_eval_in_num  $

def taka_pfp_hypergeometric_pfq(A,B,Z) {
  extern Xm_eval;
  if (iand(Xm_eval,Xm_eval_in_series)) {
    return taka_pfp(A,B,Z,10); /* first ten terms, experimental */
  }
  if (iand(Xm_eval,Xm_eval_in_num)) {
    if (iand(Xm_eval,Xm_eval_mask) == Xm_by_int) {
       if (length(A) == 2 && length(B) == 1) {
          return tam_int1_2f1(A[0],A[1],B[0],Z);
       }
    }
    if (iand(Xm_eval,Xm_eval_mask) == Xm_by_runge_kutta) {
       return taka_pfp_runge_kutta(A,B,Z);
    }
    if (iand(Xm_eval,Xm_eval_mask) == Xm_by_conn_formula) {
       return tam_conn_pfq(A,B,Z);
    }
    /* Else use series */
    return taka_pfp(A,B,Z,10);
  }
  F=base_replace(quote(hypergeometric_pfp(a,b,z)),
                 [[a,A],[b,B],[z,Z]]);
  return F;
}

def taka_pfp_hypergeometric_2f1(A,B,C,Z) {
  return taka_pfp_hypergeometric_pfq([A,B],[C],Z);
}

def taka_pfp_gamma(A) {
  extern Xm_eval;
  if (Xm_eval == 0) {
    F=base_replace(quote(hypergeometric_gamma(a)),
                 [[a,A]]);
    return F;
  }else{
    return(pari(gamma,A));
  }
}

/* They have not yet been registered in names.rr */

/* Codes from misc/200205/runge-kutta.rr */
def taka_pfp_poch(A,B,N) {
  R = 1;
  M1 = length(A);
  M2 = length(B);
  for (I=0; I<N; I++) {
    for (J=0; J<M1; J++) {
      R = R*(A[J]+I);
    }
    for (J=0; J<M2; J++) {
      R = R/(B[J]+I);
    }
    R = R/(I+1);
  }
  return(R);
}
def taka_pfp_poch1(A,B,N) {
  R = 1;
  M1 = length(A);
  M2 = length(B);
  for (I=0; I<N; I++) {
    for (J=0; J<M1; J++) {
      R = R*(A[J]+I);
    }
    for (J=0; J<M2; J++) {
      R = R/(B[J]+I);
    }
  }
  return(R);
}

def taka_pfp_poch2(A,B,N,R) {
  if (N == 0) return R;
  M1 = length(A);
  M2 = length(B);
  for (J=0; J<M1; J++) {
      R = R*(A[J]+N-1);
  }
  for (J=0; J<M2; J++) {
      R = R/(B[J]+N-1);
  }
  R = R/N;
  return(R);
}

def taka_pfp(A,B,Z,N) {
  R = 0; P = 1;
  for (I=0; I<N; I++) {
     P = taka_pfp_poch2(A,B,I,P);
     R = R+P*Z^I;
  }
  return(R);
}

/*&usage  pfp_euler_derivates(A,B,Z,N,M)
  It returns
  (F,T F, T^2 F, ..., T^M F) where T is the euler operator z d/dz.
 example: taka_pfp_euler_derivatives([1/12,5/12],[1/2],0.4,20,2);
 example: taka_pfp_euler_derivatives([1/12,5/12],[1/2],
             eval(exp(0)*1323/1331),300,2);
*/
/*  test data.
z*D[z*D[ Hypergeometric2F1[1/12,5/12,1/2,z],z],z] /. z->0.4
--> 0.07856325931358881
z*D[z*D[ Hypergeometric2F1[1/12,5/12,1/2,z],z],z] /. z->1323/1331
--> 1992.8936...  (The error is about 1992-1073.)
*/
def taka_pfp_euler_derivatives(A,B,Z,N,M) {
  Ans = newvect(M+1);
  for (I=0; I<=M; I++) Ans[I] = 0;
  P = 1;
  for (I=0; I<N; I++) {
     P = taka_pfp_poch2(A,B,I,P);
     Zi = Z^I;
     C = 1;
     for (J=0; J<=M; J++) {
       Ans[J] += P*Zi*C;  /* C = I^J */
       C   *= I;  
     }
  }
  return(vtol(Ans));
}

def taka_pfp_series(A,B,Z,N) {
  F = taka_pfp(A,B,Z,N);
  F1 = nm(F);
  F2 = dn(F);
  R = [];
  for (I=0; I<N; I++) {
    G1 = coef(F1,I,Z);
    G1 = red(G1/F2);
    R = cons(G1,R);
  }
  return R;
}

/* (th{x}+cc[1])(th{x}+cc[2])(th{x}+cc[3])(th{x}+cc[4])(th{x}+cc[5])  
     - x (th{x}+aa[1])(th{x}+aa[2])(th{x}+aa[3])(th{x}+aa[4])(th{x}+aa[5]) */
/* 
Taka_pfp_rule = [[cc[1],0],[cc[2],c1-1],[cc[3],c2-1],[cc[4],c3-1],[cc[5],c4-1],
                 [aa[1],a1],[aa[2],a2],[aa[3],a3],[aa[4],a4],[aa[5],a5]]$
*/
def taka_pfp_rule(CC,Size) {
  Rule = [[taka_pfp_iv(CC,1),0]];
  for (I=2; I<=Size; I++) {
    Rule = cons([taka_pfp_iv(CC,I),taka_pfp_iv(CC,I-1)-1],Rule);
  }
  return Rule;
}

def taka_pfp_iv(V,I) {
  return strtov(rtostr(V)+rtostr(I));
}

/* Differential equation for p F p-1 is  dY/dx=omega[] Y       */
def taka_pfp_omega(Size) {
     Aaa = taka_pfp_fun(aa,Size); Ccc = taka_pfp_fun(cc,Size); 
     Tmp = newvect(Size);
     for (K=0; K<Size-1; K++) Tmp[K] = taka_pfp_nvec(K+1,Size);
     Tmp[Size-1] = taka_pfp_nnvec(Aaa,Ccc,Size);
     Tmp = map(vtol,Tmp);
     Tmp = vtol(Tmp);
     Tmp = newmat(Size,Size,Tmp);
     return base_replace(Tmp,taka_pfp_rule(cc,Size));
}

def taka_pfp_fun(Bb,Size) {
  Tmp = 1;
  for (K=1; K<=Size; K++) {
    Tmp = Tmp*(x+taka_pfp_iv(Bb,K));
  }
  Ans = newvect(Size);
  for (K=0; K<Size; K++) {
    Ans[K] = coef(Tmp,K,x);
  }
  return Ans;
}

def taka_pfp_nvec(K,Size) {
  Tmp = newvect(Size);
  Tmp[K] = 1/x;
  return Tmp;
}

def taka_pfp_nnvec(Aaa,Ccc,Size) {
  Tmp = newvect(Size);
  for (K=0; K<Size; K++) {
    Tmp[K] = (Aaa[K]*x-Ccc[K])/(x*(1-x));
  }
  return Tmp;
}

/*
 Taka_pfp_rule1 
=[[aa1,1/5], [aa2,1/5], [aa3,1/5], [aa4,1/5],[aa5,1],
  [cc1,2/5], [cc2, 3/5],[cc3,4/5], [cc4,1]]
*/
def taka_pfp_1_omega() {
   Taka_pfp_rule1 
   =[[aa1,1/5], [aa2,1/5], [aa3,1/5], [aa4,1/5],
     [cc1,2/5], [cc2,3/5], [cc3,4/5]]$
   Omega = taka_pfp_omega(4);
   Omega = base_replace(Omega, Taka_pfp_rule1);
   return Omega;
}

/* 
  Test programs.
  Need glib, taka_plot.rr, taka_runge_kutta.rr 
*/
def taka_pfp_1() {
  /* Solving p F p-1 by the adaptive Runge-Kutta. */
   Omega = taka_pfp_1_omega();
   F0 = Omega*newmat(4,1,[[y1],[y2],[y3],[y4]]);
   F = newvect(4,[ F0[0][0], F0[1][0], 
                   F0[2][0], F0[3][0] ]);
   F = base_cancel(F);
   print (F);
   Y = [y1,y2,y3,y4];
   /* T=taka_runge_kutta_4_adaptive(F,x,Y,-10,[1,0,0,0],0.1,-1); */
   T=taka_runge_kutta_4_adaptive(F,x,Y,-10,[1,0,0,0.5],0.1,-0.1); 
   /* T=taka_runge_kutta_4_adaptive(F,x,Y,-0.5,[1,0,0,1],-0.1,-10); */
   /* T1 = number_real_part(T); 
   taka_plot_auto(T1); */ 
   taka_plot_auto(T);
}

def taka_pfp_2_omega() {
   Taka_pfp_rule1 
   =[[aa1,1/2], [aa2,1/2],
     [cc1,3/2]]$
   Omega = taka_pfp_omega(2);
   Omega = base_replace(Omega, Taka_pfp_rule1);
   return Omega;
}

def taka_pfp_2() {
  /* Solving 2 F 1 by the adaptive Runge-Kutta. */
   Omega = taka_pfp_2_omega();
   F0 = Omega*newmat(2,1,[[y1],[y2]]);
   F = newvect(2,[ F0[0][0], F0[1][0]]);
   F = base_cancel(F);
   print (F);
   Y = [y1,y2];
   T=taka_runge_kutta_4_adaptive(F,x,Y,0.1,[1.01746,0.018314],0.001,0.25);
   /* The initial value is computed by Mathematica */
   /* T1 = number_real_part(T); 
   taka_plot_auto(T1); */ 
   print(T[0][1]*3);
   taka_plot_auto(T);
}

def taka_pfp_aa_cc_rule(A,B) {
  Rule = [ ];
  P = length(A);
  Q = length(B);
  for (I=0; I<P; I++) {
    Rule = cons([taka_pfp_iv(aa,I+1),A[I]],Rule);
  }
  for (I=0; I<Q; I++) {
    Rule = cons([taka_pfp_iv(cc,I+1),B[I]],Rule);
  }
  return Rule;
}

/*&usage  pfp_runge_kutta(A,B,Z)
 example: taka_pfp_runge_kutta([1/12,5/12],[1/2],1323/1331); -->exact: 1.3659
 example: taka_pfp_runge_kutta([7/2,7/2],[31/5],1.3+1.3*@i);
             -------> exact: -0.7973029718645974 - 0.5288535335094111*I
 example: taka_pfp_runge_kutta([7/2,7/2],[31/5],13+13*@i);
             -------> exact: 0.001644777051991857 + 0.0005057784138817566*I
*/
#define SERIES_N  30
#define SERIES_LIMIT 0.5
#define INITIAL_STEP 100
def taka_pfp_runge_kutta(A,B,Z) {
  P = length(A); Q = length(B);
  if (P-1 != Q) {
    print([A,B,Z]);
    error("pfp_runge_kutta() does not work for these parameters.");
  }
  if (Z == 1) error("pfp_runge_kutta() cannot evaluate the value at z=1.");
  for (I=0; I<Q; I++) {
    if (type(B[I]) <= 1) {
      if ((ntype(B[I]) == 0) && (dn(B[I]) == 1) && (B[I] < 0)) {
        error("p F q is not defined for parameters in negative integers.");
      }
    }
  }

  /* cf. taka_pfp_1() */
  Omega = taka_pfp_omega(P);
  Ylist1 = [ ];  Ylist0 = [ ];
  for (I=1; I<=P; I++) {
    Ylist1 = cons([taka_pfp_iv("y",I)],Ylist1);
    Ylist0 = cons(taka_pfp_iv("y",I),Ylist0);
  }
  Ylist1 = reverse(Ylist1); Ylist0 = reverse(Ylist0);

  F0 = Omega*newmat(P,1,Ylist1);
  Ylist2 = [ ];
  for (I=0; I<P; I++) {
    Ylist2 = cons(F0[I][0], Ylist2);
  }
  Ylist2 = reverse(Ylist2);

  F = newvect(P,Ylist2);
  F = base_cancel(F);

  Y = Ylist0;

  if (number_abs(Z) < SERIES_LIMIT) {
    return eval(taka_pfp(A,B,eval(exp(0)*Z),SERIES_N));
  }

  if (imag(Z) == 0 && Z > 1) {
     error("pfp_runge_kutta() cannot evalute the value on [1,infty]");
  }
  /* Set the staring point X0 */
  X0 = eval(Z*0.3/number_abs(Z));
  /* Set the initial step */
  H = eval((Z-X0)/INITIAL_STEP);

  InitialValue = taka_pfp_euler_derivatives(A,B,X0,SERIES_N,P-1);
  Rule = taka_pfp_aa_cc_rule(A,B);
  F = base_replace(F,Rule);

  T = taka_runge_kutta_4_adaptive(F,x,Y,X0,InitialValue,H,Z);
  print("Stopped to use the adaptive method. Step size is "+rtostr(H));
  T2 = taka_runge_kutta_4(F,x,Y,T[1][0],cdr(T[1]),H,Z);
  return [[T2[0][0],T2[0][1]], T, T2];  
}

Loaded_taka_pfp=1$
end$