source:git/Singular/LIB/findifs.lib@867e1a3

spielwiese
Last change on this file since 867e1a3 was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to 100755
File size: 17.0 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="$Id$";
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 polynomial 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 polynomial 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 polynomial into summands
341  // w.r.t. 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 w.r.t. 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.