Return to tk_pf2.rr CVS log

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

mtg.plot3d and tk_pfn.graph accept the option fit=1.
If it is set to one, (max+min)/2 is moved to the center of the z-axis.

/*
$OpenXM: OpenXM/src/asir-contrib/packages/src/tk_pf2.rr,v 1.11 2010/08/31 06:37:13 takayama Exp $
tk_pf2 is a package for numerical analysis of Pfaffian systems (systems of differential equaitons) in two variables.
Preliminary version was misc-2009/11/fish/ri.rr
tk_pf2.rk2 and tk_pf2.rk2all are main functions of this package which numerically solve dF=AFdx+BFdy
As to examples of using these functions, see functions test*()
cf. @s/2009/11/14-my-note-ri-pfaffian-na.pdf
    @s/2009/12/09-my-note-*
*/
import("glib3.rr")$
import("taka_runge_kutta.rr")$  /* updated */
Taka_Runge_kutta_debug=0$
import("taka_plot.rr")$
import("mt_graph.rr");


module mtg;
localf setGtable;
localf setGtable0;
localf genGtable;
localf genGtable0;
localf f_gtable;
localf staticv;
localf test6a;
localf test6a2;
localf test6b;
localf test6b2;
localf test6c;
localf plot2d;
localf shrink$

static Gtable$
static Gtable_xmin$
static Gtable_ymin$
static Gtable_xh$
static Gtable_yh$
static Gtable_xn$
static Gtable_yn$

def shrink(Dom) {
  X=Dom[0][1]-Dom[0][0];
  Y=Dom[1][1]-Dom[1][0];
  X = X*0.1;
  Y = Y*0.1;
  return [[Dom[0][0]+X,Dom[0][1]-X],
          [Dom[1][0]+Y,Dom[1][1]-Y]];
}

def setGtable(L) {
  Eps = 0.00001;
  N = length(L);
  Y0 = L[0][0][1];
  Xn = 0;
  for (I=0; number_abs(Y0-L[I][0][1])<Eps; I++) {
    Xn++;
  }
  Yn = number_floor(length(L)/Xn+Eps+1);
  Hx = L[1][0][0] -L[0][0][0];
  if (number_abs(Hx) < Eps) Hx=L[2][0][0] -L[1][0][0];
  Hy = L[Xn+1][0][1]-L[0][0][1];
  if (Hx < 0) return setGtable(reverse(L)); 
  setGtable0(L,Xn+1,Yn+1,Hx,Hy);
}

def setGtable0(L,Xn,Yn,Hx,Hy) {
/*  
  extern Gtable;
  extern Gtable_xmin;
  extern Gtable_ymin;
  extern Gtable_xh;
  extern Gtable_yh;
  extern Gtable_xn;
  extern Gtable_yn;
*/
  Gtable_xn = Xn; Gtable_yn = Yn;
  Gtable=newmat(Gtable_xn,Gtable_yn);
  Gtable_xmin=L[0][0][0];
  Gtable_ymin=L[0][0][1];
  Gtable_xh =Hx;
  Gtable_yh =Hy;
  Eps = 0.00001;
  for (I=0; I<length(L); I++) {
     X=L[I][0][0];
     Y=L[I][0][1];
     P=number_floor((X-Gtable_xmin)/Gtable_xh + Eps);
     Q=number_floor((Y-Gtable_ymin)/Gtable_yh + Eps);
     if (P==0 && Q==0) print(I);
     Gtable[P][Q] = L[I][1];
  }
  return Gtable; 
}

