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