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