/* setGtable is obsolete now. */
def genGtable(L) {
  Eps = (0.1)^10;
  N = length(L);
  Xmin=Xmax=L[0][0][0];
  Ymin=Ymax=L[0][0][1];
  for (I=1; I<N; I++) {
    X=L[I][0][0];
    Y=L[I][0][1];
    if (X < Xmin) Xmin = X;
    if (X > Xmax) Xmax = X;
    if (Y < Ymin) Ymin = Y;
    if (Y > Ymax) Ymax = Y;
  }
  Hx = Hy=Eps;
  Xold=L[0][0][0]; Yold=L[0][0][1];
  for (I=1; I<N; I++) {
    X=L[I][0][0];
    Y=L[I][0][1];
    Dx= number_abs(X-Xold);
    Dy = number_abs(Y-Yold);
    if ((Dx > Hx) && (Dx < (Xmax-Xmin)/1.2)) Hx = Dx;
    if ((Dy > Hy) && (Dy < (Ymax-Ymin)/1.2)) Hy = Dy;
    Xold=X; Yold=Y;
  }
  Xn = number_floor((Xmax-Xmin)/Hx)+2;
  Yn = number_floor((Ymax-Ymin)/Hy)+2;

/*
  print("domain=",0);print([[Xmin,Xmax],[Ymin,Ymax]]);
  print("Hx,Hy=",0); print([Hx,Hy]);
*/
  if ((Hx <= Eps) || (Hy <= Eps)) {
    error("The data is too coarse. Hx or Hy is too small."); debug();
  }
  genGtable0(L,Xn+1,Yn+1,Hx,Hy,[[Xmin,Xmax],[Ymin,Ymax]]);
}

def genGtable0(L,Xn,Yn,Hx,Hy,Dom) {
/*  
  extern Gtable;
  extern Gtable_xmin;
  extern Gtable_ymin;
  extern Gtable_xh;
  extern Gtable_yh;
  extern Gtable_xn;
  extern Gtable_yn;
*/
  Gtable_xn = Xn; Gtable_yn = Yn;
  Gtable=newmat(Gtable_xn,Gtable_yn);
  Gtable_xmin=Dom[0][0];
  Gtable_ymin=Dom[1][0];
  Gtable_xh =Hx;
  Gtable_yh =Hy;
  Eps = 0.00001;
  for (I=0; I<length(L); I++) {
     X=L[I][0][0];
     Y=L[I][0][1];
     P=number_floor((X-Gtable_xmin)/Gtable_xh + Eps);
     Q=number_floor((Y-Gtable_ymin)/Gtable_yh + Eps);
     if (P==0 && Q==0) print(I);
     if ((P<Xn) && (Q < Yn)) Gtable[P][Q] = L[I][1];
  }
  return Gtable; 
}


def f_gtable(X,Y) {
/*
  extern Gtable;
  extern Gtable_xmin;
  extern Gtable_ymin;
  extern Gtable_xh;
  extern Gtable_yh;
  extern Gtable_xn;
  extern Gtable_yn;
*/
  Eps = 0.00001;
  if (X < Gtable_xmin) return -1;
  if (Y < Gtable_ymin) return -1;
  if (X > Gtable_xmin+Gtable_xh*(Gtable_xn-1)) return -1;
  if (Y > Gtable_ymin+Gtable_yh*(Gtable_yn-1)) return -1;

  Pr=(X-Gtable_xmin)/Gtable_xh;
  Qr=(Y-Gtable_ymin)/Gtable_yh;
  P=number_floor(Pr+Eps);
  Q=number_floor(Qr+Eps);
  if ((P+1 >= Gtable_xn) || (Q+1 >= Gtable_yn)) return Gtable[P][Q];  
  Wx=Pr-P; Wy =Qr-Q;
  if (Wx+Wy<=1) return (Gtable[P][Q]+Wx*(Gtable[P+1][Q]-Gtable[P][Q])
                                    +Wy*(Gtable[P][Q+1]-Gtable[P][Q]));
  else {
     Wx=1-Wx; Wy=1-Wy;
     return (Gtable[P+1][Q+1]+Wx*(Gtable[P][Q+1]-Gtable[P+1][Q+1])
                             +Wy*(Gtable[P+1][Q]-Gtable[P+1][Q+1]));
  }
     
}

def staticv() {
  A=[
   Gtable_xmin,
   Gtable_ymin,
   Gtable_xh,
   Gtable_yh,
   Gtable_xn,
   Gtable_yn,
   Gtable];
  return A;
}
def test6a() {
  L=[];
  for (Y=-2; Y<2; Y += 0.1) {
    for (X=-2; X<2; X+= 0.1) {
      L=cons([[X,Y],X^2-Y^2],L);
    }
  }
  setGtable(reverse(L)); 
  mtg.plot3d(quote(f_gtable(x,y)) | domain=[[-1.9,1.9],[-1.9,1.9]]);
}


