reduce-historical/r35/xlog/taylor.log

Codemist Standard Lisp 3.54 for DEC Alpha: May 23 1994
Dump file created: Mon May 23 10:39:11 1994
REDUCE 3.5, 15-Oct-93 ...
Memory allocation: 6023424 bytes

+++ About to read file ndotest.red


comment
        Test and demonstration file for the Taylor expansion package,
        by Rainer M. Schoepf.  Works with version 1.8b (03-Sep-93);


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 A more complicated example;


hugo := taylor(log(1/(1-x)),x,0,5);


             1   2    1   3    1   4    1   5      6
hugo := x + ---*x  + ---*x  + ---*x  + ---*x  + O(x )
             2        3        4        5

taylorcombine(exp(hugo/(1+hugo)));


         1    4      6
1 + x + ----*x  + O(x )
         12


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 The trigonometric functions;


taylor (tan x / x, x, 0, 2);


     1   2      3
1 + ---*x  + O(x )
     3


taylorcombine sin ws;


          cos(1)   2      3
sin(1) + --------*x  + O(x )
            3


taylor (cot x / x, x, 0, 4);


 -2    1     1    2     2    4      5
x   - --- - ----*x  - -----*x  + O(x )
       3     45        945


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 + (4 terms) + O(x ,y )
               2


taylorcombine (ws**2);


             2                                3  3
1 + 2*y + 2*y  + 2*x + 4*y*x + (4 terms) + 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 + (4 terms) + 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 + (4 terms) + 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  + (1 term) + 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  + (1 term) + 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  + (1 term) + O(x )
                          3           3


taylororiginal ws;


 2*x + y
e


taylorcombine (xx1 / x);


   -1              2
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 to finite number of terms;


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    2  3
x - ---*x  + y - ---*y*x  - ---*y *x + ----*y *x  + (2 terms)
     6            2          2          12

      4  4
 + O(x ,y )


taylor (e^x - 1 - x,x,0,6);


 1   2    1   3    1    4     1    5     1    6      7
---*x  + ---*x  + ----*x  + -----*x  + -----*x  + O(x )
 2        6        24        120        720


taylorcombine sqrt ws;


    1              1       2        1        3         1        4
---------*x + -----------*x  + ------------*x  + -------------*x
 sqrt(2)       6*sqrt(2)        36*sqrt(2)        270*sqrt(2)

         1         5      6
 + --------------*x  + O(x )
    2592*sqrt(2)



comment A more complicated example contributed by Stan Kameny;


zz2 := (z*(z-2*pi*i)*(z-pi*i/2)^2)/(sinh z-i);


                 3            2       2        3
        z*(2*i*pi  - 12*i*pi*z  - 9*pi *z + 4*z )
zz2 := -------------------------------------------
                     4*(sinh(z) - i)

dz2 := df(zz2,z);


                         3                      3               2  2
dz2 := ( - 2*cosh(z)*i*pi *z + 12*cosh(z)*i*pi*z  + 9*cosh(z)*pi *z

                      4                 3                    2
         - 4*cosh(z)*z  + 2*sinh(z)*i*pi  - 36*sinh(z)*i*pi*z

                        2                 3          2           3
         - 18*sinh(z)*pi *z + 16*sinh(z)*z  + 18*i*pi *z - 16*i*z

               3          2             2
         + 2*pi  - 36*pi*z )/(4*(sinh(z)  - 2*sinh(z)*i - 1))

z0 := pi*i/2;


       i*pi
z0 := ------
        2

taylor(dz2,z,z0,6);


                2
            - pi  + 16        i*pi      pi        i*pi  2
 - 2*pi + -------------*(z - ------) + ----*(z - ------)
               4*i             2        2          2

        2
    3*pi  - 80        i*pi  3    pi        i*pi  4
 + ------------*(z - ------)  - ----*(z - ------)
      120*i            2         24         2

           2
     - 5*pi  + 168        i*pi  5                      i*pi  7
 + ----------------*(z - ------)  + (1 term) + O((z - ------) )
        3360*i             2                            2



comment A problem are non-analytic terms: there are no precautions
        taken to detect or handle them;


taylor (sqrt (x), x, 0, 2);


taylor(sqrt(x),x,0,2)


taylor (e**(1/x), x, 0, 2);


        1/x
taylor(e   ,x,0,2)


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);


          1
taylor(--------,x,infinity,5)
        sin(x)


clear z;



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 - ---*x  - ---*y  + ---*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 A simple example of Taylor kernel differentiation;


hugo := taylor(e^x,x,0,5);


                 1   2    1   3    1    4     1    5      6
