Artifact 94c3e296c7919826d7c024af599d23caa070c115009bc43b03a9fa35deaa7444:
- Executable file
r37/packages/taylor/taylor.rlg
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 63428) [annotate] [blame] [check-ins using] [more...]
Mon Jun 7 18:22:24 MET DST 1999 REDUCE 3.7, 15-Apr-1999 ... 1: 1: 2: 2: 2: 2: 2: 2: 2: 2: 2: 3: 3: comment Test and demonstration file for the Taylor expansion package, by Rainer M. Schoepf. Works with version 2.2 (18-Jun-97); %%% showtime; 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 taylorcombine (xx**3); 9 2 9 3 27 4 5 1 + 3*x + ---*x + ---*x + ----*x + O(x ) 2 2 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; 65 1536 2933040 2 3 xxa := 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 The poles of these functions are correctly handled; taylor(tan x,x,pi/2,0); pi -1 pi - (x - ----) + O(x - ----) 2 2 taylor(tan x,x,pi/2,3); pi -1 1 pi 1 pi 3 pi 4 - (x - ----) + ---*(x - ----) + ----*(x - ----) + O((x - ----) ) 2 3 2 45 2 2 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); ***** Invalid substitution in Taylor kernel: 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 ) sub (y=0, ws); 4 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/(1+y^4+x^2*y^2+x^4),{x,y},0,6); 4 2 2 4 7 1 - y - y *x - x + O({x,y} ) 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 taylor(f(sqrt(x^2+y^2)),x,x0,4,y,y0,4); 2 2 2 2 f(sqrt(x0 + y0 )) + sub(y=y0,df(f(sqrt(x0 + y )),y))*(y - y0) 2 2 sub(y=y0,df(f(sqrt(x0 + y )),y,2)) 2 + -------------------------------------*(y - y0) 2 2 2 sub(y=y0,df(f(sqrt(x0 + y )),y,3)) 3 + -------------------------------------*(y - y0) 6 2 2 sub(y=y0,df(f(sqrt(x0 + y )),y,4)) 4 + -------------------------------------*(y - y0) 24 2 2 + sub(x=x0,df(f(sqrt(x + y0 )),x))*(x - x0) + (19 terms) 5 5 + O((x - x0) ,(y - y0) ) 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 4 1 + -------*x + -----------------*x + -----------------------*x + O(x ) 2 8 48 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 6 - x + ---*x + -----*x - ------- - --------*x - ---------*x + O(x ) 2 120 15120 604800 7983360 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 7 1 - ----*x + ------------------*x + ----------------------------*x + O(x ) 2 24 720 taylor (sin(x + y), x, 0, 3, y, 0, 3); 1 3 1 2 1 2 1 2 3 4 4 x - ---*x + y - ---*y*x - ---*y *x + ----*y *x + (2 terms) + O(x ,y ) 6 2 2 12 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) taylor(sin(x)/x,x,1,2); - 2*cos(1) + sin(1) 2 sin(1) + (cos(1) - sin(1))*(x - 1) + ----------------------*(x - 1) 2 3 + O((x - 1) ) taylor((sqrt(4+h)-2)/h,h,0,5); 1 1 1 2 5 3 7 4 21 5 6 --- - ----*h + -----*h - -------*h + --------*h - ---------*h + O(h ) 4 64 512 16384 131072 2097152 taylor((sqrt(x)-2)/(4-x),x,4,2); 1 1 1 2 3 - --- + ----*(x - 4) - -----*(x - 4) + O((x - 4) ) 4 64 512 taylor((sqrt(y+4)-2)/(-y),y,0,2); 1 1 1 2 3 - --- + ----*y - -----*y + O(y ) 4 64 512 taylor(x*tanh(x)/(sqrt(1-x^2)-1),x,0,3); 7 2 4 - 2 + ---*x + O(x ) 6 taylor((e^(5*x)-2*x)^(1/x),x,0,2); 3 3 3 73*e 2 3 e + 8*e *x + -------*x + O(x ) 3 taylor(sin x/cos x,x,pi/2,3); pi -1 1 pi 1 pi 3 pi 4 - (x - ----) + ---*(x - ----) + ----*(x - ----) + O((x - ----) ) 2 3 2 45 2 2 taylor(log x*sin(x^2)/(x*sinh x),x,0,2); 1 2 3 log(x)*(1 - ---*x + O(x )) 6 taylor(1/x-1/sin x,x,0,2); 1 3 - ---*x + O(x ) 6 taylor(tan x/log cos x,x,pi/2,2); pi -1 pi - (x - ----) + O(x - ----) 2 2 ------------------------------- pi - 2*x log(----------) 2 taylor(log(x^2/(x^2-a)),x,0,3); 2 - x taylor(log(--------),x,0,3) 2 a - x comment Three more complicated examples 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 3 - 18*sinh(z)*pi *z + 16*sinh(z)*z + 18*i*pi *z - 16*i*z + 2*pi 2 2 - 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 i*(pi - 16) i*pi pi i*pi 2 - 2*pi + --------------*(z - ------) + ----*(z - ------) 4 2 2 2 2 i*( - 3*pi + 80) i*pi 3 pi i*pi 4 + -------------------*(z - ------) - ----*(z - ------) 120 2 24 2 2 i*(5*pi - 168) i*pi 5 i*pi 7 + -----------------*(z - ------) + (1 term) + O((z - ------) ) 3360 2 2 zz3:=(z*(z-2*pi)*(z-pi/2)^2)/(sin z-1); 3 2 2 3 z*( - 2*pi + 9*pi *z - 12*pi*z + 4*z ) zz3 := ------------------------------------------ 4*(sin(z) - 1) dz3 := df(zz3,z); 3 2 2 3 4 dz3 := (2*cos(z)*pi *z - 9*cos(z)*pi *z + 12*cos(z)*pi*z - 4*cos(z)*z 3 2 2 3 - 2*sin(z)*pi + 18*sin(z)*pi *z - 36*sin(z)*pi*z + 16*sin(z)*z 3 2 2 3 2 + 2*pi - 18*pi *z + 36*pi*z - 16*z )/(4*(sin(z) - 2*sin(z) + 1)) z1 := pi/2; pi z1 := ---- 2 taylor(dz3,z,z1,6); 2 2 pi - 16 pi pi pi 2 3*pi - 80 pi 3 2*pi + ----------*(z - ----) + ----*(z - ----) + ------------*(z - ----) 4 2 2 2 120 2 2 pi pi 4 5*pi - 168 pi 5 pi 7 + ----*(z - ----) + -------------*(z - ----) + (1 term) + O((z - ----) ) 24 2 3360 2 2 taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,6); 5 2 1313 4 2773 6 7 1 + ---*x + ------*x - -------*x + O(x ) 3 1890 11907 comment If the expansion point is not constant, it has to be taken care of in differentation, as the following examples show; taylor(sin(x+a),x,a,8); sin(2*a) 2 cos(2*a) 3 sin(2*a) + cos(2*a)*(x - a) - ----------*(x - a) - ----------*(x - a) 2 6 sin(2*a) 4 cos(2*a) 5 9 + ----------*(x - a) + ----------*(x - a) + (3 terms) + O((x - a) ) 24 120 df(ws,a); cos(2*a) 2 sin(2*a) 3 cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a) + ----------*(x - a) 2 6 cos(2*a) 4 sin(2*a) 5 8 + ----------*(x - a) - ----------*(x - a) + (2 terms) + O((x - a) ) 24 120 taylor(cos(x+a),x,a,7); cos(2*a) 2 sin(2*a) 3 cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a) + ----------*(x - a) 2 6 cos(2*a) 4 sin(2*a) 5 8 + ----------*(x - a) - ----------*(x - a) + (2 terms) + O((x - a) ) 24 120 comment A problem are non-analytical terms: rational powers and logarithmic terms can be handled, but other types of essential singularities cannot; taylor(sqrt(x),x,0,2); 1/2 3 x + O(x ) taylor(asinh(1/x),x,0,5); 1 1 2 3 4 5 6 7 log(---) + (log(2) + ---*x - ----*x + ----*x + O(x )) x 4 32 96 taylor(e**(1/x),x,0,2); 1/x taylor(e ,x,0,2) comment Another example for non-integer powers; sub (y = sqrt (x), yy); 1/2 1 1 3/2 1 2 5/2 1 + x + ---*x + ---*x + ----*x + O(x ) 2 6 24 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 You may access the expansion with the PART operator; part(yy,0); plus part(yy,1); 1 part(yy,4); 3 y ---- 6 part(yy,6); ***** Expression taylor(1 + y + 1/2*y**2 + 1/6*y**3 + 1/24*y**4,y,0,4) does not have part 6 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 1 2 1 3 4 3 zz := ((1 + x + ---*x + ---*x + O(x )) + (1 + x + ---*x + ---*x + O(x )) 2 6 2 6 1 2 1 3 4 2 1 2 1 3 4 + (1 + x + ---*x + ---*x + O(x )) + (1 + x + ---*x + ---*x + O(x )) 2 6 2 6 1 2 1 3 4 + 1)*((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 simple cases are supported, otherwise 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); 1 2 1 4 1 6 7 x 2 x 4 6 (---*x - ----*x + -----*x + O(x )) + (x - ---*y + ----*y + O(y )) 2 24 720 2 24 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 6 y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 2 3 4 5 taylor (log y, y, 1, 5); 1 2 1 3 1 4 1 5 6 y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 2 3 4 5 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 15 poly := x - 210*x + 20615*x - 1256850*x + 53327946*x - 1672280820*x 14 13 12 + 40171771630*x - 756111184500*x + 11310276995381*x 11 10 9 - 135585182899530*x + 1307535010540395*x - 10142299865511450*x 8 7 6 + 63030812099294896*x - 311333643161390640*x + 1206647803780373360*x 5 4 - 3599979517947607200*x + 8037811822645051776*x 3 2 - 12870931245150988800*x + 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 3 20 + 8.22034512178e-18*eps - 2.39726594662e-34*eps + 1.09290580232e-50*eps 4 5 - 5.97114159465e-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 4 - 1.02053899929e-17*(x - 1.57079632679) + 0.0416666666667*(x - 1.57079632679) 5 + 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 16 1 + ---*x + ----*x + -----*x + -------*x + -------*x + (10 terms) + O(x ) 4 96 128 92160 40960 comment An example that involves intermediate non-analytical terms which cancel entirely; taylor(x^(5/2)/(log(x+1)*tan(x^(3/2))),x,0,5); 1 1 2 7 3 139 4 67 5 6 1 + ---*x - ----*x - ----*x - -----*x + ------*x + O(x ) 2 12 24 720 1440 comment Other 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 11 log(x)*(x + O(x )) taylor(log(x)*x^10,x,0,11); 10 12 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); 1 1 1/15 1 2/15 1 1/5 1 4/15 1 1/3 --- + -----*x + -----*x + ------*x + -------*x + --------*x e 2*e 8*e 48*e 384*e 3840*e 31/15 + (25 terms) + O(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 taylor(dilog(x),x,0,4); 1 2 1 3 1 4 5 log(x)*(x + ---*x + ---*x + ---*x + O(x )) 2 3 4 2 pi 1 2 1 3 1 4 5 + (----- - x - ---*x - ---*x - ----*x + O(x )) 6 4 9 16 taylor(ei(x),x,0,4); 1 2 1 3 1 4 5 log(x) - psi(1) + (x + ---*x + ----*x + ----*x + O(x )) 4 18 96 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 2 3 1 5 6 ----------*x - ------------*x + ------------*x + O(x ) sqrt(pi) 3*sqrt(pi) 5*sqrt(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; comment There are two special operators, implicit_taylor and inverse_taylor, to compute the Taylor expansion of implicit or inverse functions; implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5); 1 2 1 4 6 1 - ---*x - ---*x + O(x ) 2 8 implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20); 1 2 1 4 1 6 5 8 7 10 21 1 - ---*x - ---*x - ----*x - -----*x - -----*x + (5 terms) + O(x ) 2 8 16 128 256 implicit_taylor(x+y^3-y,x,y,0,0,8); 3 5 7 9 x + x + 3*x + 12*x + O(x ) implicit_taylor(x+y^3-y,x,y,0,1,5); 1 3 2 1 3 105 4 3 5 6 1 - ---*x - ---*x - ---*x - -----*x - ---*x + O(x ) 2 8 2 128 2 implicit_taylor(x+y^3-y,x,y,0,-1,5); 1 3 2 1 3 105 4 3 5 6 - 1 - ---*x + ---*x - ---*x + -----*x - ---*x + O(x ) 2 8 2 128 2 implicit_taylor(y*e^y-x,x,y,0,0,5); 2 3 3 8 4 125 5 6 x - x + ---*x - ---*x + -----*x + O(x ) 2 3 24 comment This is the function exp(-1/x^2), which has an essential singularity at the point 0; implicit_taylor(x^2*log y+1,x,y,0,0,3); ***** Computation of Taylor series of implicit function failed Input expression non-zero at given point inverse_taylor(exp(x)-1,x,y,0,8); 1 2 1 3 1 4 1 5 1 6 9 y - ---*y + ---*y - ---*y + ---*y - ---*y + (2 terms) + O(y ) 2 3 4 5 6 inverse_taylor(exp(x),x,y,0,5); 1 2 1 3 1 4 1 5 6 y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 2 3 4 5 inverse_taylor(sqrt(x),x,y,0,5); 2 6 y + O(y ) inverse_taylor(log(1+x),x,y,0,5); 1 2 1 3 1 4 1 5 6 y + ---*y + ---*y + ----*y + -----*y + O(y ) 2 6 24 120 inverse_taylor((e^x-e^(-x))/2,x,y,0,5); 1 3 3 5 6 y - ---*y + ----*y + O(y ) 6 40 comment In the next two cases the inverse functions have a branch point, therefore the computation fails; inverse_taylor((e^x+e^(-x))/2,x,y,0,5); ***** Computation of Taylor series of inverse function failed inverse_taylor(exp(x^2-1),x,y,0,5); ***** Computation of Taylor series of inverse function failed inverse_taylor(exp(sqrt(x))-1,x,y,0,5); 2 3 11 4 5 5 6 y - y + ----*y - ---*y + O(y ) 12 6 inverse_taylor(x*exp(x),x,y,0,5); 2 3 3 8 4 125 5 6 y - y + ---*y - ---*y + -----*y + O(y ) 2 3 24 %%% showtime; 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; %%% showtime; comment An example provided by Alan Barnes: elliptic functions; % Jacobi's elliptic functions % sn(x,k), cn(x,k), dn(x,k). % The modulus and complementary modulus are denoted by K and K!' % usually written mathematically as k and k' respectively % % epsilon(x,k) is the incomplete elliptic integral of the second kind % usually written mathematically as E(x,k) % % KK(k) is the complete elliptic integral of the first kind % usually written mathematically as K(k) % K(k) = arcsn(1,k) % KK!'(k) is the complementary complete integral % usually written mathematically as K'(k) % NB. K'(k) = K(k') % EE(k) is the complete elliptic integral of the second kind % usually written mathematically as E(k) % EE!'(k) is the complementary complete integral % usually written mathematically as E'(k) % NB. E'(k) = E(k') operator sn, cn, dn, epsilon; elliptic_rules := { % Differentiation rules for basic functions df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), df(cn(~x,~k),~x) => -sn(x,k)*dn(x,k), df(dn(~x,~k),~x) => -k^2*sn(x,k)*cn(x,k), df(epsilon(~x,~k),~x)=> dn(x,k)^2, % k-derivatives % DF Lawden Elliptic Functions & Applications Springer (1989) df(sn(~x,~k),~k) => (k*sn(x,k)*cn(x,k)^2 -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2) + x*cn(x,k)*dn(x,k)/k, df(cn(~x,~k),~k) => (-k*sn(x,k)^2*cn(x,k) +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2) - x*sn(x,k)*dn(x,k)/k, df(dn(~x,~k),~k) => k*(-sn(x,k)^2*dn(x,k) +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2) - k*x*sn(x,k)*cn(x,k), df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k) -epsilon(x,k)*cn(x,k)^2)/(1-k^2) -k*x*sn(x,k)^2, % parity properties sn(-~x,~k) => -sn(x,k), cn(-~x,~k) => cn(x,k), dn(-~x,~k) => dn(x,k), epsilon(-~x,~k) => -epsilon(x,k), sn(~x,-~k) => sn(x,k), cn(~x,-~k) => cn(x,k), dn(~x,-~k) => dn(x,k), epsilon(~x,-~k) => epsilon(x,k), % behaviour at zero sn(0,~k) => 0, cn(0,~k) => 1, dn(0,~k) => 1, epsilon(0,~k) => 0, % degenerate cases of modulus sn(~x,0) => sin(x), cn(~x,0) => cos(x), dn(~x,0) => 1, epsilon(~x,0) => x, sn(~x,1) => tanh(x), cn(~x,1) => 1/cosh(x), dn(~x,1) => 1/cosh(x), epsilon(~x,1) => tanh(x) }; elliptic_rules := {df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), df(cn(~x,~k),~x) => - sn(x,k)*dn(x,k), 2 df(dn(~x,~k),~x) => - k *sn(x,k)*cn(x,k), 2 df(epsilon(~x,~k),~x) => dn(x,k) , 2 dn(x,k) k*sn(x,k)*cn(x,k) - epsilon(x,k)*cn(x,k)*--------- k df(sn(~x,~k),~k) => ----------------------------------------------------- 2 1 - k dn(x,k) + x*cn(x,k)*---------, k 2 dn(x,k) - k*sn(x,k) *cn(x,k) + epsilon(x,k)*sn(x,k)*--------- k df(cn(~x,~k),~k) => -------------------------------------------------------- 2 1 - k dn(x,k) - x*sn(x,k)*---------, k 2 - sn(x,k) *dn(x,k) + epsilon(x,k)*sn(x,k)*cn(x,k) df(dn(~x,~k),~k) => k*---------------------------------------------------- 2 1 - k - k*x*sn(x,k)*cn(x,k), df(epsilon(~x,~k),~k) 2 sn(x,k)*cn(x,k)*dn(x,k) - epsilon(x,k)*cn(x,k) 2 => k*------------------------------------------------- - k*x*sn(x,k) , 2 1 - k sn( - ~x,~k) => - sn(x,k), cn( - ~x,~k) => cn(x,k), dn( - ~x,~k) => dn(x,k), epsilon( - ~x,~k) => - epsilon(x,k), sn(~x, - ~k) => sn(x,k), cn(~x, - ~k) => cn(x,k), dn(~x, - ~k) => dn(x,k), epsilon(~x, - ~k) => epsilon(x,k), sn(0,~k) => 0, cn(0,~k) => 1, dn(0,~k) => 1, epsilon(0,~k) => 0, sn(~x,0) => sin(x), cn(~x,0) => cos(x), dn(~x,0) => 1, epsilon(~x,0) => x, sn(~x,1) => tanh(x), 1 cn(~x,1) => ---------, cosh(x) 1 dn(~x,1) => ---------, cosh(x) epsilon(~x,1) => tanh(x)} let elliptic_rules; hugo := taylor(sn(x,k),k,0,6); 2 2 cos(x)*(cos(x) *x + cos(x)*sin(x) + sin(x) *x - 2*x) 2 hugo := sin(x) + ------------------------------------------------------*k + ( 4 5 4 2 4 cos(x) *x - 2*cos(x) *sin(x)*x + 5*cos(x) *sin(x) 3 2 3 2 3 2 - 10*cos(x) *sin(x) *x + 6*cos(x) *x - 4*cos(x) *sin(x) *x 2 3 2 2 2 + cos(x) *sin(x) + 8*cos(x) *sin(x)*x + 4*cos(x) *sin(x) 4 2 - 11*cos(x)*sin(x) *x + 30*cos(x)*sin(x) *x - 16*cos(x)*x 5 2 3 2 2 4 - 2*sin(x) *x + 8*sin(x) *x - 8*sin(x)*x )/64*k + ( 7 3 7 6 2 - 6*cos(x) *x + 17*cos(x) *x - 99*cos(x) *sin(x)*x 6 5 2 3 5 2 + 21*cos(x) *sin(x) - 18*cos(x) *sin(x) *x - 71*cos(x) *sin(x) *x 5 3 5 4 3 2 + 36*cos(x) *x - 18*cos(x) *x - 135*cos(x) *sin(x) *x 4 3 4 2 4 - 133*cos(x) *sin(x) + 324*cos(x) *sin(x)*x + 172*cos(x) *sin(x) 3 4 3 3 4 - 18*cos(x) *sin(x) *x - 13*cos(x) *sin(x) *x 3 2 3 3 2 3 3 + 72*cos(x) *sin(x) *x - 156*cos(x) *sin(x) *x - 72*cos(x) *x 3 2 5 2 2 5 + 160*cos(x) *x + 27*cos(x) *sin(x) *x - 118*cos(x) *sin(x) 2 3 2 2 2 + 176*cos(x) *sin(x) - 108*cos(x) *sin(x)*x + 32*cos(x) *sin(x) 6 3 6 4 3 - 6*cos(x)*sin(x) *x + 75*cos(x)*sin(x) *x + 36*cos(x)*sin(x) *x 4 2 3 2 - 498*cos(x)*sin(x) *x - 72*cos(x)*sin(x) *x + 888*cos(x)*sin(x) *x 3 7 2 5 2 + 48*cos(x)*x - 384*cos(x)*x + 63*sin(x) *x - 324*sin(x) *x 3 2 2 6 7 + 540*sin(x) *x - 288*sin(x)*x )/2304*k + O(k ) otto := taylor(cn(x,k),k,0,6); 2 2 sin(x)*( - cos(x) *x - cos(x)*sin(x) - sin(x) *x + 2*x) 2 otto := cos(x) + ---------------------------------------------------------*k + 4 5 2 4 3 2 2 ( - 2*cos(x) *x - 5*cos(x) *sin(x)*x - 4*cos(x) *sin(x) *x 3 2 3 2 2 3 - 7*cos(x) *sin(x) + 8*cos(x) *x + 2*cos(x) *sin(x) *x 2 4 2 4 + 2*cos(x) *sin(x)*x - 2*cos(x)*sin(x) *x - 3*cos(x)*sin(x) 2 2 2 2 5 + 8*cos(x)*sin(x) *x - 4*cos(x)*sin(x) - 8*cos(x)*x + 7*sin(x) *x 3 4 7 2 - 22*sin(x) *x + 16*sin(x)*x)/64*k + ( - 9*cos(x) *x 6 3 6 5 2 2 + 6*cos(x) *sin(x)*x - 71*cos(x) *sin(x)*x + 135*cos(x) *sin(x) *x 5 2 5 2 4 3 3 - 66*cos(x) *sin(x) - 36*cos(x) *x + 18*cos(x) *sin(x) *x 4 3 4 3 4 - cos(x) *sin(x) *x - 36*cos(x) *sin(x)*x + 18*cos(x) *sin(x)*x 3 4 2 3 4 + 297*cos(x) *sin(x) *x + 61*cos(x) *sin(x) 3 2 2 3 2 3 2 - 720*cos(x) *sin(x) *x - 208*cos(x) *sin(x) + 252*cos(x) *x 2 5 3 2 5 + 18*cos(x) *sin(x) *x + 31*cos(x) *sin(x) *x 2 3 3 2 3 - 72*cos(x) *sin(x) *x - 24*cos(x) *sin(x) *x 2 3 2 6 2 + 72*cos(x) *sin(x)*x + 56*cos(x) *sin(x)*x + 153*cos(x)*sin(x) *x 6 4 2 4 + 91*cos(x)*sin(x) - 684*cos(x)*sin(x) *x - 212*cos(x)*sin(x) 2 2 2 2 + 900*cos(x)*sin(x) *x - 32*cos(x)*sin(x) - 288*cos(x)*x 7 3 7 5 3 5 + 6*sin(x) *x - 39*sin(x) *x - 36*sin(x) *x + 318*sin(x) *x 3 3 3 3 + 72*sin(x) *x - 672*sin(x) *x - 48*sin(x)*x + 384*sin(x)*x)/2304 6 7 *k + O(k ) taylorcombine(hugo^2 + otto^2); 2 2 7 cos(x) + sin(x) + O(k ) clearrules elliptic_rules; clear sn, cn, dn, epsilon; %%% showtime; comment That's all, folks; end; 4: 4: 4: 4: 4: 4: 4: 4: 4: Time for test: 4800 ms, plus GC time: 160 ms 5: 5: Quitting Mon Jun 7 18:22:34 MET DST 1999