def test6b() {
  A=tk_pf2.test6(10);  
  mtg.plot3d(quote(f_gtable(x,y)) | domain=[[1,2.5],[1,2.5]]);
}

def test6b2(Y) {
  AA=[];
  for (X=1; X<2.5; X += 0.2) {
    AA=cons([X,f_gtable(X,Y)],AA);
  }
  taka_plot_auto(reverse(AA));
}

def test6c() {
  for (X=1; X<2.5; X += 0.1) {
   for (Y=0.1; Y<2.5; Y += 0.1) {
      if (number_abs(f_gtable(X,Y)-(X^2-Y^2))> 0.2)
         print([X,Y,f_gtable(X,Y)]);
   }
  }
}

/*
  tk_pf2.test6(10)$
  mtg.plot2d(quote(mtg.f_gtable(x,2)) | domain=[1,2]);
  mtg.plot2d(x^2-2^2 | domain=[1,2]);

  mtg.plot2d(quote(imag(eval((-2-x*@i)^(1/2))) ) | domain=[-2,2]);
  If no "eval", it does not work. There still be a parse error --> bug? 
*/
def plot2d(F) {
  Opt = getopt();
  if ((type(F) == 4) || (type(F)==5)) return taka_plot_auto(F | option_list=Opt);
  Mtsize=1;
  D=getopt(domain);
  if (type(D) <0) {
    D=[-Mtsize,Mtsize];
  }
  NN=getopt(mesh);
  if (type(NN) < 0) {
    NN=20;
  }
  H = eval(exp(0)*(D[1]-D[0])/NN);

  L=[];
  for (X=D[0]; X<=D[1]; X += H) {   
    if (type(F) == 17) { /* quote */
       Z = base_replace(F,[[x,X]]);
       Z = eval_quote(Z);
    }else{
       Z = subst(F,x,X);
    }
    Z = eval(exp(0)*Z);
    L = cons([X,Z],L);
  }
  taka_plot_auto(L | option_list=Opt);
  return([D,NN]); 
}
endmodule;

module tk_pf2;
localf test1 ;
localf test1rr;
localf rk2x ;
localf test2 ;
localf rk2y ;
localf test3 ;
localf rk2;
localf test5;
localf test5rr;
localf test5b;
localf test5c;
localf rk2all;
localf test6;
localf test6n;
localf test6rr;

/* f=x^2-y^2
  f_xy = 0
  x*f_xx-f_x=0
  y*f_yy-f_y=0
*/
def test1() {
  S=[x^2-y^2, 2*x,-2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,2],[y,0.1]]); /* target 1 */
  print("Value at (2,0.1)",0); print(TT);
  TT=base_replace(S,[[x,2],[y,3]]); /* target 2 */
  print("Value at (2,3)",0); print(TT);
  A=rk2x([[y1,y1/x0,0],[y2,0,y2/x1]],
          [x0,x1],[y0,y1,y2],
          [1,0.1], II , [2,3],0.1);
}

def test1rr() {
  S=[x^2-y^2, 2*x,-2*y];
  II=base_replace(S,[[x,3],[y,0.1]]); /* Initial */
  print("Value at (3,0.1)",0); print(II);
  TT=base_replace(S,[[x,2],[y,0.1]]); /* target 1 */
  print("Value at (2,0.1)",0); print(TT);
  TT=base_replace(S,[[x,1],[y,3]]); /* target 2 */
  print("Value at (1,3)",0); print(TT);
  A=rk2x([[y1,y1/x0,0],[y2,0,y2/x1]],
          [x0,x1],[y0,y1,y2],
          [3,0.1], II , [1,3],0.1);
}

/*
  F=[[y1,y1/x0,0],[y2,0,y2/x1]] Pfaffian for x0 and x1.
  X=[x0,x1];  independent variable
  Y=[y0,y1,y1]; dependent variable
  Xs=[s0,s1];  starting point
  Ys=          starting value of Y
  Xt=[t0,t1];  target point.
  H step
*/
def rk2x(F,X,Y,Xs,Ys,Xt,H) {
  X0=X[0]; X1=X[1];
  S0=Xs[0]; S1=Xs[1];
  T0=Xt[0]; T1=Xt[1];
  
  /* goto  x direction */
  FF = base_replace(F[0],[[X1,S1]]);
  A=tk_rk.taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H);
  T=A[0];
