Return to taka_parter.rr CVS log

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

A system of algeraic equations with parameters
derived from a diffusion-reaction equation.
See comments of the file for usages, references, and problems.

/* 
  $OpenXM: OpenXM/src/asir-contrib/packages/src/taka_parter.rr,v 1.1 2006/10/04 23:18:30 takayama Exp $

  This pacakge requires asir-contrib functions.

  cf. this04/misc-2004/A1/reaction-diff/parter.rr 
*/

module taka_parter;
localf set_globals;
localf parter;
localf bif;
localf parter_shape;
localf parter_shape1;
localf demo1;
localf get_real_roots;
localf p3;
localf p3a;
localf p2;
localf demo1_2;
localf parter_shape1c;

/*&usage begin: taka_parter.parter(N)
 It returns the system of algebraic equations parter(N), which contains paramters d and a.
 example: taka_parter.parter(4);
end: */

/*&usage begin: taka_parter.demo1(N)
 It draws a graph of the parameter a and the real values of the unknown u_1. 
 example: taka_parter.set_globals();  taka_parter.demo1(4);
end: */

/*  Consider the reaction diffusion equation.
    u_t = d u_xx + a(1-u) u,
    u(t,0)=u(t,1) = 0,  (Question: consider also this equation with the Neumann condition)

   The parter(N) is a system of algebraic equations with parameters a and d
   and unknowns u_1, ..., u_N obtained by discretizing d u_xx + a(1-u)u = 0.

   We are interested in solutions satisfying 0< u_1, ..., u_N < 1.

  References
  1.  Parter, S.V., Mildly nonlinear elliptic partial differential equations and
  their numerical solutions, I, Numer. Math., 7 (1965), 338--366.  
  2. Mimura,  $BHyJ,J}Dx<0$H:9J,J}Dx<0(B $B$b;2>H(B. [$B?tCM2r@O$HHs@~7A8=>](B]
    p.60

  Theorem:
  a = xi_1,  xi_1= (4d/h^2)*sin^2(h*pi/2) is the bufurcation point.  
  We have a real solution 0 < u_i < 1 when a > xi_1.

  Problem 0: How far can we compute the shape basis of parter(N)?
             How far can we compute the primary ideal decomposition of parter(N)?
  Problem 1: Can we rediscover this theorem with a help of mathematical software?
  Problem 2: Get a series expansion of solutions in terms of a and d around the bifurcation point.
*/

/* set orders */
def set_globals() {
  extern Glib_math_coordinate;
  ord([u_1,a,d])$
  pari(allocatemem,10^6)$
  Glib_math_coordinate=1$
}


def parter(N) {
  H = 1/(N+1);
  Eq = [ ];
  V = [ ];
  for (I=N; I >= 1; I--) {
     U0 = util_v("u",[I-1]);
     U1 = util_v("u",[I]);
     U2 = util_v("u",[I+1]);
     if (I==1) U0 = 0;
     if (I==N) U2 = 0;
     Eq = cons( d*(U0 - 2*U1+U2) + H*a*(1-U1)*U1, Eq);
     V = cons(U1,V);
  }
  return [Eq,V];
}

/* Bifurcation value of a. */
def bif(N) {
   H = 1/(N+1);
   return (4*d/H^2)*sin(H*@pi/2)^2;
}

/* Return the shape base of Parter's system. */
def parter_shape(N) {
  Eqs = newvect(N+2);
  H = 1/(N+1);
  /* u_i = Eqs[i]  $B$,4X78<0(B. Eqs[i] $B$O(B u_1 $B$GI=8=$5$l$F$$$k(B. */
  /* $B$H$/$K(B u_{N+1} = 0 $B$J$N$G(B, Eqs[N+1] = 0 $B$,(B u_1 $B$NK~$?$9J}Dx<0(B. */
  Eqs[1] = util_v("u",[1]);
  for (I=1; I<=N; I++) {
     U0 = util_v("u",[I-1]);
     U1 = util_v("u",[I]);
     if (I == 1) U0=0;
     U2 = 2*U1-U0 - (H^2/d)*a*U1*(1-U1);
     if (I != 1) U2 = subst(U2,U0,Eqs[I-1]);
     U2 = subst(U2,U1,Eqs[I]);
     Eqs[I+1] = U2;
  }
  return Eqs;
}