hugo := 1 + x + ---*x  + ---*x  + ----*x  + -----*x  + O(x )
                 2        6        24        120

df(hugo^2,x);


             2    8   3    4   4      5
2 + 4*x + 4*x  + ---*x  + ---*x  + O(x )
                  3        3


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      7
---*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 The following example calculates the perturbation expansion
        of the root x = 20 of the following polynomial in terms of
        EPS, in ROUNDED mode;


poly := for r := 1 : 20 product (x - r);


         20        19          18            17             16
poly := x   - 210*x   + 20615*x   - 1256850*x   + 53327946*x

                       15                14                 13
         - 1672280820*x   + 40171771630*x   - 756111184500*x

                           12                    11
         + 11310276995381*x   - 135585182899530*x

                             10                      9
         + 1307535010540395*x   - 10142299865511450*x

                              8                       7
         + 63030812099294896*x  - 311333643161390640*x

                                6                        5
         + 1206647803780373360*x  - 3599979517947607200*x

                                4                         3
         + 8037811822645051776*x  - 12870931245150988800*x

                                 2
         + 13803759753640704000*x  - 8752948036761600000*x

         + 2432902008176640000


on rounded;



tpoly := taylor (poly, x, 20, 4);


                                                                2
tpoly := 1.21649393692e+17*(x - 20) + 4.31564847287e+17*(x - 20)

                                      3                             4
          + 6.68609351672e+17*(x - 20)  + 6.10115975015e+17*(x - 20)

                      5
          + O((x - 20) )


taylorrevert (tpoly, x, eps);


                                                  2
20 + 8.22034512177e-18*eps - 2.39726594661e-34*eps

                        3                       4        5
 + 1.09290580231e-50*eps  - 5.9711415946e-67*eps  + O(eps )


comment Some more examples using rounded mode;


taylor(sin x/x,x,0,4);


                    2                     4      5
1 - 0.166666666667*x  + 0.00833333333333*x  + O(x )


taylor(sin x,x,pi/2,4);


                                                                   2
1 + 6.12323399574e-17*(x - 1.57079632679) - 0.5*(x - 1.57079632679)

                                        3
 - 1.02053899929e-17*(x - 1.57079632679)

                                      4                        5
 + 0.0416666666667*(x - 1.57079632679)  + O((x - 1.57079632679) )


taylor(tan x,x,pi/2,4);


                      -1
 - (x - 1.57079632679)   + 0.333333333333*(x - 1.57079632679)

                                      3                        5
 + 0.0222222222222*(x - 1.57079632679)  + O((x - 1.57079632679) )


off rounded;




comment An example that involves computing limits of type 0/0 if
        expansion is done via differentiation;



taylor(sqrt((e^x - 1)/x),x,0,15);


     1       5    2     1    3     79     4      3     5
1 + ---*x + ----*x  + -----*x  + -------*x  + -------*x  + (10 terms)
     4       96        128        92160        40960

      16
 + O(x  )



comment Examples involving non-analytical terms;


taylor(log(e^x-1),x,0,5);


           1       1    2     1     4      5
log(x) + (---*x + ----*x  - ------*x  + O(x ))
           2       24        2880


taylor(e^(1/x)*(e^x-1),x,0,5);


 1/x       1   2    1   3    1    4     1    5      6
e   *(x + ---*x  + ---*x  + ----*x  + -----*x  + O(x ))
           2        6        24        120


taylor(log(x)*x^10,x,0,5);


         10      15
log(x)*(x   + O(x  ))


taylor(log(x)*x^10,x,0,11);


         10      21
log(x)*(x   + O(x  ))


taylor(log(x-a)/((a-b)*(a-c)) + log(2(x-b))/((b-c)*(b-a))
     + log(x-c)/((c-a)*(c-b)),x,infinity,2);


           log(2)            1   1        1
 - ---------------------- - ---*---- + O(----)
                 2           2    2        3
    a*b - a*c - b  + b*c         x        x


ss := (sqrt(x^(2/5) +1) - x^(1/3)-1)/x^(1/3);


             2/5         1/3
       sqrt(x    + 1) - x    - 1
ss := ---------------------------
                  1/3
                 x

taylor(exp ss,x,0,2);


            2/5           3
      sqrt(x    + (1 + O(x )))
 exp(--------------------------)
                 1/3
                x
---------------------------------
              3
       1 + O(x )           3
  exp(-----------)*(e + O(x ))
          1/3
         x


taylor(exp sub(x=x^15,ss),x,0,2);


 1      1         1    2      3
--- + -----*x + -----*x  + O(x )
 e     2*e       8*e