/*  print("Mid point ",0); print(T); */

  /* goto y direction */
  Ys = cdr(T);
  FF = base_replace(F[1],[[X0,T0]]);
  A=tk_rk.taka_runge_kutta_4a(FF,X1,Y,S1,Ys,T1,H);
  T=A[0];
  return cons([T0,T[0]],cdr(T));
}

def test2() {
  S=[x^2-y^2, 2*x,-2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,1],[y,3]]); /* target 1 */
  print("Value at (1,3)",0); print(TT);
  TT=base_replace(S,[[x,2],[y,3]]); /* target 2 */
  print("Value at (2,3)",0); print(TT);
  A=rk2y([[y1,y1/x0,0],[y2,0,y2/x1]],
          [x0,x1],[y0,y1,y2],
          [1,0.1], II , [2,3],0.1);
}

def rk2y(F,X,Y,Xs,Ys,Xt,H) {
  X0=X[0]; X1=X[1];
  S0=Xs[0]; S1=Xs[1];
  T0=Xt[0]; T1=Xt[1];
  
  /* goto y direction */
  FF = base_replace(F[1],[[X0,S0]]);
  A=tk_rk.taka_runge_kutta_4a(FF,X1,Y,S1,Ys,T1,H);
  T=A[0];
/*  print("Mid point ",0); print(T); */

  /* goto  x direction */
  Ys = cdr(T);
  FF = base_replace(F[0],[[X1,T1]]);
  A=tk_rk.taka_runge_kutta_4a(FF,X0,Y,S0,Ys,T0,H);
  T=A[0];
  return cons([T[0],T1],cdr(T));
}


def test3() {
  S=[x^2-y^2, -2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,1],[y,3]]); /* target 1 */
  print("Value at (1,3)",0); print(TT);
  TT=base_replace(S,[[x,2],[y,3]]); /* target 2 */
  print("Value at (2,3)",0); print(TT);
  A=rk2y([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [2,3],0.1);
}


/*
  A=[[y1,y1/x0,0],[y2,0,y2/x1]] Pfaffian for x0 and x1.
  T=[x0,x1];  independent variable
  F=[y0,y1,y1]; dependent variable
  Ts=[s0,s1];  starting point
  Fs=          starting value of Y
  Tt=[t0,t1];  target point.
  H step
  Ts < Tt  or   Ts > Tt

  options  mesh  H/mesh is used for the steps in rk2x and rk2y. The data for these finer steps are discarded.
                Increase H and Mesh to get smaller list size as the return of rk2().
           eps   eps is used for workaround of the case that number_floor(2) < 2 holds.
*/
def rk2(A,T,F,Ts,Fs,Tt,H) {
  Mesh = getopt(mesh);
  if (type(Mesh)<0) Mesh=5;
  Eps = getopt(eps);
  if (type(Eps)<0) Eps=0.0001;
  X=T[0]; Y=T[1];
  Sx=Ts[0]; Sy=Ts[1];
  Tx=Tt[0]; Ty=Tt[1];
  Xsig = 1; Ysig = 1;
  if (Ty-Sy < 0) Ysig = -1;
  if (Tx-Sx < 0) Xsig = -1;
  Dy=Ysig*(Ty-Sy); 
  Dx=Xsig*(Tx-Sx); 
  Ans=[cons(Ts,Fs)];
  if ((Dx == 0) && (Dy==0)) return Ans;
  if (Dy < Dx) {
    Hx = H;
    Hy = H*Dy/Dx;
    Ux = Sx+Xsig*Hx; Uy=Sy+Ysig*Hy;
    while ( (Ux-Tx)*Xsig < Eps) {
      NN=rk2x(A,T,F,Ts,Fs,[Ux,Uy],H/Mesh);
      Ans=cons(NN,Ans);
      Ts=car(NN);
      Fs=cdr(NN);
      Ux=Ux+Xsig*Hx; Uy=Uy+Ysig*Hy;
    }
    NN=rk2x(A,T,F,Ts,Fs,Tt,H);
    return cons(NN,Ans);
  }else {
    Hy = H;
    Hx = H*Dx/Dy;
    Ux = Sx+Xsig*Hx; Uy=Sy+Ysig*Hy;
    while ( (Uy-Ty)*Ysig < Eps) {
      print(Uy);
      NN=rk2y(A,T,F,Ts,Fs,[Ux,Uy],H/Mesh);
      Ans=cons(NN,Ans);
      Ts=car(NN);
      Fs=cdr(NN);
      Ux=Ux+Xsig*Hx; Uy=Uy+Ysig*Hy;
    }
    NN=rk2y(A,T,F,Ts,Fs,Tt,H);
    return cons(NN,Ans);
  }
}