/* Algebraic equation of u_1. */
def parter_shape1(N) {
  E = parter_shape(N);
  F = nm(E[N+1]);
  return F;
}


def demo1(N) {
  F = parter_shape1(N);
  F = subst(F,d,1);

  D = bif(N);
  D = subst(D,d,1);
  D = eval(D);

  glib_window(-D*2,-10,D*2,10);
  print("Bifurcation point a="+rtostr(D));
  for (A = -D*2; A < D*2; A += eval(D/5)) {
     R = pari(roots,subst(F,a,A));
     /* print([A,R]); */
     if (R == 0) {
       print(subst(F,a,A));     
     } else {
       R1 = get_real_roots(R);
       print([A,R1]);
       N = length(R1);
       for (I=0; I<N; I++) {
         if ((R1[I] < 10) && (R1[I] > -10)) { 
           glib_putpixel(A,R1[I]);
         }
       }
     }
  }
}

def get_real_roots(R) {
  N = length(R);
  A = [ ];
  for (I=0; I<N; I++) {
     if (number_imaginary_part(R[I]) == 0) {
         A = cons(R[I],A);
     }
  }
  return A;
}

/* parter_shape1(3) = (degree 3)*(degree 4) */
def p3() {
  F=a^3*u_1^3+(-2*a^3+64*d*a^2)*u_1^2+(a^3-80*d*a^2+1536*d^2*a)*u_1+16*d*a^2-1024*d^2*a+8192*d^3;
  F1 = diff(F,u_1);
  R = res(u_1,F,F1);
  G = a^4*u_1^4+(-2*a^4+64*d*a^3)*u_1^3+(a^4-80*d*a^3+1536*d^2*a^2)*u_1^2+(16*d*a^3-1024*d^2*a^2+16384*d^3*a)*u_1-4096*d^3*a+131072*d^4;
  G1 = diff(G,u_1);
  RR = res(u_1,G,G1);
  print("---9.37----");
  print(pari(roots,subst(R,d,1)));
  print(pari(roots,subst(RR,d,1)));
  return [R,RR];
}

def p3a() {
  F= parter_shape1(3);
  F1 = diff(F,u_1);
  R = res(u_1,F,F1);
  print(subst(subst(F,d,1),a,9.372583));
  print(pari(roots,subst(R,d,1)));
}
/*
  c_8 u_1^8+...+c_1 u_1 = 0
  We are interested in the asymptotic expansion of solutions around c_1=0. 
*/

def p2() {
  N = 2;
  R = parter_shape1(N);
  print(R);
/*  demo1(2); */
  D = coef(R,1);
  print(fctr(D));
  M = newmat(3,3,[[-2,1,0],[1,-2,1],[0,1,-2]]);
  M = M*d/((1/3)^2);
  L = newmat(3,3,[[ell,0,0],[0,ell,0],[0,0,ell]]);
  print(fctr(det(M-L))); 
  return R;
}

def demo1_2() {
  N = 2;
  F = parter_shape1(N);
  F = subst(F,d,1);

  D = bif(N);
  D = subst(D,d,1);
  D = eval(D);

  glib_window(-D*2,-10,D*2,10);
  print("Bifurcation point a="+rtostr(D));
  for (A = -D*2; A < D*2; ) {
     R = pari(roots,subst(F,a,A));
     /* print([A,R]); */
     if (R == 0) {
       print(subst(F,a,A));     
     } else {
       R1 = get_real_roots(R);
       print([A,R1]);
       N = length(R1);
       for (I=0; I<N; I++) {
         glib_putpixel(A,R1[I]);
       }
     }
     Skip=0.2;
     if ((A > -Skip) && (A < Skip)) {
        A = Skip;
     }
     if ((A > -1) && (A < 1)) {
        A += 0.01;
     }else { 
        A += 0.1;
     }
  }
}

/* 
   parter 5 (degree 1) or (degree 2)
   parter 6  d^54*(degree 3)*(degree 3)
   Conj: R has solution at a=(4d/h^2)*sin^2(h*pi/2), h = 1/(N+1).
*/
def parter_shape1c(N) {
  E = parter_shape(N);
  F = nm(E[N+1]);
  R=coef(F,1);
  print(fctr(R));
  R1=subst(R,d,1);
  print(pari(roots,R1));
  return(fctr(R));
}


endmodule;

end$