source: git/Singular/LIB/findifs.lib @ 1288ef

spielwiese
Last change on this file since 1288ef was 1288ef, checked in by Motsak Oleksandr <motsak@…>, 15 years ago
resurection git-svn-id: file:///usr/local/Singular/svn/trunk@11634 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 17.0 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: findifs.lib,v 1.6 2009-04-06 17:43:33 motsak Exp $";
3category="Applications";
4info="
5LIBRARY: findifs.lib  Tools for the finite difference schemes
6AUTHORS: Viktor Levandovskyy,     levandov@risc.uni-linz.ac.at
7
8THEORY: We provide the presentation of difference operators in a polynomial,
9        semi-factorized and a nodal form. Running @code{findifs_example();}
10        will show how we generate finite difference schemes of linear PDE's
11        from given approximations.
12
13PROCEDURES:
14decoef(P,n);  decompose poly P into summands with respect to the number n
15difpoly2tex(S,P[,Q]); present the difference scheme in the nodal form
16
17AUXILIARY PROCEDURES:
18findifs_example();  containes a guided explanation of our approach
19exp2pt(P[,L]); convert a polynomial M into the TeX format, in nodal form
20mon2pt(P[,L]); convert a monomial M into the TeX format, in nodal form
21texcoef(n);  converts the number n into TeX
22npar(n); search for 'n' among the parameters and returns its number
23magnitude(P); compute the square of the magnitude of a complex expression
24replace(s,what,with);  replace in s all the substrings with a given string
25xchange(w,a,b); exchange two substrings of a string
26par2tex(s);     convert special characters to TeX in s
27";
28
29
30LIB "latex.lib";
31LIB "poly.lib";
32
33// 1. GLOBAL ASSUME: in the ring we have first Tx, then Tt: [FIXED, not needed anymore]!
34// 2. map vars other than Tx,Tt to parameters instead or just ignore them [?]
35// 3. clear the things with brackets
36// 4. todo: content resp lcmZ, gcdZ
37
38proc xchange(string where, string a, string b)
39"USAGE:  xchange(w,a,b);  w,a,b strings
40RETURN:  string
41PURPOSE:  exchanges substring 'a' with a substring 'b' in the string w
42NOTE:
43EXAMPLE: example xchange; shows examples
44"{
45  // replaces a<->b in where
46  // assume they are of the same size [? seems to work]
47  string s = "H";
48  string t;
49  t = replace(where,a,s);
50  t = replace(t,b,a);
51  t = replace(t,s,b);
52  return(t);
53}
54example
55{
56 " EXAMPLE:"; echo=2;
57 ring r = (0,dt,dh,A),Tt,dp;
58 poly p = (Tt*dt+dh+1)^2+2*A;
59 string s = texpoly("",p);
60 s;
61 string t = xchange(s,"dh","dt");
62 t;
63}
64
65proc par2tex(string s)
66"USAGE:  par2tex(s);  s a string
67RETURN:  string
68PURPOSE: converts special characters to TeX in s
69NOTE: the convention is the following:
70      'Tx' goes to 'T_x',
71      'dx' to '\\tri x' (the same for dt, dy, dz),
72      'theta', 'ro', 'A', 'V' are converted to greek letters.
73EXAMPLE: example par2tex; shows examples
74"{
75// s is a tex string with a poly
76// replace theta with \theta
77// A with \lambda
78// dt with \tri t
79// dh with \tri h
80// Tx with T_x, Ty with T_y
81// Tt with T_t
82// V with \nu
83// ro with \rho
84// dx with \tri x
85// dy with \tri y
86  string t = s;
87  t = replace(t,"Tt","T_t");
88  t = replace(t,"Tx","T_x");
89  t = replace(t,"Ty","T_y");
90  t = replace(t,"dt","\\tri t");
91  t = replace(t,"dh","\\tri h");
92  t = replace(t,"dx","\\tri x");
93  t = replace(t,"dy","\\tri y");
94  t = replace(t,"theta","\\theta");
95  t = replace(t,"A","\\lambda");
96  t = replace(t,"V","\\nu");
97  t = replace(t,"ro","\\rho");
98  return(t);
99}
100example
101{
102 " EXAMPLE:"; echo=2;
103 ring r = (0,dt,theta,A),Tt,dp;
104 poly p = (Tt*dt+theta+1)^2+2*A;
105 string s = texfactorize("",p);
106 s;
107 par2tex(s);
108 string T = texfactorize("",p*(-theta*A));
109 par2tex(T);
110}
111
112proc replace(string s, string what, string with)
113"USAGE:  replace(s,what,with);  s,what,with strings
114RETURN:  string
115PURPOSE: replaces in 's' all the substrings 'what' with substring 'with'
116NOTE:
117EXAMPLE: example replace; shows examples
118"{
119  // clear: replace in s, "what" with "with"
120  int ss = size(s);
121  int cn = find(s,what);
122  if ( (cn==0) || (cn>ss))
123  {
124    return(s);
125  }
126  int gn = 0; // global counter
127  int sw = size(what);
128  int swith = size(with);
129  string out="";
130  string tmp;
131  gn  = 0;
132  while(cn!=0)
133  {
134//    "cn:";    cn;
135//    "gn";    gn;
136    tmp = "";
137    if (cn>gn)
138    {
139      tmp = s[gn..cn-1];
140    }
141//    "tmp:";tmp;
142//    out = out+tmp+" "+with;
143    out = out+tmp+with;
144//    "out:";out;
145    gn  = cn + sw;
146    if (gn>ss)
147    {
148      // ( (gn>ss) || ((sw>1) && (gn >= ss)) )
149      // no need to append smth
150      return(out);
151    }
152//     if (gn == ss)
153//     {
154
155//     }
156    cn  = find(s,what,gn);
157  }
158  // and now, append the rest of s
159  //  out = out + " "+ s[gn..ss];
160  out = out + s[gn..ss];
161  return(out);
162}
163example
164{
165 " EXAMPLE:"; echo=2;
166 ring r = (0,dt,theta),Tt,dp;
167 poly p = (Tt*dt+theta+1)^2+2;
168 string s = texfactorize("",p);
169 s;
170 s = replace(s,"Tt","T_t"); s;
171 s = replace(s,"dt","\\tri t"); s;
172 s = replace(s,"theta","\\theta"); s;
173}
174
175proc exp2pt(poly P, list #)
176"USAGE:  exp2pt(P[,L]);  P poly, L an optional list of strings
177RETURN:  string
178PURPOSE: convert a polynomial M into the TeX format, in nodal form
179ASSUME: coefficients must not be fractional
180NOTE: an optional list L contains a string, which will replace the default
181value 'u' for the discretized function
182EXAMPLE: example exp2pt; shows examples
183"{
184// given poly in vars [now Tx,Tt are fixed],
185// create Tex expression for points of lattice
186// coeffs must not be fractional
187  string varnm = "u";
188  if (size(#) > 0)
189  {
190    if (typeof(#[1])=="string")
191    {
192      varnm = string(#[1]);
193    }
194  }
195  //  varnm;
196  string rz,mz;
197  while (P!=0)
198  {
199    mz = mon2pt(P,varnm);
200    if (mz[1]=="-")
201    {
202      rz = rz+mz;
203    }
204    else
205    {
206      rz = rz + "+" + mz;
207    }
208    P = P-lead(P);
209  }
210  rz  = rz[2..size(rz)];
211  return(rz);
212}
213example
214{
215  " EXAMPLE:"; echo=2;
216  ring r = (0,dh,dt),(Tx,Tt),dp;
217  poly M = (4*dh*Tx^2+1)*(Tt-1)^2;
218  print(exp2pt(M));
219  print(exp2pt(M,"F"));
220}
221
222proc mon2pt(poly M, string V)
223"USAGE:  mon2pt(M,V);  M poly, V a string
224RETURN:  string
225PURPOSE: convert a monomial M into the TeX format, nodal form
226EXAMPLE: example mon2pt; shows examples
227"{
228  // searches for Tx, then Tt
229  // monomial to the lattice point conversion
230  // c*X^a*Y^b --> c*U^{n+a}_{j+b}
231  number cM = leadcoef(M);
232  intvec e = leadexp(M);
233  //  int a = e[2]; // convention: first Tx, then Tt
234  //  int b = e[1];
235  int i;
236  int a , b, c = 0,0,0;
237  int ia,ib,ic = 0,0,0;
238  int nv = nvars(basering);
239  string s;
240  for (i=1; i<=nv ; i++)
241  {
242    s = string(var(i));
243    if (s=="Tt") { a = e[i]; ia = i;}
244    if (s=="Tx") { b = e[i]; ib = i;}
245    if (s=="Ty") { c = e[i]; ic = i;}
246  }
247  //  if (ia==0) {"Error:Tt not found!"; return("");}
248  //  if (ib==0) {"Error:Tx not found!"; return("");}
249  //  if (ic==0)   {"Error:Ty not found!"; return("");}
250  //  string tc = texobj("",c); // why not texpoly?
251  string tc = texcoef(cM);
252  string rs;
253  if (cM==-1)
254  {
255    rs = "-";
256  }
257  if (cM^2 != 1)
258  {
259// we don't need 1 or -1 as coeffs
260//    rs = clTex(tc)+" ";
261//    rs = par2tex(rmDol(tc))+" ";
262    rs = par2tex(tc)+" ";
263  }
264// a = 0 or b = 0
265  rs = rs + V +"^{n";
266  if (a!=0)
267  {
268    rs = rs +"+"+string(a);
269  }
270  rs = rs +"}_{j";
271  if (b!=0)
272  {
273    rs = rs +"+"+string(b);
274  }
275  if (c!=0)
276  {
277    rs = rs + ",k+";
278    rs = rs + string(c);
279  }
280  rs = rs +"}";
281  return(rs);
282}
283example
284{
285 "EXAMPLE:"; echo=2;
286  ring r = (0,dh,dt),(Tx,Tt),dp;
287  poly M = (4*dh^2-dt)*Tx^3*Tt;
288  print(mon2pt(M,"u"));
289  poly N = ((dh-dt)/(dh+dt))*Tx^2*Tt^2;
290  print(mon2pt(N,"f"));
291  ring r2 = (0,dh,dt),(Tx,Ty,Tt),dp;
292  poly M  = (4*dh^2-dt)*Tx^3*Ty^2*Tt;
293  print(mon2pt(M,"u"));
294}
295
296proc npar(number n)
297"USAGE:  npar(n);  n a number
298RETURN:  int
299PURPOSE: searches for 'n' among the parameters and returns its number
300EXAMPLE: example npar; shows examples
301"{
302  // searches for n amongst parameters
303  // and returns its number
304  int i,j=0,0;
305  list L = ringlist(basering);
306  list M = L[1][2]; // pars
307  string sn = string(n);
308  sn = sn[2..size(sn)-1];
309  for (i=1; i<=size(M);i++)
310  {
311    if (M[i] == sn)
312    {
313      j = i;
314    }
315  }
316  if (j==0)
317  {
318    "Incorrect parameter";
319  }
320  return(j);
321}
322example
323{
324 "EXAMPLE:"; echo=2;
325  ring r = (0,dh,dt,theta,A),t,dp;
326  npar(dh);
327  number T = theta;
328  npar(T);
329  npar(dh^2);
330}
331
332proc decoef(poly P, number n)
333"USAGE:  decoef(P,n);  P a poly, n a number
334RETURN:  ideal
335PURPOSE:  decompose poly P into summands with respect
336to the  presence of a number n in the coefficients
337NOTE:    n is usually a parameter with no power
338EXAMPLE: example decoef; shows examples
339"{
340  // decomposes poly into summands
341  // wrt the presence of a number n in coeffs
342  // returns ideal
343  def br = basering;
344  int i,j=0,0;
345  int pos = npar(n);
346  if ((pos==0) || (P==0))
347  {
348    return(0);
349  }
350  pos = pos + nvars(basering);
351// map all pars except to vars, provided no things are in denominator
352  number con = content(P);
353  con = numerator(con);
354  P = cleardenom(P); //destroys content!
355  P = con*P; // restore the numerator part of the content
356  list M = ringlist(basering);
357  list L = M[1..4];
358  list Pars = L[1][2];
359  list Vars = L[2] + Pars;
360  L[1] = L[1][1]; // characteristic
361  L[2] = Vars;
362  // for non-comm things: don't need nc but graded algebra
363  //  list templ;
364  //  L[5] = templ;
365  //  L[6] = templ;
366  def @R = ring(L);
367  setring @R;
368  poly P = imap(br,P);
369  poly P0 = subst(P,var(pos),0);
370  poly P1 = P - P0;
371  ideal I = P0,P1;
372  setring br;
373  ideal I = imap(@R,I);
374  kill @R;
375  // check: P0+P1==P
376  poly Q = I[1]+I[2];
377  if (P!=Q)
378  {
379    "Warning: problem in decoef";
380  }
381  return(I);
382// substract the pure part from orig and check if n is remained there
383}
384example
385{
386 " EXAMPLE:"; echo=2;
387  ring r = (0,dh,dt),(Tx,Tt),dp;
388  poly P = (4*dh^2-dt)*Tx^3*Tt + dt*dh*Tt^2 + dh*Tt;
389  decoef(P,dt);
390  decoef(P,dh);
391}
392
393proc texcoef(number n)
394"USAGE:  texcoef(n);  n a number
395RETURN:  string
396PURPOSE:  converts the number n into TeX format
397NOTE: if n is a polynomial, texcoef adds extra brackets.
398@* Performs some space substitutions
399EXAMPLE: example texcoef; shows examples
400"{
401  // makes tex from n
402  // and uses substitutions
403  // if n is a polynomial, adds brackets
404  number D  = denominator(n);
405  int DenIsOne = 0;
406  if ( D==number(1) )
407  {
408    DenIsOne = 1;
409  }
410  string sd = texpoly("",D);
411  sd = rmDol(sd);
412  sd = par2tex(sd);
413  number N  = numerator(n);
414  string sn = texpoly("",N);
415  sn = rmDol(sn);
416  sn = par2tex(sn);
417  string sout="";
418  int i;
419  int NisPoly = 0;
420  if (DenIsOne)
421  {
422    sout = sn;
423    for(i=1; i<=size(sout); i++)
424    {
425      if ( (sout[i]=="+") || (sout[i]=="-") )
426      {
427        NisPoly = 1;
428      }
429    }
430    if (NisPoly)
431    {
432      sout = "("+sout+")";
433    }
434  }
435  else
436  {
437    sout = "\\frac{"+sn+"}{"+sd+"}";
438  }
439  return(sout);
440}
441example
442{
443 " EXAMPLE:"; echo=2;
444  ring r = (0,dh,dt),(Tx,Tt),dp;
445  number n1,n2,n3 = dt/(4*dh^2-dt),(dt+dh)^2, 1/dh;
446  n1; texcoef(n1);
447  n2; texcoef(n2);
448  n3; texcoef(n3);
449}
450
451static proc rmDol(string s)
452{
453  // removes $s and _no_ (s on appearance
454  int i = size(s);
455  if (s[1] == "$") { s = s[2..i]; i--;}
456  if (s[1] == "(") { s = s[2..i]; i--;}
457  if (s[i] == "$") { s = s[1..i-1]; i--;}
458  if (s[i] == ")") { s = s[1..i-1];}
459  return(s);
460}
461
462proc difpoly2tex(ideal S, list P, list #)
463"USAGE:  difpoly2tex(S,P[,Q]);  S an ideal, P and optional Q are lists
464RETURN:  string
465PURPOSE: present the difference scheme in the nodal form
466ASSUME: ideal S is the result of 'decoef' procedure
467NOTE: a list P may be empty or may contain parameters, which will not
468appear in denominators
469@* an optional list Q represents the part of the scheme, depending
470on other function, than the major part
471EXAMPLE: example difpoly2tex; shows examples
472"
473{
474  // S = sum s_i = orig diff poly or
475  // the result of decoef
476  // P = list of pars (numbers) not to be divided with, may be empty
477  // # is an optional list of polys, repr. the part dep. on "f", not on "u"
478  //  S = simplify(S,2); // destroys the leadcoef
479  // rescan S and remove 0s from it
480  int i;
481  ideal T;
482  int ss = ncols(S);
483  int j=1;
484  for(i=1; i<=ss; i++)
485  {
486    if (S[i]!=0)
487    {
488      T[j]=S[i];
489      j++;
490    }
491  }
492  S  = T;
493  ss = j-1;
494  int GotF = 1;
495  list F;
496  if (size(#)>0)
497  {
498    F = #;
499    if ( (size(F)==1) && (F[1]==0) )
500    {
501      GotF = 0;
502    }
503  }
504  else
505  {
506    GotF = 0;
507  }
508  int sf = size(F);
509
510  ideal SC;
511  int  sp = size(P);
512  intvec np;
513  int GotP = 1;
514  if (sp==0)
515  {
516    GotP = 0;
517  }
518  if (sp==1)
519  {
520    if (P[1]==0)
521    {
522      GotP = 0;
523    }
524  }
525  if (GotP)
526  {
527    for (i=1; i<=sp; i++)
528    {
529      np[i] = npar(P[i])+ nvars(basering);
530    }
531  }
532  for (i=1; i<=ss; i++)
533  {
534    SC[i] = leadcoef(S[i]);
535  }
536  if (GotF)
537  {
538    for (i=1; i<=sf; i++)
539    {
540      SC[ss+i] = leadcoef(F[i]);
541    }
542  }
543  def br = basering;
544// map all pars except to vars, provided no things are in denominator
545  list M = ringlist(basering);
546  list L = M[1..4]; // erase nc part
547  list Pars = L[1][2];
548  list Vars = L[2] + Pars;
549  L[1] = L[1][1]; // characteristic
550  L[2] = Vars;
551
552  def @R = ring(L);
553  setring @R;
554  ideal SC = imap(br,SC);
555  if (GotP)
556  {
557    for (i=1; i<=sp; i++)
558    {
559      SC = subst(SC,var(np[i]),1);
560    }
561  }
562  poly q=1;
563  q = lcm(q,SC);
564  setring br;
565  poly  q = imap(@R,q);
566  number lq = leadcoef(q);
567  //  lq;
568  number tmp;
569  string sout="";
570  string vname = "u";
571  for (i=1; i<=ss; i++)
572  {
573    tmp  = leadcoef(S[i]);
574    S[i] = S[i]/tmp;
575    tmp = tmp/lq;
576    sout = sout +"+ "+texcoef(tmp)+"\\cdot ("+exp2pt(S[i])+")";
577  }
578  if (GotF)
579  {
580    vname = "p"; //"f";
581    for (i=1; i<=sf; i++)
582    {
583      tmp  = leadcoef(F[i]);
584      F[i] = F[i]/tmp;
585      tmp = tmp/lq;
586      sout = sout +"+ "+texcoef(tmp)+"\\cdot ("+exp2pt(F[i],vname)+")";
587    }
588  }
589  sout = sout[3..size(sout)]; //rm first +
590  return(sout);
591}
592example
593{
594  "EXAMPLE:"; echo=2;
595  ring r = (0,dh,dt,V),(Tx,Tt),dp;
596  poly M = (4*dh*Tx+dt)^2*(Tt-1) + V*Tt*Tx;
597  ideal I = decoef(M,dt);
598  list L; L[1] = V;
599  difpoly2tex(I,L);
600  poly G = V*dh^2*(Tt-Tx)^2;
601  difpoly2tex(I,L,G);
602}
603
604
605proc magnitude(poly P)
606"USAGE:  magnitude(P);  P a poly
607RETURN:  poly
608PURPOSE: compute the square of the magnitude of a complex expression
609ASSUME: i is the variable of a basering
610EXAMPLE: example magnitude; shows examples
611"
612{
613  // check whether i is present among the vars TODO
614  // assume i^2+1=0;
615  poly re = subst(P,i,0);
616  poly im = (P - re)/i;
617  return(re^2+im^2);
618}
619example
620{
621  "EXAMPLE:"; echo=2;
622  ring r = (0,d),(g,i,sin,cos),dp;
623  poly P = d*i*sin - g*cos +d^2;
624  magnitude(P);
625}
626
627
628static proc clTex(string s)
629// removes beginning and ending $'s
630{
631  string t;
632  if (size(s)>2)
633  {
634    // why -3?
635     t = s[2..(size(s)-3)];
636  }
637  return(t);
638}
639
640static proc ElimModComp(module M, int num)
641{
642  // experimental, todo
643  // check whether a curr.ord is elim-comp
644  string s = ordstr(basering);
645  int ns = size(s);
646  //  s[ns] = "c", "C";
647  if (s[1]="c")
648  {
649    // num last components
650  }
651  if (s[1]="C")
652  {
653    // num first components
654  }
655}
656
657
658static proc simfrac(poly up, poly down)
659{
660  // simplifies a fraction up/down
661  // into the form up/down = RT[1] + RT[2]/down
662  list LL = division(up,down);
663  list RT;
664  RT[1] =  LL[1][1,1]; // integer part
665  RT[2] = L[2][1]; // new numerator
666  return(RT);
667}
668
669proc findifs_example()
670"USAGE:  findifs_example();
671RETURN:  nothing (demo)
672PURPOSE:  demonstration of our approach and this library
673EXAMPLE: example findifs_example; shows examples
674"
675{
676
677  "* Equation: u_tt - A^2 u_xx -B^2 u_yy = 0; A,B are constants";
678  "* we employ three central differences";
679  "* the vector we act on is (u_xx, u_yy, u_tt, u)^T";
680  "* Set up the ring: ";
681  "ring r = (0,A,B,dt,dx,dy),(Tx,Ty,Tt),(c,dp);";
682  ring r = (0,A,B,dt,dx,dy),(Tx,Ty,Tt),(c,dp);
683  "* Set up the matrix with equation and approximations: ";
684  "matrix M[4][4] =";
685    "    // direct equation:";
686  "    -A^2, -B^2, 1, 0,";
687  "    // central difference u_tt";
688  "    0, 0,  -dt^2*Tt, (Tt-1)^2,";
689  "    // central difference u_xx";
690  "    -dx^2*Tx, 0, 0, (Tx-1)^2,";
691  "    // central difference u_yy";
692  "    0, -dy^2*Ty, 0, (Ty-1)^2;";
693  matrix M[4][4] =
694    // direct equation:
695    -A^2, -B^2, 1, 0,
696    // central difference u_tt
697    0, 0,  -dt^2*Tt, (Tt-1)^2,
698    // central difference u_xx
699    -dx^2*Tx, 0, 0, (Tx-1)^2,
700    // central difference u_yy
701    0, -dy^2*Ty, 0, (Ty-1)^2;
702//=========================================
703// CHECK THE CORRECTNESS OF EQUATIONS AS INPUT:
704ring rX = (0,A,B,dt,dx,dy,Tx,Ty,Tt),(Uxx, Uyy,Utt, U),(c,Dp);
705matrix M = imap(r,M);
706vector X = [Uxx, Uyy, Utt, U];
707 "* Print the differential form of equations: ";
708print(M*X);
709// END CHECK
710//=========================================
711 setring r;
712 "* Perform the elimination of module components:";
713 " module R = transpose(M);";
714 " module S = std(R);";
715 " poly   p = S[4,1];" ;
716 module R = transpose(M);
717 module S = std(R);
718 poly   p = S[4,1]; // by elimination of module components
719 list L; L[1]=A;L[2] = B;
720 ideal I = decoef(p,dt); // make splitting wrt the appearance of dt
721 "* Create the nodal of the scheme in TeX format: ";
722 " ideal I = decoef(p,dt);";
723 " difpoly2tex(I,L);";
724 difpoly2tex(I,L); // the nodal form of the scheme in TeX
725 "* Preparations for the semi-factorized form: ";
726 poly pi1 = subst(I[2],B,0);
727 poly pi2 = I[2] - pi1;
728 " poly pi1 = subst(I[2],B,0);";
729 " poly pi2 = I[2] - pi1;";
730 "* Show the semi-factorized form of the scheme: 1st summand";
731 " factorize(I[1]); ";
732 factorize(I[1]); // semi-factorized form of the scheme: 1st summand
733 "* Show the semi-factorized form of the scheme: 2nd summand";
734 " factorize(pi1);";
735 factorize(pi1); // semi-factorized form of the scheme: 2nd summand
736 "* Show the semi-factorized form of the scheme: 3rd summand";
737 " factorize(pi1);";
738 factorize(pi2); // semi-factorized form of the scheme: 3rd summand
739}
740example
741{
742  "EXAMPLE:"; echo=1;
743 findifs_example();
744}
Note: See TracBrowser for help on using the repository browser.