def test5() {
  S=[x^2-y^2, -2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,5],[y,3]]); /* target 2 */
  print("Value at (5,3)",0); print(TT);
  A=rk2([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [5,3],0.1);
}

def test5rr() {
  S=[x^2-y^2, -2*y];
  II=base_replace(S,[[x,3],[y,2.5]]); /* Initial */
  print("Value at (3,2.5)",0); print(II);
  TT=base_replace(S,[[x,1],[y,2]]); /* target 2 */
  print("Value at (1,2)",0); print(TT);
  A=rk2([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [3,2.5], II , [1,2],0.1);
}

def test5b() {
  S=[x^2-y^2, -2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,2],[y,3]]); /* target 2 */
  print("Value at (2,3)",0); print(TT);
  A=rk2([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [2,3],0.1);
}
def test5c() {
  S=[x^2-y^2, -2*y];
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  TT=base_replace(S,[[x,1],[y,3]]); /* target 2 */
  print("Value at (1,3)",0); print(TT);
  A=rk2([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [1,3],0.1);
}

def rk2all(A,T,F,Ts,Fs,Tt,H) {
  X=T[0]; Y=T[1];
  Sx=Ts[0]; Sy=Ts[1];
  Tx=Tt[0]; Ty=Tt[1];
  AnsX=rk2(A,T,F,Ts,Fs,[Tx,Sy],H);
  AnsY=rk2(A,T,F,Ts,Fs,[Sx,Ty],H);
  AnsY=reverse(AnsY);
  Ans =reverse(AnsX);
  AnsY = cdr(AnsY);
  while (AnsY != []) {
    AA = car(AnsY);
    /* solve to x direction */
    TT = rk2(A,T,F,AA[0],cdr(AA),[Tx,AA[0][1]],H);
    AnsY = cdr(AnsY);
    Ans = append(Ans,reverse(TT));
  }
  return Ans;
}

/* test6(4); --> 0.5 --> bad in test6c() */
/* test6(10);   --> good except x=1 */
def test6(N) {
  S=[x^2-y^2, -2*y];
  H=deval(2/N);
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  A=rk2all([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [3,3],H);
  /* setGtable0(A,N+1,number_floor(deval(N*3/2))+2,H,H); */
/*  return A; */
  mtg.setGtable(A); 
  return A;
}

/* tk_pf2.test6n(10) for testing... */
def test6n(N) {
  S=[x^2-y^2, -2*y];
  H=deval(2/N);
  II=base_replace(S,[[x,1],[y,0.1]]); /* Initial */
  print("Value at (1,0.1)",0); print(II);
  A=rk2all([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [1,0.1], II , [3,3],H);
  mtg.genGtable(A); 
  Dom=[[1,3],[0.1,3]];
  mtg.plot3d(quote(mtg.f_gtable(x,y)) | domain=Dom);
  sleep(2000);
  A2=test6(N);
  mtg.plot3d(quote(mtg.f_gtable(x,y)) | domain=Dom);
  return [A,A2];
}

def test6rr(N) {
  S=[x^2-y^2, -2*y];
  H=deval(2/N);
  II=base_replace(S,[[x,3],[y,3]]); /* Initial */
  print("Value at (3,3)",0); print(II);
  A=rk2all([[(2*y0-x1*y1)/x0,0],[y1,y1/x1]],
          [x0,x1],[y0,y1],
          [3,3], II , [1,1],H);
  /* setGtable0(A,N+1,number_floor(deval(N*3/2))+2,H,H); */
/*  return A; */
  mtg.setGtable(A); 
  return A;
}

endmodule;


end$