comment In the following we demonstrate the possibiblity to compute the
        expansion of a function which is given by a simple first order
        differential equation: the function myexp(x) is exp(-x^2);


operator myexp,myerf;


let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1,
     df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0};


taylor(myexp(x),x,0,5);


     2    1   4      6
1 - x  + ---*x  + O(x )
          2

taylor(myerf(x),x,0,5);


 2*sqrt(pi)       2*sqrt(pi)   3    sqrt(pi)   5      6
------------*x - ------------*x  + ----------*x  + O(x )
     pi              3*pi             5*pi

clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1,
       df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0};


clear myexp,erf;



showtime;


Time: 13350 ms  plus GC time: 534 ms


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.

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;


operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy;


operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy;


operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,
         gxxxx,gxxxy,gxxyy,gxyyy,gyyyy;


operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy;


operator gxxxxyy,gxxxyyy,gxxyyyy;



operator_diff_rules := {
  df(a(~x,~y),~x) => ax(x,y),
  df(a(~x,~y),~y) => ay(x,y),
  df(f(~x,~y),~x) => fx(x,y),
  df(f(~x,~y),~y) => fy(x,y),
  df(gg(~x,~y),~x) => gx(x,y),
  df(gg(~x,~y),~y) => gy(x,y),

  df(ax(~x,~y),~x) => axx(x,y),
  df(ax(~x,~y),~y) => axy(x,y),
  df(ay(~x,~y),~x) => axy(x,y),
  df(ay(~x,~y),~y) => ayy(x,y),
  df(fx(~x,~y),~x) => fxx(x,y),
  df(fx(~x,~y),~y) => fxy(x,y),
  df(fy(~x,~y),~x) => fxy(x,y),
  df(fy(~x,~y),~y) => fyy(x,y),
  df(gx(~x,~y),~x) => gxx(x,y),
  df(gx(~x,~y),~y) => gxy(x,y),
  df(gy(~x,~y),~x) => gxy(x,y),
  df(gy(~x,~y),~y) => gyy(x,y),

  df(axx(~x,~y),~x) => axxx(x,y),
  df(axy(~x,~y),~x) => axxy(x,y),
  df(ayy(~x,~y),~x) => axyy(x,y),
  df(ayy(~x,~y),~y) => ayyy(x,y),
  df(fxx(~x,~y),~x) => fxxx(x,y),
  df(fxy(~x,~y),~x) => fxxy(x,y),
  df(fxy(~x,~y),~y) => fxyy(x,y),
  df(fyy(~x,~y),~x) => fxyy(x,y),
  df(fyy(~x,~y),~y) => fyyy(x,y),
  df(gxx(~x,~y),~x) => gxxx(x,y),
  df(gxx(~x,~y),~y) => gxxy(x,y),
  df(gxy(~x,~y),~x) => gxxy(x,y),
  df(gxy(~x,~y),~y) => gxyy(x,y),
  df(gyy(~x,~y),~x) => gxyy(x,y),
  df(gyy(~x,~y),~y) => gyyy(x,y),

  df(axyy(~x,~y),~x) => axxyy(x,y),
  df(axxy(~x,~y),~x) => axxxy(x,y),
  df(ayyy(~x,~y),~x) => axyyy(x,y),
  df(fxxy(~x,~y),~x) => fxxxy(x,y),
  df(fxyy(~x,~y),~x) => fxxyy(x,y),
  df(fyyy(~x,~y),~x) => fxyyy(x,y),
  df(gxxx(~x,~y),~x) => gxxxx(x,y),
  df(gxxy(~x,~y),~x) => gxxxy(x,y),
  df(gxyy(~x,~y),~x) => gxxyy(x,y),
  df(gyyy(~x,~y),~x) => gxyyy(x,y),
  df(gyyy(~x,~y),~y) => gyyyy(x,y),

  df(axxyy(~x,~y),~x) => axxxyy(x,y),
  df(axyyy(~x,~y),~x) => axxyyy(x,y),
  df(fxxyy(~x,~y),~x) => fxxxyy(x,y),
  df(fxyyy(~x,~y),~x) => fxxyyy(x,y),
  df(gxxxy(~x,~y),~x) => gxxxxy(x,y),
  df(gxxyy(~x,~y),~x) => gxxxyy(x,y),
  df(gxyyy(~x,~y),~x) => gxxyyy(x,y),
  df(gyyyy(~x,~y),~x) => gxyyyy(x,y),

  df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y),
  df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y),
  df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)
};


