# source:git/Singular/LIB/findiff.lib@f0a801b

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