share/contrib/Grobner/grobner.demo

   1 /*  -*-  Mode: Maxima -*- */
   2
   3 /*
   4 **
   5 ** $Id: grobner.demo,v 1.2 2003-05-03 11:40:00 starseeker Exp $
   6 ** Copyright (C) 1999, 2002 Marek Rychlik <rychlik@u.arizona.edu>
   7 **
   8 ** This program is free software; you can redistribute it and/or modify
   9 ** it under the terms of the GNU General Public License as published by
  10 ** the Free Software Foundation; either version 2 of the License, or
  11 ** (at your option) any later version.
  12 **
  13 ** This program is distributed in the hope that it will be useful,
  14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 ** GNU General Public License for more details.
  17 **
  18 ** You should have received a copy of the GNU General Public License
  19 ** along with this program; if not, write to the Free Software
  20 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  21 **
  22 */
  23 showtime:true;
  24
  25 /* POLY_MONOMIAL_ORDER switch represents the monomial order that will globally be in effect
  26   for the succeeding operations. */
  27
  28 poly_monomial_order:'lex;
  29
  30 /* POLY_EXPAND parses polynomials to internal form and back. It can be used to test
  31   whether an expression correctly parses to the internal representation.
  32   The following examples illustrate that indexed and transcendental function variables
  33   are allowed. */
  34
  35 poly_expand(x,[x,y]);
  36 poly_expand(x+y,[x,y]);
  37 poly_expand(x-y,[x,y]);
  38 poly_expand((x-y)*(x+y),[x,y]);
  39 poly_expand((x+y)^2,[x,y]);
  40 poly_expand((x+y)^5,[x,y]);
  41 poly_expand(x/y-1,[x]);
  42 poly_expand(x^2/sqrt(y)-x*exp(y)-1,[x]);
  43 poly_expand(sin(x)-sin(x)^2-1,[sin(x)]);
  44 poly_expand((x[2]/sin(y[3])-1)^5,[x[2]]),poly_return_term_list:true;
  45
  46 /* POLY_ADD, POLY_SUBTRACT, POLY_MULTIPLY and POLY_EXPT are the arithmetical operations on polynomials.
  47    These are performed using the internal representation, but the results are converted back to the
  48    Maxima general form */
  49
  50 poly_add(x^2*y+z,x-z,[x,y,z]);
  51 poly_subtract(x^2*y+z,x-z,[x,y,z]);
  52 poly_multiply(x^2*y+z,x-z,[x,y,z]) - (x^2*y+z)*(x-z), expand;
  53 poly_expt(x-y, 3, [x,y]) - (x-y)^3, expand;
  54
  55 /* POLY_CONTENT extracts the GCD of its coefficients */
  56 poly_content(21*x+35*y,[x,y]);
  57
  58 /* POLY_PRIMITIVE_PART divides a polynomial by the GCD of its coefficients */
  59 poly_primitive_part(21*x+35*y,[x,y]);
  60
  61 /* POLY_S_POLYNOMIAL computest the syzygy polynomial (S-polynomial) of two polynomials */
  62 poly_s_polynomial(x+y,x-y,[x,y]);
  63
  64
  65 /* POLY_NORMAL_FORM finds the normal form of a polynomial with respect to a set of polynomials. */
  66 poly_normal_form(x^2+y^2,[x-y,x+y],[x,y]);
  67 poly_pseudo_divide(2*x^2+3*y^2,[7*x-y^2,11*x+y],[x,y]);
  68 poly_exact_divide((x+y)^2,x+y,[x,y]);
  69
  70 /* POLY_BUCHBERGER performs the Buchberger algorithm on a list of polynomials and returns
  71    the resulting Grobner basis */
  72 poly_buchberger([x^2-y*x,x^2+y+x*y^2],[x,y]);
  73
  74 /* POLY_REDUCTION reduces a set of polynomials, so that
  75   each polynomial is fully reduced with respect to the other polynomials */
  76
  77 poly_reduction([x^2-x*y,x*y^2+y+x^2,x*y^2+x*y+y,x*y-y^2,y^3+y^2+y],[x,y]);
  78
  79 /* POLY_MINIMIZATION selects a subset of a set of polynomials, so that no leading monomial is divisible by
  80   another leading monomial */
  81
  82 poly_minimization([x^2-x*y,x*y^2+y+x^2,x*y^2+x*y+y,x*y-y^2,y^3+y^2+y],[x,y]);
  83
  84 /* POLY_REDUCED_GROBNER returns a reduced Grobner basis */
  85 poly_reduced_grobner([x^2-y*x,x^2+y+x*y^2],[x,y]);
  86
  87 /* POLY_NORMALIZE divides a polynomial by its leading coefficient */
  88 poly_normalize(2*x+y,[x,y]);
  89
  90 /* POLY_NORMALIZE_LIST applies POLY_NORMALIZE to each polynomial in the list */
  91
  92 poly_normalize_list([2*x+y,3*x^2+7],[x,y]);
  93
  94 /* POLY_DEPENDS_P tests whether a polynomial depends on a variable */
  95
  96 poly_depends_p(x^2+y,x,[x,y,z]);
  97 poly_depends_p(x^2+y,z,[x,y,z]);
  98
  99
 100 /* POLY_ELIMINATION_IDEAL returns the grobner basis of the K-th elimination ideal of an
 101    ideal specified as a list of generating polynomials (not necessarily Grobner basis */
 102
 103 poly_elimination_ideal([x+y,x-y],0,[x,y]);
 104 poly_elimination_ideal([x+y,x-y],1,[x,y]);
 105 poly_elimination_ideal([x+y,x-y],2,[x,y]);
 106
 107 /* POLY_IDEAL_INTERSECTION returns the intersection of two ideals */
 108 poly_ideal_intersection([x^2+y,x^2-y],[x,y^2],[x,y]);
 109
 110 /* POLY_LCM and POLY_GCD are the Grobner versions of LCM and GCD */
 111
 112 poly_lcm(x*y^2-x,x^2*y+x,[x,y]);
 113 poly_gcd(x*y^2-x,x^2*y+x,[x,y]);
 114
 115 /* POLY_GROBNER_MEMBER tests whether a polynomial belongs to an ideal generated by a list of polynomials,
 116    which is assumed to be a Grobner basis. Equivalent to NORMAL_FORM being 0. */
 117
 118 poly_grobner_member(x+y,[x,y],[x,y]);
 119
 120 /* POLY_GROBNER_EQUAL tests whether two Grobner bases generate the same ideal.
 121    This is equivalent to checking that every polynomial of the first basis reduces to 0
 122    modulo the second basis and vice versa. Note that in the example below the
 123    first list is not a Grobner basis, and thus the result is FALSE. */
 124
 125 poly_grobner_equal([x+y,x-y],[x,y],[x,y]);
 126
 127 /* POLY_GROBNER_SUBSETP tests whether an ideal generated by the first list of polynomials
 128    is contained in the ideal generated by the second list. For this test to always succeed,
 129    the second list must be a Grobner basis */
 130
 131 poly_grobner_subsetp([x+y,x-y],[x,y],[x,y]);
 132
 133 /* POLY_POLYSATURATION_EXTENSION implements the famous Rabinowitz trick. */
 134 poly_polysaturation_extension([x,y],[x^2,y^3],[x,y],[u,v]);
 135
 136 poly_saturation_extension([x,y],[x^2,y^3],[x,y],[u,v]);
 137
 138 /* POLY_IDEAL_POLYSATURATION1 for a given ideal I and polynomials f, g, ..., finds
 139    the colon ideal I : f^inf : g^inf : ... */
 140 poly_ideal_polysaturation1([x,y],[x^2,y^3],[x,y]);
 141
 142 /* POLY_IDEAL_SATURATION for given ideals I and J computes the ideal I : J^inf. */
 143 poly_ideal_saturation([x,y],[x^2,y^3],[x,y]);
 144
 145 /* POLY_IDEAL_POLYSATURATION for a given ideal I and a sequence of ideals J1, J2, J3, ...,
 146   finds the ideal I : J1^inf : J2^inf : J3^inf : ... */
 147 poly_ideal_polysaturation([x,y],[[x^2],[y^3]],[x,y]);
 148 poly_ideal_polysaturation([x^4-y^4], [[x-y],[x^2+y^2, x+y]],[x,y]);
 149
 150 /* POLY_COLON_IDEAL finds the reduced Grobner basis of the colon ideal I:J, i.e. the set of polynomials h
 151    such that there is a polynomial F in J for which H*F is in I */
 152
 153 poly_colon_ideal([x^2*y],[y],[x,y]);
 154
 155 /* POLY_BUCHBERGER_CRITERION verifies whether a given set of polynomials is a Grobner basis with respect
 156    to the current term order */
 157 poly_buchberger_criterion([x,y],[x,y]);
 158 poly_buchberger_criterion([x-y,x+y],[x,y]);
 159
 160 /* Grobner basis associated with Enneper minimal surface */
 161 poly_grobner([x-3*u-3*u*v^2+u^3,y-3*v-3*u^2*v+v^3,z-3*u^2+3*v^2],[u,v,x,y,z]);
 162 poly_reduced_grobner([x-3*u-3*u*v^2+u^3,y-3*v-3*u^2*v+v^3,z-3*u^2+3*v^2],[u,v,x,y,z]);
 163
 164 /* Cyclic roots of degree 5 */
 165 poly_reduced_grobner([x+y+z+u+v,x*y+y*z+z*u+u*v+v*x,x*y*z+y*z*u+z*u*v+u*v*x+v*x*y,x*y*z*u+y*z*u*v+z*u*v*x+u*v*x*y+v*x*y*z,x*y*z*u*v-1],[u,v,x,y,z]);
 166
 167 /* The next example demonstrates the use of the switch
 168   POLY_RETURN_TERM_LIST, which, if set to TRUE, makes the results to
 169   appear as lists of terms listed in the current monomial order rather
 170   than a general form expression */
 171
 172 block([orders:[lex,grlex,grevlex,invlex]],
 173 for i:1 thru length(orders) do (
 174   print(ev([orders[i], poly_expand((x^2+x+y)^3,[x,y])], poly_monomial_order=orders[i]))
 175   )
 176 ), poly_return_term_list=true;
 177
 178 /* Grobner bases can be computed over coefficient ring of maxima general expressions */
 179 poly_grobner([x*y-1,x+y],[x]);
 180
 181 /* A tough example learned from Cox */
 182 poly_grobner([x^5+y^4+z^3-1,x^3+y^3+z^2-1], [x,y,z]);
 183
 184 /* An even tougher example of Cox */
 185 poly_grobner([x^5+y^4+z^3-1,x^3+y^3+z^2-1], [x,y,z]);
 186
 187 /* We can also perform Grobner basis calculations modulo prime */
 188 poly_grobner([x^5+y^4+z^3-1,x^3+y^3+z^2-1], [x,y,z]), modulus=3;
 189
 190 /* We can also explicitly ask for the Grobner basis to be calculated using only
 191   integer coefficients. An error will result if this assertion is not satisfied. */
 192 poly_grobner([x^5+y^4+z^3-1,x^3+y^3+z^2-1], [x,y,z]), poly_coefficient_ring='ring_of_integers;
 193
 194 /* The following several tests demonstrate the use of jet variables useful in processing differential equations */
 195
 196
 197 /* Clear some variables */
 198 kill(ode,t,x,y,u,v,r);
 199
 200 /* Set up dependencies */
 201 depends([x,y,u,v,r],t);
 202
 203 /* These are equations representing mathematical pendulum */
 204 ode:[x^2+y^2-c,'diff(x,t)-u,'diff(y,t)-v,'diff(u,t)+r*x,'diff(v,t)+r*y+1];
 205
 206 jet_vars(k):=apply(append,reverse(makelist(['diff(x,t,j),'diff(y,t,j),'diff(u,t,j),'diff(v,t,j),'diff(r,t,j)],j,0,k+1)));
 207
 208 /* Define k-fold prolongation */
 209 prolongate(k):=apply(append,makelist(diff(ode,t,j),j,0,k));
 210
 211 /* Define Grobner basis of k-fold prolongation */
 212 gb(k):=poly_reduced_grobner(prolongate(k),jet_vars(k));
 213
 214 /* Define the l-th projection of the k-th prolongation */
 215 projection(l, k):=poly_elimination_ideal(prolongate(k),5*l,jet_vars(k));
 216
 217 /* Compute some projections */
 218 projection(0, 0);
 219 projection(1, 1);