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