operator_diff_rules := {df(a(~x,~y),~x) => ax(x,y),

   df(a(~x,~y),~y) => ay(x,y),

   df(f(~x,~y),~x) => fx(x,y),

   df(f(~x,~y),~y) => fy(x,y),

   df(gg(~x,~y),~x) => gx(x,y),

   df(gg(~x,~y),~y) => gy(x,y),

   df(ax(~x,~y),~x) => axx(x,y),

   df(ax(~x,~y),~y) => axy(x,y),

   df(ay(~x,~y),~x) => axy(x,y),

   df(ay(~x,~y),~y) => ayy(x,y),

   df(fx(~x,~y),~x) => fxx(x,y),

   df(fx(~x,~y),~y) => fxy(x,y),

   df(fy(~x,~y),~x) => fxy(x,y),

   df(fy(~x,~y),~y) => fyy(x,y),

   df(gx(~x,~y),~x) => gxx(x,y),

   df(gx(~x,~y),~y) => gxy(x,y),

   df(gy(~x,~y),~x) => gxy(x,y),

   df(gy(~x,~y),~y) => gyy(x,y),

   df(axx(~x,~y),~x) => axxx(x,y),

   df(axy(~x,~y),~x) => axxy(x,y),

   df(ayy(~x,~y),~x) => axyy(x,y),

   df(ayy(~x,~y),~y) => ayyy(x,y),

   df(fxx(~x,~y),~x) => fxxx(x,y),

   df(fxy(~x,~y),~x) => fxxy(x,y),

   df(fxy(~x,~y),~y) => fxyy(x,y),

   df(fyy(~x,~y),~x) => fxyy(x,y),

   df(fyy(~x,~y),~y) => fyyy(x,y),

   df(gxx(~x,~y),~x) => gxxx(x,y),

   df(gxx(~x,~y),~y) => gxxy(x,y),

   df(gxy(~x,~y),~x) => gxxy(x,y),

   df(gxy(~x,~y),~y) => gxyy(x,y),

   df(gyy(~x,~y),~x) => gxyy(x,y),

   df(gyy(~x,~y),~y) => gyyy(x,y),

   df(axyy(~x,~y),~x) => axxyy(x,y),

   df(axxy(~x,~y),~x) => axxxy(x,y),

   df(ayyy(~x,~y),~x) => axyyy(x,y),

   df(fxxy(~x,~y),~x) => fxxxy(x,y),

   df(fxyy(~x,~y),~x) => fxxyy(x,y),

   df(fyyy(~x,~y),~x) => fxyyy(x,y),

   df(gxxx(~x,~y),~x) => gxxxx(x,y),

   df(gxxy(~x,~y),~x) => gxxxy(x,y),

   df(gxyy(~x,~y),~x) => gxxyy(x,y),

   df(gyyy(~x,~y),~x) => gxyyy(x,y),

   df(gyyy(~x,~y),~y) => gyyyy(x,y),

   df(axxyy(~x,~y),~x) => axxxyy(x,y),

   df(axyyy(~x,~y),~x) => axxyyy(x,y),

   df(fxxyy(~x,~y),~x) => fxxxyy(x,y),

   df(fxyyy(~x,~y),~x) => fxxyyy(x,y),

   df(gxxxy(~x,~y),~x) => gxxxxy(x,y),

   df(gxxyy(~x,~y),~x) => gxxxyy(x,y),

   df(gxyyy(~x,~y),~x) => gxxyyy(x,y),

   df(gyyyy(~x,~y),~x) => gxyyyy(x,y),

   df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y),

   df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y),

   df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)}


let operator_diff_rules;



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)*gxy(x,y)*gy(x,y)

         + 2*a(x,y)*fxy(x,y)*gx(x,y)*gy(x,y)

         + a(x,y)*fy(x,y)*gx(x,y)*gxy(x,y)

         + a(x,y)*fy(x,y)*gxx(x,y)*gy(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)*gxy(x,y)*gy(x,y)

           + 2*a(x,y)*fxy(x,y)*gx(x,y)*gy(x,y)

           + a(x,y)*fy(x,y)*gx(x,y)*gxy(x,y)

           + a(x,y)*fy(x,y)*gxx(x,y)*gy(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



clear diff(~f,~arg);


clearrules operator_diff_rules;


clear diff,a,f,gg;


clear ax,ay,fx,fy,gx,gy;


clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy;


clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy;


clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy;


clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy;


clear gxxxxyy,gxxxyyy,gxxyyyy;



taylorprintterms := 5;


taylorprintterms := 5


off taylorautoexpand,taylorkeeporiginal;



comment That's all, folks;


showtime;


Time: 11283 ms  plus GC time: 300 ms

end;
(TIME:  taylor 24633 25467)


End of Lisp run after 24.64+1.58 seconds