# source:git/Singular/LIB/findiff.lib@66d68c

spielwiese
Last change on this file since 66d68c was 66d68c, checked in by Hans Schoenemann <hannes@…>, 13 years ago
• Property mode set to `100644`
File size: 41.1 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Teaching";
4info="
5LIBRARY:  findiff.lib    procedures to compute finite difference schemes for
6                         linear differential equations
7AUTHOR:                  Christian Dingler
8
9OVERVIEW:
10@texinfo
12 library requires the program @code{qepcad} to be installed.
15 @end texinfo
16
17PROCEDURES:
18 visualize(f);            shows a scheme in index-notation
19 u(D[,#]);                gives some vector; depends on @derivatives
20 scheme([v1,..,vn]);      computes the finite difference scheme defined by v1,..,vn
21 laxfrT(Ut,U,space);      Lax-Friedrich-approximation for the time-direction
22 laxfrX(Ux,U,space);      Lax-Friedrich-approximation for the space-direction
23 forward(U1,U2,VAR);      forward-approximation
24 backward(U1,U2,VAR);     backward-approximation
25 central1st(U1,U2,VAR);   central-approximation of first order
26 central2nd(U1,U2,VAR);   central-approximation of second order
27 trapezoid(U1,U2,VAR);    trapezoid-approximation
28 midpoint(U1,U2,VAR);     midpoint-approximation
29 pyramid(U1,U2,VAR);      pyramid-approximation
30 setinitials(variable,der[,#]);    constructs and sets the basering for further computations
31 errormap(f);        performs the Fouriertransformation of a poly
32 matrixsystem(M,A);    gives the scheme of a pde-system as one matrix
33 timestep(M);        gives the several timelevels of a scheme derived from a pde-system
34 fouriersystem(M,A);  performs the Fouriertransformation of a matrix scheme
35 PartitionVar(f,n);      partitions a poly into the var(n)-part and the rest
36 ComplexValue(f);        computes the complex value of f, var(1) being the imaginary unit
37 VarToPar(f);            substitute var(i) by par(i)
38 ParToVar(f);            substitute par(i) by var(i)
41";
42
43LIB "ring.lib";
44LIB "general.lib";
45LIB "standard.lib";
46LIB "linalg.lib";
47LIB "matrix.lib";
48LIB "poly.lib";
49LIB "teachstd.lib";
50LIB "qhmoduli.lib";
51///////////////////////////////////////////////////////////////////////
52static proc getit(module M)
53{
54  int nderiv=pos(U,@derivatives);
55  def M2=groebner(M);
56  module N1=gen(nderiv);
57  def N2=intersect(M2,N1);
58  def S=N2[1][nderiv];
59  return(S);
60}
61///////////////////////////////////////////////////////////////////////
62proc visualize(poly f)
63"USAGE:   visualize(f); f of type poly.
64RETURN:  type string; translates the polynomial form of a finite difference scheme into an indexed one as often seen in literature
65EXAMPLE: example visualize; shows an example
66"
67{
68  def n=size(f);
69  string str;
70  intvec v;
71  if (n>0)
72  {
73    int i;
74    int j;
75    for(i=1;i<=n;i++)
76    {
78      for(j=1;j<=size(@variables);j++)
79      {
80        v[j]=w[j+1];
81      }
82      if(i==1)
83      {
85      }
86      else
87      {
89      }
90      kill w;
92    }
93  }
94  return(str);
95}
96example
97{
98"EXAMPLE:";echo=2;
99 list D="Ux","Ut","U";
100 list P="a";
101 list V="t","x";
102 setinitials(V,D,P);
103 scheme(u(Ut)+a*u(Ux),trapezoid(Ux,U,x),backward(Ut,U,t));
104 visualize(_);
105}
106
107///////////////////////////////////
108static proc imageideal()
109{
110
111  def n=size(@variables)-1;
112  ideal IDEAL=var(1),var(2);
113  int j;
114  for(j=1;j<=n;j++)
115  {
116    ideal II=var(2+j+n)+var(1)*var(2+2*n+j);
117    IDEAL=IDEAL+II;
118    kill II;
119  }
120  return(IDEAL);
121}
122/////////////////////////////////////
123proc u(D,list #)
124"USAGE: u(D[,#]); D a string that occurs in the list of @derivatives, # an optional list of integers.
125RETURN: type vector; gives the vector, that corresponds with gen(n)*m, where m is the monomial defined by #
126EXAMPLE: example u; shows an example
127"
128{
129  def n=size(#);
130  def nv=nvars(basering)-1;
131  int nn;
132  if(nv<=n)
133  {
134    nn=nv;
135  }
136  else
137  {
138    nn=n;
139  }
140  int index=pos(D,@derivatives);
141  poly g=1;
142  if(nn>=1)
143  {
144    int j;
145    for(j=1;j<=nn;j++)
146    {
147      int nnn;
148      nnn=int(#[j]);
149      g=var(1+j)**nnn*g;
150      kill nnn;
151    }
152    return(gen(index)*g);
153  }
154  else
155  {
156    return(gen(index)*g);
157  }
158}
159example
160{
161  "EXAMPLE:";echo=2;
162  list D="Ux","Uy","Ut","U";
163  list P="a","b";
164  list V="t","x","y";
165  setinitials(V,D,P);
166  u(Ux);
167  u(Ux,2,3,7);
168  u(Uy)+u(Ut)-u(Ux);
169  u(U)*234-dx*dt*dy*3*u(Uy);
170}
171
172/////////////////////////////////////////////////////////
173static proc pos(string D,list L)
174{
175  int j;
176  int index=-1;
177  def n=size(L);
178  for(j=1;j<=n;j++)
179  {
180    if(D==L[j])
181    {
182      index=j;
183    }
184  }
185  return(index);
186}
187///////////////////////////////////
188static proc re(list L)
189{
190  def n=size(L);
191  int j;
192  for(j=1;j<=n;j++)
193  {
194    string s="string "+print(L[j],"%s")+"="+"nameof("+print(L[j],"%s")+")"+";";
195    execute(s);
196    kill s;
197    exportto(Top,`L[j]`);
198  }
199}
200///////////////////////////////////////////////
201proc scheme(list #)
202"USAGE:  scheme([v1,..,vn]); v1,..,vn of type vector
203RETURN:  poly
204PURPOSE: performs substitutions by the means of Groebner basis computation
205of the submodule, generated by the input vectors, then intersects the
206intermediate result with the suitable component in order to get a finite
207difference scheme
208NOTE:  works only for a single PDE, for the case of a system use @code{matrixsystem}
209EXAMPLE: example scheme; shows an example
210"
211{
212  def N=size(#);
213  if(N==0)
214  {
215    if(defined(M)==1)
216    {
217      kill M;
218      module M;
219    }
220    else
221    {
222      module M;
223    }
224  }
225  else
226  {
227    int j;
228    if(defined(M)==1)
229    {
230      kill M;
231      module M;
232      for(j=1;j<=N;j++)
233      {
234        M=M+#[j];
235      }
236    }
237    else
238    {
239      module M;
240      for(j=1;j<=N;j++)
241      {
242        M=M+#[j];
243      }
244    }
245  }
246  def S=getit(M);
247  matrix mat[1][1]=S;
248  list l=timestep(mat);
249  poly f=l[2][1,1];
250  return(f);
251}
252example
253{
254  "EXAMPLE:";echo=2;
255  list D="Ux","Ut","U";
256  list P="a";
257  list V="t","x";
258  setinitials(V,D,P);
259  def s1=scheme(u(Ut)+a*u(Ux),backward(Ux,U,x),forward(Ut,U,t));
260  s1;
261}
262////////////////////////
263static proc diffpar(poly ff)
264{
265  def gg=print(ff,"%s");
266  def str="d"+gg;
267  return(`str`);
268}
269////////////////////////
270proc laxfrT(string Ut, string U, poly space)
271"USAGE:   laxfrT(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
272RETURN:   type vector; gives a predefined approximation of the Lax-Friedrich-approximation for the derivation in the timevariable as often used in literature;
273NOTE:     see also laxfrX, setinitials, scheme; Warning: laxfrT is not to be interchanged with laxfrX
274EXAMPLE:  example laxfrT; shows an example
275"
276{
277  poly time=var(2);
278  poly dtime=diffpar(time);
279  poly dspace=diffpar(space);
280  def v=dtime*space*u(Ut)-time*space*u(U)+1/2*(space**2*u(U)+u(U));
281  return(v);
282}
283example
284{
285  "EXAMPLE:";echo=2;
286  list D="Ux","Ut","U";
287  list P="a";
288  list V="t","x";
289  setinitials(V,D,P);
290  laxfrT(Ux,U,x);
291}
292////////////////////////
293proc laxfrX(string Ux, string U, poly space)
294"USAGE: laxfrX(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
295RETURN: type vector; gives a predefined approximation of the Lax-Friedrich-approximation for the derivation in one of the spatial variables as often used in literature;
296NOTE:   see also laxfrT, setinitials, scheme; Warning: laxfrX is not to be interchanged with laxfrT
297EXAMPLE: example laxfrX; shows an example
298"
299{
300  poly dspace = diffpar(space);
301  def v=2*dspace*space*u(Ux)-(space**2-1)*u(U);
302  return(v);
303}
304example
305{
306  "EXAMPLE:";echo=2;
307  list D="Ux","Ut","U";
308  list P="a";
309  list V="t","x";
310  setinitials(V,D,P);
311  laxfrX(Ux,U,x);
312}
313////////////////////////
314proc forward(string U1,string U2,poly VAR)
315"USAGE:   forward(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
316RETURN:   type vector; gives a predefined approximation of the forward approximation as often used in literature;
318EXAMPLE:  example forward; shows an example
319"
320{
321  if(pos(U1,@derivatives)<pos(U2,@derivatives))
322    {
323      def V1=U1;
324      def V2=U2;
325    }
326  else
327    {
328      def V1=U2;
329      def V2=U1;
330    }
331  def v=diffpar(VAR)*u(V1)+u(V2)-VAR*u(V2);
332  return(v);
333}
334example
335{
336  "EXAMPLE:";echo=2;
337  list D="Ut","Ux","Uy","U";
338  list V="t","x","y";
339  list P="a","b";
340  setinitials(V,D,P);
341  forward(Ux,U,x);
342  forward(Uy,U,y);
343  forward(Ut,U,t);
344}
345///////////////////////
346proc backward(string U1,string U2,poly VAR)
347"USAGE:   backward(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
348RETURN:   type vector; gives a predefined approximation of the backward approximation as often used in literature;
350EXAMPLE:  example backward; shows an example
351"
352{
353  if(pos(U1,@derivatives)<pos(U2,@derivatives))
354    {
355      def V1=U1;
356      def V2=U2;
357    }
358  else
359    {
360      def V1=U2;
361      def V2=U1;
362    }
363  def v=diffpar(VAR)*VAR*u(V1)+u(V2)-VAR*u(V2);
364  return(v);
365}
366example
367{
368  "EXAMPLE:";echo=2;
369  list D="Ut","Ux","Uy","U";
370  list V="t","x","y";
371  list P="a","b";
372  setinitials(V,D,P);
373  backward(Ux,U,x);
374  backward(Uy,U,y);
375  backward(Ut,U,t);
376}
377/////////////////////////////
378proc central1st(string U1,string U2,poly VAR)
379"USAGE:   central1st(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
380RETURN:   type vector; gives a predefined approximation of the first-order-central-approximation as often used in literature;
382EXAMPLE:  example central1st; shows an example
383"
384{
385  if(pos(U1,@derivatives)<pos(U2,@derivatives))
386    {
387      def V1=U1;
388      def V2=U2;
389    }
390  else
391    {
392      def V1=U2;
393      def V2=U1;
394    }
395  def v=2*diffpar(VAR)*VAR*u(V1)+u(V2)-VAR**2*u(V2);
396  return(v);
397}
398example
399{
400  "EXAMPLE:";echo=2;
401  list D="Ut","Ux","Uy","U";
402  list V="t","x","y";
403  list P="a","b";
404  setinitials(V,D,P);
405  central1st(Ux,U,x);
406  central1st(Uy,U,y);
407}
408////////////////////////////////
409proc central2nd(string U1,string U2,poly VAR)
410"USAGE:   central2nd(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
411RETURN:   type vector; gives a predefined approximation of the second-order-central-approximation as often used in literature;
413EXAMPLE:  example central2nd; shows an example
414"
415{
416  if(pos(U1,@derivatives)<pos(U2,@derivatives))
417    {
418      def V1=U1;
419      def V2=U2;
420    }
421  else
422    {
423      def V1=U2;
424      def V2=U1;
425    }
426  def v=diffpar(VAR)**2*VAR*u(V1)-(VAR**2*u(V2)-2*VAR*u(V2)+u(V2));
427  return(v);
428}
429example
430{
431  "EXAMPLE:";echo=2;
432  list D="Uxx","Ux","Utt","Ut","U";
433  list P="lambda";
434  list V="t","x";
435  setinitials(V,D,P);
436  central2nd(Uxx,U,x);
437  central2nd(Utt,U,t);
438}
439/////////////////////////////////
440proc trapezoid(string U1,string U2,poly VAR)
441"USAGE:   trapezoid(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
442RETURN:   type vector; gives a predefined approximation of the trapezoid-approximation as often used in literature;
444EXAMPLE:  example trapezoid; shows an example
445"
446{
447  if(pos(U1,@derivatives)<pos(U2,@derivatives))
448    {
449      def V1=U1;
450      def V2=U2;
451    }
452  else
453    {
454      def V1=U2;
455      def V2=U1;
456    }
457  def v=1/2*diffpar(VAR)*(VAR+1)*u(V1)+(1-VAR)*u(V2);
458  return(v);
459}
460example
461{
462  "EXAMPLE:";echo=2;
463  list D="Uxx","Ux","Utt","Ut","U";
464  list P="lambda";
465  list V="t","x";
466  setinitials(V,D,P);
467  trapezoid(Uxx,Ux,x);
468  trapezoid(Ux,U,x);
469}
470///////////////////////////////////
471proc midpoint(string U1,string U2,poly VAR)
472"USAGE:   midpoint(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
473RETURN:   type vector; gives a predefined approximation of the midpoint-approximation as often used in literature;
475EXAMPLE:  example midpoint; shows an example
476"
477{
478  if(pos(U1,@derivatives)<pos(U2,@derivatives))
479    {
480      def V1=U1;
481      def V2=U2;
482    }
483  else
484    {
485      def V1=U2;
486      def V2=U1;
487    }
488  def v=2*diffpar(VAR)*VAR*u(V1)+(1-VAR**2)*u(V2);
489  return(v);
490}
491example
492{
493  "EXAMPLE:";echo=2;
494  list D="Uxx","Ux","Utt","Ut","U";
495  list P="lambda";
496  list V="t","x";
497  setinitials(V,D,P);
498  midpoint(Ux,U,x);
499}
500//////////////////////////////////////
501proc pyramid(string U1,string U2,poly VAR)
502"USAGE:   pyramid(U1,U2,var); U1, U2 are the names of occuring derivatives, var is a variable in the basering;
503RETURN:   type vector; gives a predefined approximation of the pyramid-approximation as often used in literature;
505EXAMPLE:  example pyramid; shows an example
506"
507{
508  if(pos(U1,@derivatives)<pos(U2,@derivatives))
509    {
510      def V1=U1;
511      def V2=U2;
512    }
513  else
514    {
515      def V1=U2;
516      def V2=U1;
517    }
518  def v=1/3*diffpar(VAR)*(VAR**2+VAR+1)*u(V1)+(VAR-VAR**3)*u(V2);
519  return(v);
520}
521example
522{
523  "EXAMPLE:";echo=2;
524  list D="Uxx","Ux","Utt","Ut","U";
525  list P="lambda";
526  list V="t","x";
527  setinitials(V,D,P);
528  pyramid(Ux,U,x);
529}
530//////////////////////////////////////////////
531proc setinitials(list variable, list der,list #)
532"USAGE:  setinitials(V,D[,P]); V,D,P are lists with strings as elements
533RETURN:  no return value: sets the dependence order of the occuring derivatives,
534constructs the suitable ring to compute in containing user chosen parameters, sets new basering
535NOTE:    P is optional, used to introduce some additional parameters into the ring. The Sine and
536Cosine values needed for the fourier transformation are symbolically introduced under the names
537string(c)+nameof(variable), i.e. if x is any spatial variable then cx:=cosine(dx*ksi), when
538regarding the fourier transform after ksi (for sine respectively). Artificial parameters I,T,Px,Py
539are introduced for the later eigenvalue analysis. Variables can be transformed into parameters
540of similar name
541EXAMPLE: example setinitials; shows an example
542"
543{
544  def LV=variable;
545  def @variables=variable;
546  def @derivatives=der;
547  exportto(Top,@variables);
548  exportto(Top,@derivatives);
549  re(der);
550  int j;
551  string dvar="d"+print(LV[1],"%s");
552  dvar=dvar+","+"d"+print(LV[2],"%s");
553  string pvar="P"+print(LV[2],"%s");
554  string COS="C"+print(LV[2],"%s");
555  string SIN="S"+print(LV[2],"%s");
556  for(j=3;j<=size(LV);j++)
557    {
558      dvar=dvar+","+"d"+print(LV[j],"%s");
559      pvar=pvar+","+"P"+print(LV[j],"%s");
560      COS=COS+","+"C"+print(LV[j],"%s");
561      SIN=SIN+","+"S"+print(LV[j],"%s");
562    }
563  string scf="(0,"+"I,"+"T,"+pvar+","+COS+","+SIN+","+print(#,"%s")+","+dvar+")";  //coefficient_field
564  string svars="(i";
565  kill j;
566  int j;
567  for(j=1;j<=size(LV);j++)
568    {
569      svars=svars+","+print(LV[j],"%s");
570    }
571  string cosine;
572  string sine;
573  kill j;
574  int j;
575  cosine="c"+print(LV[2],"%s");
576  sine="s"+print(LV[2],"%s");
577  for(j=3;j<=size(LV);j++)
578    {
579      cosine=cosine+","+"c"+print(LV[j],"%s");
580      sine=sine+","+"s"+print(LV[j],"%s");
581    }
582  kill j;
583  string strvar=cosine+","+sine+")";
584  svars=svars+","+strvar;            ////variables
585  string sord="(c,lp)";         ////ordering
586  string sring="ring Q="+scf+","+svars+","+sord+";";
587  execute(sring);
588  ideal Id=i**2+1;
589  int j;
590  for(j=1;j<=size(LV)-1;j++)
591    {
592      ideal II=var(2+j+size(LV)-1)**2+var(2+j+2*(size(LV)-1))**2-1;
593      Id=Id+II;
594      kill II;
595    }
596  if(defined(basering)==1)
597    {
598      kill basering;
599    }
600  qring R=std(Id);
601  setring R;
602  // comment by VL: it's better to return this ring! causes many changes
603  // across the library
604  export(R);
605  exportto(Top,basering);
606}
607example
608{
609  "EXAMPLE:"; echo = 2;
610  list D="Ut","Ux","Uy","U";
611  list V="t","x","y";
612  list P="alpha","beta","gamma";
613  setinitials(V,D,P);////does not show the ring, since there is no output
614  basering;///does show the ring
615}
616////////////////////////////
617proc errormap(poly f)
618"USAGE:    errormap(f); f of type poly
619RETURN:   type poly; performs the fouriertransformation of a single polynomial
620EXAMPLE:  example errormap; shows an example
621"
622{
623  ideal Id=imageideal();
624  map phi=R,Id;
625  def g=phi(f);
626  g=reduce(g,std(0));
627  return(g);
628}
629example
630{
631  "EXAMPLE";echo=2;
632  list D="Ux","Ut","U";
633  list P="a";
634  list V="t","x";
635  setinitials(V,D,P);
636  scheme(u(Ut)+a*u(Ux),central1st(Ux,U,x),backward(Ut,U,t));
637  errormap(_);
638}
639/////////////////////////////////////
640static proc stepmatrix(int n, poly f)
641{
642  int spavars=size(@variables)-1;
643  int range=n*spavars;
644  if(f==0)
645    {
646      return(unitmat(range));
647    }
648  matrix M[range][range];
649  int length=size(f);
650  intvec max=maxexp(f);
651  int i;
652  intvec shiftback;
653  intvec vzero;
654  intvec vmax;
655  intvec shiftforward;
656  for(i=1;i<=size(max);i++)
657    {
658     shiftback[i]=int(floor(max[i]/2));
659     vzero[i]=0;
660     vmax[i]=n-1;
661     shiftforward[i]=0;
662    }
663  kill i;
664  int i;
665  for(i=1;i<=range;i++)
666    {
667      poly g=f;
668    }
669  kill i;
670}
671//////////////////////////////////////
672static proc floor(n)
673{
674  number h1=numerator(n);
675  number h2=denominator(n);
676  return((h1- (h1 mod h2))/h2);
677}
678/////////////////////////////////////
679static proc maxexp(poly f)
680{
681  int length=size(f);
682  poly g=f;
683  intvec v;
684  int i;
685  for(i=1;i<size(@variables);i++)
686    {
688    }
689  while(g!=0)
690    {
691      int j;
692      for(j=1;j<size(@variables);j++)
693  {
696  }
697      kill j;
698    }
699  return(v);
700}
701////////////////////////////////////
702static proc maximal(int n1,int n2)
703{
704  if(n1>n2)
705    {
706      return(n1);
707    }
708  else
709    {
710      return(n2);
711    }
712}
713////////////////////////////////////
714static proc minimal(int n1, int n2)
715{
716  return(-maximal(-n1,-n2));
717}
718////////////////////////////////////
719static proc MatrixEntry(int n, intvec v)
720{
721  int j;
722  int entry;
723  int spavar=size(@variables)-1;
724  for(j=1;j<=spavar;j++)
725    {
726      entry=entry+v[j]*n**(spavar-j);
727    }
728  entry=entry+1;
729  return(entry);
730}
731//////////////////////////////////
732static proc CompareVec(intvec ToTest, intvec Reference)//1 if ToTest>=Reference, 0 else
733{
734  int i;
735  for(i=1;i<=size(@variables)-1;i++)
736    {
737      if(ToTest[i+2]<Reference[i])
738  {
739    return(0);
740  }
741    }
742  return(1);
743}
744/////////////////////////////////
745static proc MaxVecZero(intvec ToTest, intvec Reference) //KGV, size=size(input)
746{
747  int length=size(ToTest);
748  int i;
749  intvec Maximum;
750  for(i=1;i<=length;i++)
751    {
752      Maximum[i]=maximal(ToTest[i],Reference[i]);
753    }
754  return(Maximum);
755}
756//////////////////////////////////
757proc matrixsystem(list Matrices,list Approx)
758"USAGE:   matrixsytem(M,A); where the M=Mt,M1,..,Mn is a list with square matrices of the same dimension as entries, and A=At,A1,..,An gives the corresponding approximations for the several variables (t,x1,..,xn) as vector. Intended to solve Mt*U_t + M1*U_x1+..+Mn*U_xn=0 as a linear sytem of partial differential equations numerically by a finite difference scheme;
759RETURN:   type matrix; gives back the matrices B1,B2 that represent the finite difference scheme, partitioned into different time levels in the form:  B1*U(t=N)=B2*U(t<N), where N is the maximal occurring degree (timelevel) of t.
760EXAMPLE:  example matrixsystem; shows an example
761"
762{
763  if(size(Matrices)>size(@variables) or size(Matrices)!=size(Approx))
764    {
765      ERROR("Check number of variables: it must hold  #(matrices)<= #(spatial variables)+1  !!! ");
766    }
767  if(size(Matrices)!=size(Approx))
768    {
769      ERROR("Every variable needs EXACTLY ONE approximation rule, i.e. #(first argument) =#(second argument) ! ");
770    }
772  int N=size(Matrices);
773  int i;
774  for(i=2;i<=N;i++)
775    {
777    }
778  kill i;
779  poly LCM=lcm(Mon);
780  matrix M[nrows(Matrices[1])][ncols(Matrices[1])];
781  int i;
782  for(i=1;i<=size(Matrices);i++)
783    {
785    }
786  kill i;
787  return(M);
788}
789example
790{
791  "EXAMPLE:";echo=2;
792  list D="Ut","Ux","Uy","U";
793  list V="t","x","y";
794  list P="a","b";
795  setinitials(V,D,P);
796  list Mat=unitmat(2),unitmat(2);
797  list Appr=forward(Ut,U,t),forward(Ux,U,x);
798  matrixsystem(Mat,Appr);
799}
800//////////////////////////////////
801proc timestep(matrix M)
802"USAGE:   timestep(M); M a square matrix with polynomials over the basering as entries;
803RETURN:   type list; gives two  matrices M1,M2 that are the splitting of M with respect to the degree of the variable t in the entries, where the first list-entry M1 consists of the polynomials of the highest timelevel and M2 of the lower levels in the form:  M=0  =>  M1=M2,    i.e. M1-M2=M
804NOTE:     intended to be used for the finite-difference-scheme-construction and partition into the several time steps
805EXAMPLE:  example timestep; shows an example
806"
807{
808  int N=nrows(M);
809  int i;
810  int maxdegT;
811  for(i=1;i<=N;i++)
812    {
813      int j;
814      for(j=1;j<=N;j++)
815  {
816    poly f=M[i,j];
817    int k;
818    for(k=1;k<=size(f);k++)
819      {
821    {
823    }
825      }
826    kill f;
827    kill k;
828  }
829      kill j;
830    }
831  kill i;
832  matrix highT[nrows(M)][nrows(M)];
833  vector leftside=0;
834  int GenIndex=0;
835  int i;
836  for(i=1;i<=N;i++)
837    {
838      int j;
839      for(j=1;j<=N;j++)
840  {
841    poly f=M[i,j];
842    int k;
843    for(k=1;k<=size(f)+1;k++)
844      {
846    {
847      GenIndex++;
849      leftside=leftside+highT[i,j]*gen(GenIndex);
850    }
852      }
853    kill k;
854    kill f;
855  }
856      kill j;
857    }
858  kill i;
859  matrix tUpper=highT;
860  matrix tLower=-1*(M-tUpper);
861  tUpper=tUpper/content(leftside);
862  tLower=tLower/content(leftside);
863  list L=tUpper,tLower;
864  return(L);
865}
866example
867{
868  "EXAMPLE:"; echo=2;
869  list D="Ut","Ux","Uy","U";
870  list V="t","x","y";
871  list P="a","b";
872  setinitials(V,D,P);
873  list Mat=unitmat(2),unitmat(2);
874  list Apr=forward(Ut,U,t),forward(Ux,U,x);
875  matrixsystem(Mat,Apr);
876  timestep(_);
877}
878//////////////////////////////////
879proc fouriersystem(list Matrices, list Approx)
880"USAGE:   fouriersystem(M,A); M a list of matrices, A a list of approximations;
881RETURN:   type list; each entry is some matrix obtained by performing the substitution of the single approximations into the system of pde's, partitioning the equation into the several timesteps and fouriertransforming these parts
882EXAMPLE:  example fouriersystem; shows an example
883"
884{
885  matrix M=matrixsystem(Matrices,Approx);
886  matrix T1=timestep(M)[1];
887  matrix T0=timestep(M)[2];
888  int i;
889  for(i=1;i<=nrows(M);i++)
890    {
891      int j;
892      for(j=1;j<=nrows(M);j++)
893  {
894    T1[i,j]=errormap(T1[i,j]);
895    T1[i,j]=VarToPar(T1[i,j]);
896    T0[i,j]=errormap(T0[i,j]);
897    T0[i,j]=VarToPar(T0[i,j]);
898  }
899      kill j;
900    }
901  kill i;
902  ideal EV1=eigenvals(T1)[1];
903  ideal EV0=eigenvals(T0)[1];
904  list L=list(T1,T0),list(EV1,EV0);
905  def N1=size(EV1);
906  def N0=size(EV0);
907  list CV1;
908  list CV0;
909  int i;
910  for(i=1;i<=N1;i++)
911  {
912  CV1[i]=VarToPar(EV1[i]);
913  if(content(CV1[i])==CV1[i])
914  {
915  CV1[i]=content(CV1[i]);
916  CV1[i]=VarToPar(ComplexValue(numerator(CV1[i])))/VarToPar(ComplexValue(denominator(CV1[i])));
917  }
918  }
919  kill i;
920  int i;
921  for(i=1;i<=N0;i++)
922  {
923  CV0[i]=VarToPar(EV0[i]);
924  if(content(CV0[i])==CV0[i])
925  {
926  CV0[i]=content(CV0[i]);
927  CV0[i]=VarToPar(ComplexValue(numerator(CV0[i])))/VarToPar(ComplexValue(denominator(CV0[i])));
928  }
929  }
930  kill i;
931  list CV=list(CV1,CV0);
932  L=L,CV;
933  return(L);
934}
935example
936{
937  "EXAMPLE:"; echo = 2;
938  list D="Ut","Ux","Uy","U";
939  list V="t","x","y";
940  list P="a","b";
941  setinitials(V,D,P);
942  matrix M[2][2]=0,-a,-a,0;
943  list Mat=unitmat(2),M;
944  list Appr=forward(Ut,U,t),trapezoid(Ux,U,x);
945  def s=fouriersystem(Mat,Appr);s;
946}
947//////////////////////////////////
948proc PartitionVar(poly f,int n)
949"USAGE:   PartitionVar(f); f a poly in the basering;
950RETURN:   type poly; gives back a list L=f1,f2 obtained by the partition of f into two parts f1,f2 with deg_var_n(f1) >0  deg_var_n(f2)=0
951EXAMPLE:  example PartitionVar; shows an example
952"
953{
954  if(n>=nvars(basering))
955    {
956      ERROR("this variable does not exist in the current basering");
957    }
958  int i;
959  poly partition=0;
960  poly g=f;
961  for(i=1;i<=size(f);i++)
962    {
964  {
966  }
968    }
969  list L=partition,f-partition;
970  return(L);
971}
972example
973{
974  "EXAMPLE:"; echo = 2;
975  list D="Ut","Ux","Uy","U";
976  list V="t","x","y";
977  list P="a","b";
978  setinitials(V,D,P);////does not show the ring, since there is no output
979  basering;///does show the ring
980  poly f=t**3*cx**2-cy**2*dt+i**3*sx;
981  PartitionVar(f,1); ////i is the first variable
982}
983//////////////////////////////////
984proc ComplexValue(poly f)
985"USAGE:   ComplexValue(f); f a poly in the basering;
986RETURN:   type poly; gives back the formal complex-value of f, where var(1) is redarded as the imaginary unit. Does only make sence, if the proc <setinitials> is executed before -> nvars <= npars
987EXAMPLE:  example ComplexValue; shows an example
988"
989{
990  poly g=ParToVar(f);
991  def L=PartitionVar(g,1);
992  poly f1=subst(L[1],var(1),1);
993  poly f2=L[2];
994  poly result=reduce(f1**2+f2**2,std(0));
995  return(result);
996}
997example
998{
999  "EXAMPLE:"; echo = 2;
1000  list D="Ut","Ux","Uy","U";
1001  list V="t","x","y";
1002  list P="a","b";
1003  setinitials(V,D,P);////does not show the ring, as there is  no output
1004  basering;///does show the ring
1005  poly f=t**3*cx**2-cy**2*dt+i**3*sx;
1006  f;
1007  VarToPar(f);
1008}
1009//////////////////////////////////
1010proc VarToPar(poly f)
1011"USAGE:   VarToPar(f); f a poly in the basering;
1012RETURN:   type poly; gives back the poly obtained by substituting var(i) by par(i), for all variables. Does only make sence, if the proc <setinitials> is executed before -> nvars <= npars;
1013EXAMPLE:  example VarToPar; shows an example
1014"
1015{
1016  int N=nvars(basering);
1017  int i;
1018  def g=f;
1019  for(i=1;i<=N;i++)
1020    {
1021      g=subst(g,var(i),par(i));
1022    }
1023  return(g);
1024}
1025example
1026{
1027  "EXAMPLE:"; echo = 2;
1028  list D="Ut","Ux","Uy","U";
1029  list V="t","x","y";
1030  list P="a","b";
1031  setinitials(V,D,P);////does not show the ring, as there is  no output
1032  basering;///does show the ring
1033  poly f=t**3*cx**2-cy**2*dt+i**3*sx;
1034  f;
1035  VarToPar(f);
1036}
1037/////////////////////////////////////
1038proc ParToVar(poly f)
1039"USAGE:   ParToVar(f); f a poly in the basering;
1040RETURN:   type poly; gives back the poly obtained by substituting par(i) by var(i), for the first nvars(basering parameters. Does only make sence, if setinitials is executed before -> nvars <= npars. Is the opposite action to VarToPar, see example ParToVar;
1041EXAMPLE:  example ParToVar; shows an example
1042"
1043{
1044  int N=nvars(basering);
1045  int i;
1046  number g=number(VarToPar(f));
1047  number denom=denominator(g);
1048  g=denom*g;
1049  def gg=subst(g,par(1),var(1));
1050  for(i=2;i<=N;i++)
1051    {
1052      gg=subst(gg,par(i),var(i));
1053    }
1054  return(gg/denom);
1055}
1056example
1057{
1058  "EXAMPLE:"; echo = 2;
1059  list D="Ut","Ux","Uy","U";
1060  list V="t","x","y";
1061  list P="a","b";
1062  setinitials(V,D,P);////does not show the ring, as there is  no output
1063  basering;///does show the ring
1064  poly f=t**3*cx**2-cy**2*dt+i**3*sx/dt*dx;
1065  f;
1066  def g=VarToPar(f);
1067  g;
1068  def h=ParToVar(g);
1069  h==f;
1070}
1071/////////////////////////////////////////
1073"USAGE:   qepcad(f); f a poly in the basering;
1074RETURN:   type list; gives back some constraints that are equivalent to f<1 (computed by QEPCAD);
1075EXAMPLE:  example qepcad; shows an example
1076"
1077{
1078  // how to test, whether QEPCAD is installed?
1082  if(denominator(content(f))==1)
1083  {
1084    poly g=f-1;
1085  }
1086  else
1087  {
1088    if(f==content(f))
1089    {
1090      poly g=f*denominator(content(f))-1*denominator(content(f));
1091      g=ParToVar(g);
1092      g=reduce(g,std(0));
1093    }
1094    else
1095    {
1096      poly g=cleardenom(f)-1/content(f);
1097      g=ParToVar(g);
1098      g=reduce(g,std(0));
1099    }
1100  }
1105  string s1="[trial]";    //description
1106  string s2=varlist();    //the variables
1107  string s3=nfreevars();  //number of free variables
1108  string s4=aquantor()+"["+writepoly(g)+rel()+"]."; //the input prenex formula
1109  string s5=projection();
1110  string s6=projection();
1111  string s7=choice();
1112  string s8=solution();
1113  write(l1,s1,s2,s3,s4,s5,s6,s7,s8);
1116  print(output,"%s");
1117  if(size(output)==0)
1118  {
1119    return("Try manually");    //maybe too few cells
1120  }
1121  if(find(output,"FALSE")!=0)
1122  {
1123    return("FALSE");
1124  }
1125  if(find(output,"WARNING")!=0)
1126  {
1127    return("WARNING! Try manually");
1128  }
1129  else
1130  {
1131    string strpolys=findthepoly(output);
1132    list lpolys=listpolynew(strpolys);
1133    return(lpolys);
1134  }
1137
1138}
1139example
1140{
1141  "EXAMPLE:"; echo = 2;
1142  list D="Ux","Ut","U";
1143  list P="a";
1144  list V="t","x";
1145  setinitials(V,D,P);
1146  def s1=scheme(u(Ut)+a*u(Ux),laxfrX(Ux,U,x),laxfrT(Ut,U,x));
1147  s1;
1148  def s2=errormap(s1);
1149  s2;
1150  def s3=ComplexValue(s2);s3;
1152}
1153///////////////////////////////////////////
1154proc createQCfilter()
1155{
1156  // writes the following to the file qepcadfilter.pl
1157  // is there already such a file? remove it!
1160  write(l, "#!/usr/bin/perl");
1161  write(l, "\$flag = 0;");
1162  write(l, "\$res = \"\";");
1163  write(l,"while(<>)");
1164  write(l,"{ if (\$_ =~ /Warning|WARNING|warning|Error|error|ERROR/) { print \$_; }");
1165  write(l,"elsif (\$_ =~ /An\ equivalent/) { \$flag = 1; }");
1166  write(l,"elsif (\$flag == 1 && \$_ ne \"\n\") { print \$_; \$flag = 0; } }");
1167}
1168
1169///////////////////////////////////////////
1172RETURN:  list
1173PURPOSE: gives back some constraints that are equivalent to the
1174eigenvalues of the matrices in the list l being < 1 (computed by QEPCAD)
1175EXAMPLE:  example qepcadsystem; shows an example
1176"
1177{
1178  // how to test, whether QEPCAD is installed?
1186  string s1="[trial]";    //description
1187  string s2=varlist();    //the variables
1188  string s3=nfreevars();  //number of free variables
1189  string thepolys;
1190  int n2=size(l[2]);
1191  int count;
1192  int i;
1193  list lpolys;
1194  int j;
1195  for(j=1;j<=n2;j++)
1196  {
1197    count++;
1198    poly g2=ParToVar(l[2][j]);
1199    if(denominator(content(g2))==1)
1200    {
1201      lpolys[count]=writepoly(ParToVar(reduce(g2-1,std(0))))+rel();
1202    }
1203    else
1204    {
1205      if(g2==content(g2))
1206      {
1207        g2=g2*denominator(content(g2))-1*denominator(content(g2));
1208        g2=ParToVar(g2);
1209        g2=reduce(g2,std(0));
1210        lpolys[count]=writepoly(g2)+rel();
1211      }
1212      else
1213      {
1214        lpolys[count]=writepoly(reduce(ParToVar(cleardenom(g2)-1/content(g2)),std(0)))+rel();
1215      }
1216    }
1217    kill g2;
1218  }
1219  kill j;
1220  int k;
1221  for(k=1;k<=size(lpolys);k++)
1222  {
1223    thepolys=thepolys+lpolys[k];
1224    if(k<size(lpolys))
1225    {
1226      thepolys=thepolys+print(" /","s%")+print("\\ ","s%");
1227    }
1228  }
1229  string s4=aquantor()+"["+thepolys+"]."; //the input prenex formula
1230  string s5=projection();
1231  string s6=projection();
1232  string s7=choice();
1233  string s8=solution();
1234  write(l1,s1,s2,s3,s4,s5,s6,s7,s8);
1237  print(output,"%s");
1238  if(size(output)==0)
1239  {
1240    ERROR("Try manually");    //maybe too few cells
1241  }
1242  if(find(output,"FALSE")!=0)
1243  {
1244    ERROR("FALSE");
1245  }
1246  if(find(output,"WARNING")!=0)
1247  {
1248    ERROR("WARNING! Try manually");
1249  }
1250  else
1251  {
1252    string strpolys=findthepoly(output);
1253    list llpolys=listpolynew(strpolys);
1254    return(llpolys);
1255  }
1258}
1259example
1260{
1261  "EXAMPLE:"; echo = 2;
1262  list D="Ut","Ux","Uy","U";
1263  list V="t","x","y";
1264  list P="a","b";
1265  setinitials(V,D,P);
1266  matrix M[2][2]=0,-a,-a,0;
1267  list Mat=unitmat(2),M;
1268  list Appr=forward(Ut,U,t),forward(Ux,U,x);
1269  //matrixsystem(Mat,Appr);
1270  //timestep(_);
1271  fouriersystem(Mat,Appr);
1273}
1274///////////////////////////////////////////
1275static proc substbracketstar(string s)
1276{
1277  int i;
1278  int k;
1279  int index=1;
1280  string finish=s;
1281  for(k=1;k<=size(s);k++)
1282  {
1283    if(finish[1]=="(" or finish[1]=="*" or finish[1]==" ")
1284    {
1285      kill finish;
1286      index=index+1;
1287      string finish=s[index..size(s)];
1288    }
1289  }
1290  for(i=1;i<=size(finish);i++)
1291  {
1292    if(finish[i]=="*" or finish[i]=="(" or finish[i]== ")")
1293    {
1294      finish[i]=" ";
1295    }
1296  }
1297  return(finish);
1298}
1299
1300////////////////////////////////////
1301static proc distribution(string SUM, string MON)
1302{
1303  string sum=substbracketstar(SUM);
1304  string mon=substbracketstar(MON);
1305  string result;
1306  list signs;
1307  list p;
1308  int i;
1309  int j;
1310  int init;
1311  for(i=2;i<=size(sum);i++)
1312  {
1313    if(sum[i]=="+" or sum[i]=="-")
1314    {
1315      j++;
1316      p[j]=i;
1317    }
1318  }
1319  if(j==0)
1320  {
1321    if(sum[1]!="-")
1322    {
1323      result=sum+" "+" "+mon;
1324      result="+"+" "+result;
1325    }
1326    else
1327    {
1328      result=sum+" "+mon;
1329    }
1330  }
1331  else
1332  {
1333    int l;
1334    int anfang;
1335    if(sum[1]=="-")
1336    {
1337      result="-"+" "+result;
1338      anfang=2;
1339    }
1340    else
1341    {
1342      result="+"+" "+result;
1343      if(sum[1]=="+")
1344      {
1345        anfang=2;
1346      }
1347      else
1348      {
1349        anfang=1;
1350      }
1351    }
1352    for(l=1;l<=j;l++)
1353    {
1354      string a;
1355      int k;
1356      for(k=anfang;k<=p[l]-1;k++)
1357      {
1358        a=a+sum[k];
1359      }
1360      result=result+" "+a+" "+mon+" "+sum[p[l]];
1361      anfang=p[l]+1;
1362      kill a;
1363      kill k;
1364    }
1365    if(p[j]<size(sum))
1366    {
1367      int kk;
1368      string aa;
1369      for(kk=anfang;kk<=size(sum);kk++)
1370      {
1371        aa=aa+sum[kk];
1372      }
1373      result=result+" "+aa+" "+mon;
1374      kill aa;
1375      kill kk;
1376    }
1377    else
1378    {
1379      int kkk;
1380      string aaa;
1381      for(kkk=anfang;kkk<size(sum);kkk++)
1382      {
1383        aaa=aaa+sum[kkk];
1384      }
1385      result=result+" "+aaa+" "+mon;
1386      kill aaa;
1387      kill kkk;
1388    }
1389  }
1390  return(result);
1391}
1392/////////////////////////////////////////////////////////////////
1393static proc writepoly(poly f)
1394{
1395  poly g=f;
1396  string lc;
1397  string lm;
1398  string strpoly;
1399  string intermediate;
1400  int n=size(f);
1401  int i;
1402  for(i=1;i<=n;i++)
1403  {
1404
1407    intermediate=distribution(lc,lm);
1408    strpoly=strpoly+" "+intermediate;
1410  }
1411  return(strpoly);
1412}
1413///////////////////////////////////////////////////////////////
1414static proc varlist()
1415{
1416  poly p1=par(2+size(@variables));
1417  p1=ParToVar(p1);
1418  string name=print(p1,"%s");
1419  string p="("+name+",";
1420  int i;
1421  for(i=3+size(@variables);i<=npars(basering);i++)
1422  {
1423    p1=ParToVar(par(i));
1424    name=substbracketstar(print(p1,"%s"));
1425    p=p+name+",";
1426  }
1427  p1=ParToVar(par(2));
1428  name=print(p1,"%s");
1429  p=p+name+")";
1430  return(p);
1431}
1432///////////////////////////////////////////////////////
1433static proc normalization()
1434{
1435  int n1=npars(basering)-size(@variables);
1436  int n2=npars(basering)-n1;
1437  string assumption="assume["+" "+substbracketstar(print(par(n1+1),"%s"))+" >0";
1438  int i;
1439  for(i=2;i<=n2;i++)
1440  {
1441    string s=" /\\"+" "+substbracketstar(print(par(n1+i),"%s"))+" >0";
1442    assumption=assumption+" "+s;
1443    kill s;
1444  }
1445  assumption=assumption+" ]"+newline+"go";
1446  return(assumption);
1447}
1448///////////////////////////////////////////////////////
1449static proc projection()
1450{
1451  string s="go";
1452  return(s);
1453}
1454///////////////////////////////////////////////////////
1455static proc choice()
1456{
1457  string s="go";
1458  return(s);
1459}
1460///////////////////////////////////////////////////////
1461static proc solution()
1462{
1463  string s="go";
1464  return(s);
1465}
1466///////////////////////////////////////////////////////
1467static proc nfreevars()
1468{
1469  int n=npars(basering)-size(@variables)-1;
1470  string no=print(n,"%s");
1471  return(no);
1472}
1473/////////////////////////////////////////////////////////
1474static proc aquantor()
1475{
1476  string s="(A "+print(var(2),"%s")+")";
1477  return(s);
1478}
1479////////////////////////////////////////////////////////
1480static proc rel()
1481{
1482  return("<0");
1483}
1484/////////////////////////////////////////////////////////
1485static proc rm()
1486{
1487  system("sh","rm -f dummy");
1489}
1490////////////////////////////////////////////////////////
1491static proc findthepoly(string word)
1492{
1493  int init=1;
1494  int index=find(word,newline);
1495  if(index==size(word) or index== 0)
1496  {
1497    return(word);
1498  }
1499  else
1500  {
1501    while(index<size(word))
1502    {
1503      init=index;
1504      index=find(word,newline,index+1);
1505    }
1506  }
1507  string thepoly=word[(init+1)..(size(word)-1)];
1508  return(thepoly);
1509}
1510//////////////////////////////////////////////////
1511static proc listpolynew(string thepoly)
1512{
1513  string p=thepoly;
1514  intvec v;
1515  p=TestAndDelete(p,"[")[1];
1516  p=TestAndDelete(p,"]")[1];    //remove the brackets
1517  string AND="/"+"\\";
1518  string OR="\\"+"/";
1519  list thesigns=AND,OR,"~","<==>","==>","<==";///these can occur in QEPCAD
1520  int i;
1521  for(i=1;i<=size(thesigns);i++)
1522  {
1523    list l=TestAndDelete(p,thesigns[i]);
1524    p=l[1];
1526    kill l;
1527  }
1528  intvec w=rightorder(v);
1529  int N=size(v);
1530  list lstrings;
1531  if(size(w)==1 and w[1]==0)
1532  {
1533    N=0;
1534    lstrings[1]=p;
1535  }
1536  else
1537  {
1538    string s1=p[1..w[1]-1];
1539    lstrings[1]=s1;
1540    int j;
1541    for(j=1;j<N;j++)
1542    {
1543      string s2=p[w[j]..w[j+1]-1];
1544      lstrings[j+1]=s2;
1545      kill s2;
1546    }
1547    string s3=p[w[N]..size(p)];
1548    lstrings[N+1]=s3;
1549  }
1550  list cutstrings;
1551  int k;
1552  list Id;
1553  for(k=1;k<=N+1;k++)
1554  {
1555    cutstrings[k]=cutoffrel(lstrings[k]);
1556    Id[k]=translatenew(cutstrings[k]);
1557  }
1558  return(Id);
1559
1560
1561}
1562/////////////////////////////////////////////////
1563static proc translatenew(string word)
1564{
1565  int n=size(word);
1566  string subword=word[1..n];
1567  string strpoly="poly f=";
1568  list linterim;
1569  int i;
1570  for(i=1;i<=n;i++)
1571  {
1572    if(i<n and subword[i]==" " and subword[i+1]!="+" and subword[i+1]!="-" and subword[i+1]!=" " and subword[i-1]!="+" and subword[i-1]!="-" and subword[i-1]!=" " and i>1)
1573    {
1574      linterim[i]="*";
1575    }
1576    else
1577    {
1578      linterim[i]=word[i];
1579    }
1580    strpoly=strpoly+print(linterim[i],"%s");
1581  }
1582  strpoly=strpoly+";";
1583  execute(strpoly);
1584  return(f);
1585}
1586//////////////////////////////////////////////////
1587static proc rightorder(intvec v)
1588////////////orders the entries of an intvec from small to bigger
1589{
1590  list ldown=intveclist(v);
1591  list lup;
1592  intvec v1;
1593  int counter;
1594  int min;
1595  int i;
1596  for(i=1;i<=size(v);i++)
1597  {
1598    min=Min(listintvec(ldown));
1599    lup[i]=min;
1600    counter=listfind(ldown,min);
1601    ldown=delete(ldown,counter);
1602  }
1603  intvec result=listintvec(lup);
1604  return(result);
1605}
1606/////////////////////////////////////////////////
1607static proc listfind(list l, int n)
1609{
1610  int i;
1611  intvec counter;
1612  int j=1;
1613  for(i=1;i<=size(l);i++)
1614  {
1615    if(n==l[i])
1616    {
1617      counter[j]=i;
1618      j++;
1619    }
1620  }
1621  int result=Min(counter);
1622  return(result);
1623}
1624//////////////////////////////////////////////////
1625static proc cutoffrel(string str)
1626/////////cut off the relations "/=,<= etc." from output-string
1627{
1628  list rels="<=",">=","<",">","/=","=";  //these are all of qepcad's relations
1629  int i;
1630  list l;
1631  for(i=1;i<=size(rels);i++)
1632  {
1633    l[i]=find(str,rels[i]);
1634  }
1635  intvec v=listintvec(l);
1637  int position=Min(v);
1638  string result;
1639  if(position==0)
1640  {
1641    result=str;
1642  }
1643  else
1644  {
1645    result=str[1..position-1];
1646  }
1647  return(result);
1648}
1649//////////////////////////////////////////////////
1650static proc addintvec(intvec v1, intvec v2)
1651///////concatenates two intvecs and deletes occurring zeroentries, output intvec
1652{
1653  list lresult;
1654  int i;
1655  list l1;
1656  list l2;
1657  for(i=1;i<=size(v1);i++)
1658  {
1659    l1[i]=v1[i];
1660  }
1661  kill i;
1662  int k;
1663  for(k=1;k<=size(v2);k++)
1664  {
1665    l2[k]=v2[k];
1666  }
1667  lresult=l1+l2;
1668  intvec result;
1669  int j;
1670  int counter=1;
1671  for(j=1;j<=size(lresult);j++)
1672  {
1673    if(lresult[j]!=0)
1674    {
1675      result[counter]=lresult[j];
1676      counter++;
1677    }
1678  }
1679  return(result);
1680
1681}
1682/////////////////////////////////////////////////
1683static proc intveclist(intvec v)
1684//////returns the intvec as list
1685{
1686  list l;
1687  int i;
1688  for(i=size(v);i>=1;i--)
1689  {
1690    l[i]=v[i];
1691  }
1692  return(l);
1693}
1694//////////////////////////////////////////////////
1695static proc listintvec(list l)
1696//////////returns the list as intvec: only integer-entries allowed
1697{
1698  intvec v;
1699  int i;
1700  for(i=size(l);i>=1;i--)
1701  {
1702    v[i]=l[i];
1703  }
1704  return(v);
1705}
1706//////////////////////////////////////////////////
1707static proc TestAndDelete(string str, string booleansign)
1708{
1709  int decide=find(str,booleansign);
1710  intvec v;
1711  v[1]=decide;
1712  while(decide!=0 and decide<size(str))//if booleansign ist not contained in str
1713  {
1714    int n=size(v);
1715    decide=find(str,booleansign,decide+1);
1716    if(decide!=0)
1717    {
1718      v[n+1]=decide;
1719    }
1720    kill n;
1721  }
1722  if(size(v)==1 and v[1]==0)/////i.e. booleansign NOT in str
1723  {
1724    list result=str,v;
1725    return(result);
1726  }
1727  else
1728  {
1729    list l;
1730    string p;
1731    int i;
1732    int j;
1733    for(i=1;i<=size(str);i++)///copy the string str into l
1734    {
1735      l[i]=str[i];
1736    }
1737    kill i;
1738    for(j=1;j<=size(v);j++)///replace booleansign by empty char
1739    {
1740      int k;
1741      for(k=0;k<=size(booleansign)-1;k++)
1742      {
1743        l[v[j]+k]=" ";
1744      }
1745      kill k;
1746    }
1747    int s;
1748    for(s=1;s<=size(l);s++)///copy back
1749    {
1750      p=p+string(l[s]);
1751    }
1752    kill s;
1753    kill l;
1754    list result=p,v;
1755    return(result);
1756  }
1757}
1758//////////////////////////////////////////////////
1759static proc nchoice(int n1,int n2)
1760{
1761  int N2=maximal(n1,n2);
1762  int N1=minimal(n1,n2);
1763  if(N1==0)
1764  {
1765    return(N2);
1766  }
1767  else
1768  {
1769    return(N1);
1770  }
1771}
Note: See TracBrowser for help on using the repository browser.