Artifact 4553d568adbec47cfebf88e697e8ad84cef5e5327da753263de3a28658a83d51:
- File
r34.1/xlog/taylor.log
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 30542) [annotate] [blame] [check-ins using] [more...]
Sat May 30 16:25:21 PDT 1992 REDUCE 3.4.1, 15-Jul-92 ... 1: 1: 2: 2: (TAYLOR) 3: 3: Time: 187 ms 4: 4: comment Test and demonstration file for the Taylor expansion package, by Rainer M. Schoepf. Works with version 1.4b (16-Apr-92); showtime; Time: 0 ms on errcont; % disable interruption on errors comment Simple Taylor expansion; xx := taylor (e**x, x, 0, 4); 1 2 1 3 1 4 5 XX := 1 + X + ---*X + ---*X + ----*X + O(X ) 2 6 24 yy := taylor (e**y, y, 0, 4); 1 2 1 3 1 4 5 YY := 1 + Y + ---*Y + ---*Y + ----*Y + O(Y ) 2 6 24 comment Basic operations, i.e. addition, subtraction, multiplication, and division are possible: this is not done automatically if the switch TAYLORAUTOCOMBINE is OFF. In this case it is necessary to use taylorcombine; taylorcombine (xx**2); 2 4 3 2 4 5 1 + 2*X + 2*X + ---*X + ---*X + O(X ) 3 3 taylorcombine (ws - xx); 3 2 7 3 5 4 5 X + ---*X + ---*X + ---*X + O(X ) 2 6 8 comment The result is again a Taylor kernel; if taylorseriesp ws then write "OK"; OK comment It is not possible to combine Taylor kernels that were expanded with respect to different variables; taylorcombine (xx**yy); 1 2 1 3 1 4 5 (1 + X + ---*X + ---*X + ----*X + O(X )) 2 6 24 1 2 1 3 1 4 5 **(1 + Y + ---*Y + ---*Y + ----*Y + O(Y )) 2 6 24 comment But we can take the exponential or the logarithm of a Taylor kernel; taylorcombine (e**xx); 2 5*E 3 5*E 4 5 E + E*X + E*X + -----*X + -----*X + O(X ) 6 8 taylorcombine log ws; 1 2 1 3 1 4 5 1 + X + ---*X + ---*X + ----*X + O(X ) 2 6 24 comment We may try to expand about another point; taylor (xx, x, 1, 2); 65 8 5 2 3 ---- + ---*(X - 1) + ---*(X - 1) + O((X - 1) ) 24 3 4 comment Arc tangent is one of the functions this package knows of; xxa := taylorcombine atan ws; XXA := 65 1536 2933040 2 3 ATAN(----) + ------*(X - 1) - ----------*(X - 1) + O((X - 1) ) 24 4801 23049601 comment Sine another one; taylor (tan x / x, x, 0, 2); 1 2 3 1 + ---*X + O(X ) 3 taylorcombine sin ws; TAN(1) 1 2 3 ------------------- + ---------------------*X + O(X ) 2 2 SQRT(TAN(1) + 1) 3*SQRT(TAN(1) + 1) comment Expansion with respect to more than one kernel is possible; xy := taylor (e**(x+y), x, 0, 2, y, 0, 2); 1 2 3 3 XY := 1 + Y + ---*Y + X + Y*X + ... + O(X ,Y ) 2 taylorcombine (ws**2); 2 3 3 1 + 2*Y + 2*Y + 2*X + 4*Y*X + ... + O(X ,Y ) comment We take the inverse and convert back to REDUCE's standard representation; taylorcombine (1/ws); 2 3 3 1 - 2*X + 2*X - 2*Y + 4*Y*X + ... + O(X ,Y ) taylortostandard ws; 2 2 2 2 2 2 4*X *Y - 4*X *Y + 2*X - 4*X*Y + 4*X*Y - 2*X + 2*Y - 2*Y + 1 comment Some examples of Taylor kernel divsion; xx1 := taylor (sin (x), x, 0, 4); 1 3 5 XX1 := X - ---*X + O(X ) 6 taylorcombine (xx/xx1); -1 2 1 2 3 X + 1 + ---*X + ---*X + O(X ) 3 3 taylorcombine (1/xx1); -1 1 3 X + ---*X + O(X ) 6 tt1 := taylor (exp (x), x, 0, 3); 1 2 1 3 4 TT1 := 1 + X + ---*X + ---*X + O(X ) 2 6 tt2 := taylor (sin (x), x, 0, 3); 1 3 4 TT2 := X - ---*X + O(X ) 6 tt3 := taylor (1 + tt2, x, 0, 3); 1 3 4 TT3 := 1 + X - ---*X + O(X ) 6 taylorcombine(tt1/tt2); -1 2 2 X + 1 + ---*X + O(X ) 3 taylorcombine(tt1/tt3); 1 2 1 3 4 1 + ---*X - ---*X + O(X ) 2 6 taylorcombine(tt2/tt1); 2 1 3 4 X - X + ---*X + O(X ) 3 taylorcombine(tt3/tt1); 1 2 1 3 4 1 - ---*X + ---*X + O(X ) 2 6 comment Here's what I call homogeneous expansion; xx := taylor (e**(x*y), {x,y}, 0, 2); 3 XX := 1 + Y*X + O({X,Y} ) xx1 := taylor (sin (x+y), {x,y}, 0, 2); 3 XX1 := Y + X + O({X,Y} ) xx2 := taylor (cos (x+y), {x,y}, 0, 2); 1 2 1 2 3 XX2 := 1 - ---*Y - Y*X - ---*X + O({X,Y} ) 2 2 temp := taylorcombine (xx/xx2); 1 2 1 2 3 TEMP := 1 + ---*Y + 2*Y*X + ---*X + O({X,Y} ) 2 2 taylorcombine (ws*xx2); 3 1 + Y*X + O({X,Y} ) comment The following shows a principal difficulty: since xx1 is symmetric in x and y but has no constant term it is impossible to compute 1/xx1; taylorcombine (1/xx1); ***** Not a unit in argument to INVTAYLOR comment Substitution in Taylor expressions is possible; sub (x=z, xy); 1 2 3 3 1 + Y + ---*Y + Z + Y*Z + ... + O(Z ,Y ) 2 comment Expression dependency in substitution is detected; sub (x=y, xy); ***** Substitution of dependent variables Y Y comment It is possible to replace a Taylor variable by a constant; sub (x=4, xy); 13 2 3 13 + 13*Y + ----*Y + O(Y ) 2 sub (x=4, xx1); 3 4 + Y + O(Y ) comment This package has three switches: TAYLORKEEPORIGINAL, TAYLORAUTOEXPAND, and TAYLORAUTOCOMBINE; on taylorkeeporiginal; temp := taylor (e**(x+y), x, 0, 5); Y Y Y Y Y E 2 E 3 E 4 6 TEMP := E + E *X + ----*X + ----*X + ----*X + ... + O(X ) 2 6 24 taylorcombine (log (temp)); 6 Y + X + O(X ) taylororiginal ws; X + Y taylorcombine (temp * e**x); Y Y Y X Y Y E 2 E 3 E 4 6 E *(E + E *X + ----*X + ----*X + ----*X + ... + O(X )) 2 6 24 on taylorautoexpand; taylorcombine ws; Y Y Y Y Y 2 4*E 3 2*E 4 6 E + 2*E *X + 2*E *X + ------*X + ------*X + ... + O(X ) 3 3 taylororiginal ws; 2*X + Y E taylorcombine (xx1 / x); -1 -1 2 X + Y*X + 1 + O({X,Y} ) on taylorautocombine; xx / xx2; 1 2 1 2 3 1 + ---*Y + 2*Y*X + ---*X + O({X,Y} ) 2 2 ws * xx2; 3 1 + Y*X + O({X,Y} ) comment Another example that shows truncation if Taylor kernels of different expansion order are combined; comment First we increase the number of terms to be printed; taylorprintterms := all; TAYLORPRINTTERMS := ALL p := taylor (x**2 + 2, x, 0, 10); 2 11 P := 2 + X + O(X ) p - x**2; 2 11 2 (2 + X + O(X )) - X p - taylor (x**2, x, 0, 5); 6 2 + O(X ) taylor (p - x**2, x, 0, 6); 7 2 + O(X ) off taylorautocombine; taylorcombine(p-x**2); 11 2 + O(X ) taylorcombine(p - taylor(x**2,x,0,5)); 6 2 + O(X ) comment Switch back; taylorprintterms := 6; TAYLORPRINTTERMS := 6 comment Some more examples; taylor ((1 + x)**n, x, 0, 3); 2 N*(N - 1) 2 N*(N - 3*N + 2) 3 4 1 + N*X + -----------*X + ------------------*X + O(X ) 2 6 taylor (e**(-a*t) * (1 + sin(t)), t, 0, 4); 3 2 A*(A - 2) 2 - A + 3*A - 1 3 1 + ( - A + 1)*T + -----------*T + ------------------*T 2 6 3 2 A*(A - 4*A + 4) 4 5 + -------------------*T + O(T ) 24 operator f; taylor (1 + f(t), t, 0, 3); SUB(T=0,DF(F(T),T,2)) 2 (F(0) + 1) + SUB(T=0,DF(F(T),T))*T + -----------------------*T 2 SUB(T=0,DF(F(T),T,3)) 3 4 + -----------------------*T + O(T ) 6 clear f; taylor (sqrt(1 + a*x + sin(x)), x, 0, 3); 2 3 2 A + 1 - A - 2*A - 1 2 3*A + 9*A + 9*A - 1 3 1 + -------*X + -----------------*X + -----------------------*X 2 8 48 4 + O(X ) taylorcombine (ws**2); 1 3 4 1 + (A + 1)*X - ---*X + O(X ) 6 taylor (sqrt(1 + x), x, 0, 5); 1 1 2 1 3 5 4 7 5 6 1 + ---*X - ---*X + ----*X - -----*X + -----*X + O(X ) 2 8 16 128 256 taylor ((cos(x) - sec(x))^3, x, 0, 5); 6 0 + O(X ) taylor ((cos(x) - sec(x))^-3, x, 0, 5); -6 1 -4 11 -2 347 6767 2 15377 4 - X + ---*X + -----*X - ------- - --------*X - ---------*X 2 120 15120 604800 7983360 6 + O(X ) taylor (sqrt(1 - k^2*sin(x)^2), x, 0, 6); 2 2 2 2 4 2 K 2 K *( - 3*K + 4) 4 K *( - 45*K + 60*K - 16) 6 1 - ----*X + ------------------*X + ----------------------------*X 2 24 720 7 + O(X ) taylor (sin(x + y), x, 0, 3, y, 0, 3); 1 3 1 2 1 2 1 3 2 4 4 Y - ---*Y + X - ---*Y *X - ---*Y*X + ----*Y *X + ... + O(X ,Y ) 6 2 2 12 comment A problem are non-analytic terms: there are no precautions taken to detect or handle them; taylor (sqrt (x), x, 0, 2); ***** Zero divisor ***** Error during expansion (possible singularity!) taylor (e**(1/x), x, 0, 2); ***** Zero divisor ***** Error during expansion (possible singularity!) comment Even worse: you can substitute a non analytical kernel; sub (y = sqrt (x), yy); 1 2 1 3 1 4 1 + SQRT(X) + ---*SQRT(X) + ---*SQRT(X) + ----*SQRT(X) 2 6 24 5 + O(SQRT(X) ) comment Expansion about infinity is possible in principle...; taylor (e**(1/x), x, infinity, 5); 1 1 1 1 1 1 1 1 1 1 1 + --- + ---*---- + ---*---- + ----*---- + -----*---- + O(----) X 2 2 6 3 24 4 120 5 6 X X X X X xi := taylor (sin (1/x), x, infinity, 5); 1 1 1 1 1 1 XI := --- - ---*---- + -----*---- + O(----) X 6 3 120 5 6 X X X y1 := taylor(x/(x-1), x, infinity, 3); 1 1 1 1 Y1 := 1 + --- + ---- + ---- + O(----) X 2 3 4 X X X z := df(y1, x); 1 1 1 1 Z := - ---- - 2*---- - 3*---- + O(----) 2 3 4 5 X X X X comment ...but far from being perfect; taylor (1 / sin (x), x, infinity, 5); ***** Zero divisor ***** Error during expansion (possible singularity!) comment The template of a Taylor kernel can be extracted; taylortemplate yy; {{Y,0,4}} taylortemplate xxa; {{X,1,2}} taylortemplate xi; {{X,INFINITY,5}} taylortemplate xy; {{X,0,2},{Y,0,2}} taylortemplate xx1; {{{X,Y},0,2}} comment Here is a slightly less trivial example; exp := (sin (x) * sin (y) / (x * y))**2; 2 2 SIN(X) *SIN(Y) EXP := ----------------- 2 2 X *Y taylor (exp, x, 0, 1, y, 0, 1); 2 2 1 + O(X ,Y ) taylor (exp, x, 0, 2, y, 0, 2); 1 2 1 2 1 2 2 3 3 1 - ---*Y - ---*X + ---*Y *X + O(X ,Y ) 3 3 9 tt := taylor (exp, {x,y}, 0, 2); 1 2 1 2 3 TT := 1 - ---*Y - ---*X + O({X,Y} ) 3 3 comment An example that uses factorization; on factor; ff := y**5 - 1; 4 3 2 FF := (Y + Y + Y + Y + 1)*(Y - 1) zz := sub (y = taylor(e**x, x, 0, 3), ff); 1 2 1 3 4 4 ZZ := ((1 + X + ---*X + ---*X + O(X )) 2 6 1 2 1 3 4 3 + (1 + X + ---*X + ---*X + O(X )) 2 6 1 2 1 3 4 2 + (1 + X + ---*X + ---*X + O(X )) 2 6 1 2 1 3 4 + (1 + X + ---*X + ---*X + O(X )) + 1) 2 6 1 2 1 3 4 *((1 + X + ---*X + ---*X + O(X )) - 1) 2 6 on exp; zz; 1 2 1 3 4 5 (1 + X + ---*X + ---*X + O(X )) - 1 2 6 comment The following shows the (limited) capabilities to integrate Taylor kernels. Only a toplevel Taylor kernel is supported, in all other cases a warning is printed and the Taylor kernels are converted to standard representation; zz := taylor (sin x, x, 0, 5); 1 3 1 5 6 ZZ := X - ---*X + -----*X + O(X ) 6 120 ww := taylor (cos y, y, 0, 5); 1 2 1 4 6 WW := 1 - ---*Y + ----*Y + O(Y ) 2 24 int (zz, x); 1 2 1 4 1 6 6 ---*X - ----*X + -----*X + O(X ) 2 24 720 int (ww, x); X 2 X 4 6 X - ---*Y + ----*Y + O(Y ) 2 24 int (zz + ww, x); *** Converting Taylor kernels to standard representation 5 3 4 2 X*(X - 30*X + 360*X + 30*Y - 360*Y + 720) ----------------------------------------------- 720 comment And here we present Taylor series reversion. We start with the example given by Knuth for the algorithm; taylor (t - t**2, t, 0, 5); 2 6 T - T + O(T ) taylorrevert (ws, t, x); 2 3 4 5 6 X + X + 2*X + 5*X + 14*X + O(X ) tan!-series := taylor (tan x, x, 0, 5); 1 3 2 5 6 TAN-SERIES := X + ---*X + ----*X + O(X ) 3 15 taylorrevert (tan!-series, x, y); 1 3 1 5 6 Y - ---*Y + ---*Y + O(Y ) 3 5 atan!-series:=taylor (atan y, y, 0, 5); 1 3 1 5 6 ATAN-SERIES := Y - ---*Y + ---*Y + O(Y ) 3 5 tmp := taylor (e**x, x, 0, 5); 1 2 1 3 1 4 1 5 6 TMP := 1 + X + ---*X + ---*X + ----*X + -----*X + O(X ) 2 6 24 120 taylorrevert (tmp, x, y); 1 2 1 3 1 4 1 5 Y - 1 - ---*(Y - 1) + ---*(Y - 1) - ---*(Y - 1) + ---*(Y - 1) 2 3 4 5 6 + O((Y - 1) ) taylor (log y, y, 1, 5); 1 2 1 3 1 4 1 5 Y - 1 - ---*(Y - 1) + ---*(Y - 1) - ---*(Y - 1) + ---*(Y - 1) 2 3 4 5 6 + O((Y - 1) ) comment An application is the problem posed by Prof. Stanley: we prove that the finite difference expression below corresponds to the given derivative expression; operator diff,a,f,gg; % We use gg to avoid conflict with high energy % physics operator. for all f,arg let diff(f,arg) = df(f,arg); derivative!_expression := diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) + diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ; DERIVATIVE_EXPRESSION := 2*A(X,Y)*DF(F(X,Y),X,Y)*DF(GG(X,Y),X)*DF(GG(X,Y),Y) + A(X,Y)*DF(F(X,Y),X)*DF(GG(X,Y),X,Y)*DF(GG(X,Y),Y) + A(X,Y)*DF(F(X,Y),X)*DF(GG(X,Y),X)*DF(GG(X,Y),Y,2) + A(X,Y)*DF(F(X,Y),Y)*DF(GG(X,Y),X,Y)*DF(GG(X,Y),X) + A(X,Y)*DF(F(X,Y),Y)*DF(GG(X,Y),X,2)*DF(GG(X,Y),Y) + DF(A(X,Y),X)*DF(F(X,Y),Y)*DF(GG(X,Y),X)*DF(GG(X,Y),Y) + DF(A(X,Y),Y)*DF(F(X,Y),X)*DF(GG(X,Y),X)*DF(GG(X,Y),Y) finite!_difference!_expression := +a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) +a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) +a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) +a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) -gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) -gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) -gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) -a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) -gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) +gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2) -gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) -a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) -a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) -a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) +gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) +a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) +a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) -gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) -a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) -a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) +a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) +f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) +a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) +a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) +a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) +a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) -gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) -gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) -gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) -a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) -gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) +gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2) -gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) -a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) -a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) -a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) +gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) +a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) +a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) -gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) -a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) -a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) +a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) +f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) +f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) +f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) +f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) -f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) -f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) -f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) -f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) +f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2) +f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2) +a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) +a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) +a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) +a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) -gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) -gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) -gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) -a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) -gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) +gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2) -gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) -a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) -a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) -a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) +gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) +a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) +a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) -gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) -a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) -a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) +a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) +f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) +a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) +a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) +a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) +a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) -f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) -gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) -gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) -gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) -a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) -gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) +gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2) -gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) -a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) -a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) -a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) +gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) +a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) +a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) -gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) -a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) -a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) +gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) +a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) +f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2) +f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) +f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) +f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) +f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) -f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) -f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) -f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) -f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) +f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2) +f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2) +f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2) +f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2) +a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2) -f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2) -a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$ comment We define abbreviations for the partial derivatives; operator ax,ay,fx,fy,gx,gy; for all x,y let df(a(x,y),x) = ax(x,y); for all x,y let df(a(x,y),y) = ay(x,y); for all x,y let df(f(x,y),x) = fx(x,y); for all x,y let df(f(x,y),y) = fy(x,y); for all x,y let df(gg(x,y),x) = gx(x,y); for all x,y let df(gg(x,y),y) = gy(x,y); operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; for all x,y let df(ax(x,y),x) = axx(x,y); for all x,y let df(ax(x,y),y) = axy(x,y); for all x,y let df(ay(x,y),x) = axy(x,y); for all x,y let df(ay(x,y),y) = ayy(x,y); for all x,y let df(fx(x,y),x) = fxx(x,y); for all x,y let df(fx(x,y),y) = fxy(x,y); for all x,y let df(fy(x,y),x) = fxy(x,y); for all x,y let df(fy(x,y),y) = fyy(x,y); for all x,y let df(gx(x,y),x) = gxx(x,y); for all x,y let df(gx(x,y),y) = gxy(x,y); for all x,y let df(gy(x,y),x) = gxy(x,y); for all x,y let df(gy(x,y),y) = gyy(x,y); operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; for all x,y let df(axx(x,y),x) = axxx(x,y); for all x,y let df(axy(x,y),x) = axxy(x,y); for all x,y let df(ayy(x,y),x) = axyy(x,y); for all x,y let df(ayy(x,y),y) = ayyy(x,y); for all x,y let df(fxx(x,y),x) = fxxx(x,y); for all x,y let df(fxy(x,y),x) = fxxy(x,y); for all x,y let df(fxy(x,y),y) = fxyy(x,y); for all x,y let df(fyy(x,y),x) = fxyy(x,y); for all x,y let df(fyy(x,y),y) = fyyy(x,y); for all x,y let df(gxx(x,y),x) = gxxx(x,y); for all x,y let df(gxx(x,y),y) = gxxy(x,y); for all x,y let df(gxy(x,y),x) = gxxy(x,y); for all x,y let df(gxy(x,y),y) = gxyy(x,y); for all x,y let df(gyy(x,y),x) = gxyy(x,y); for all x,y let df(gyy(x,y),y) = gyyy(x,y); operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy, gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; for all x,y let df(axyy(x,y),x) = axxyy(x,y); for all x,y let df(axxy(x,y),x) = axxxy(x,y); for all x,y let df(ayyy(x,y),x) = axyyy(x,y); for all x,y let df(fxxy(x,y),x) = fxxxy(x,y); for all x,y let df(fxyy(x,y),x) = fxxyy(x,y); for all x,y let df(fyyy(x,y),x) = fxyyy(x,y); for all x,y let df(gxxx(x,y),x) = gxxxx(x,y); for all x,y let df(gxxy(x,y),x) = gxxxy(x,y); for all x,y let df(gxyy(x,y),x) = gxxyy(x,y); for all x,y let df(gyyy(x,y),x) = gxyyy(x,y); for all x,y let df(gyyy(x,y),y) = gyyyy(x,y); operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; for all x,y let df(axxyy(x,y),x) = axxxyy(x,y); for all x,y let df(axyyy(x,y),x) = axxyyy(x,y); for all x,y let df(fxxyy(x,y),x) = fxxxyy(x,y); for all x,y let df(fxyyy(x,y),x) = fxxyyy(x,y); for all x,y let df(gxxxy(x,y),x) = gxxxxy(x,y); for all x,y let df(gxxyy(x,y),x) = gxxxyy(x,y); for all x,y let df(gxyyy(x,y),x) = gxxyyy(x,y); for all x,y let df(gyyyy(x,y),x) = gxyyyy(x,y); operator gxxxxyy,gxxxyyy,gxxyyyy; for all x,y let df(gxxxyy(x,y),x) = gxxxxyy(x,y); for all x,y let df(gxxyyy(x,y),x) = gxxxyyy(x,y); for all x,y let df(gxyyyy(x,y),x) = gxxyyyy(x,y); texp := taylor (finite!_difference!_expression, dx, 0, 1, dy, 0, 1); TEXP := A(X,Y)*FX(X,Y)*GX(X,Y)*GYY(X,Y) + A(X,Y)*FX(X,Y)*GY(X,Y)*GXY(X,Y) + A(X,Y)*FY(X,Y)*GX(X,Y)*GXY(X,Y) + A(X,Y)*FY(X,Y)*GY(X,Y)*GXX(X,Y) + 2*A(X,Y)*GX(X,Y)*GY(X,Y)*FXY(X,Y) + AX(X,Y)*FY(X,Y)*GX(X,Y)*GY(X,Y) 2 2 + AY(X,Y)*FX(X,Y)*GX(X,Y)*GY(X,Y) + O(DX ,DY ) comment You may also try to expand further but this needs a lot of CPU time. Therefore the following line is commented out; %texp := taylor (finite!_difference!_expression, dx, 0, 2, dy, 0, 2); factor dx,dy; result := taylortostandard texp; RESULT := A(X,Y)*FX(X,Y)*GX(X,Y)*GYY(X,Y) + A(X,Y)*FX(X,Y)*GY(X,Y)*GXY(X,Y) + A(X,Y)*FY(X,Y)*GX(X,Y)*GXY(X,Y) + A(X,Y)*FY(X,Y)*GY(X,Y)*GXX(X,Y) + 2*A(X,Y)*GX(X,Y)*GY(X,Y)*FXY(X,Y) + AX(X,Y)*FY(X,Y)*GX(X,Y)*GY(X,Y) + AY(X,Y)*FX(X,Y)*GX(X,Y)*GY(X,Y) derivative!_expression - result; 0 comment That's all, folks; showtime; Time: 9945 ms plus GC time: 459 ms end; 5: 5: Time: 17 ms 6: 6: Quitting Sat May 30 16:25:34 PDT 1992