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