1comment 2 Test and demonstration file for the Taylor expansion package, 3 by Rainer M. Schoepf. Works with version 2.2a (01-Apr-2000); 4 5%%% showtime; 6 7on errcont; % disable interruption on errors 8 9comment Simple Taylor expansion; 10 11xx := taylor (e**x, x, 0, 4); 12yy := taylor (e**y, y, 0, 4); 13 14comment Basic operations, i.e. addition, subtraction, multiplication, 15 and division are possible: this is not done automatically if 16 the switch TAYLORAUTOCOMBINE is OFF. In this case it is 17 necessary to use taylorcombine; 18 19taylorcombine (xx**2); 20taylorcombine (ws - xx); 21taylorcombine (xx**3); 22 23comment The result is again a Taylor kernel; 24 25if taylorseriesp ws then write "OK"; 26 27comment It is not possible to combine Taylor kernels that were 28 expanded with respect to different variables; 29 30taylorcombine (xx**yy); 31 32comment But we can take the exponential or the logarithm 33 of a Taylor kernel; 34 35taylorcombine (e**xx); 36taylorcombine log ws; 37 38comment A more complicated example; 39 40hugo := taylor(log(1/(1-x)),x,0,5); 41taylorcombine(exp(hugo/(1+hugo))); 42 43comment We may try to expand about another point; 44 45taylor (xx, x, 1, 2); 46 47comment Arc tangent is one of the functions this package knows of; 48 49xxa := taylorcombine atan ws; 50 51comment The trigonometric functions; 52 53taylor (tan x / x, x, 0, 2); 54 55taylorcombine sin ws; 56 57taylor (cot x / x, x, 0, 4); 58 59comment The poles of these functions are correctly handled; 60 61taylor(tan x,x,pi/2,0); 62 63taylor(tan x,x,pi/2,3); 64 65comment Expansion with respect to more than one kernel is possible; 66 67xy := taylor (e**(x+y), x, 0, 2, y, 0, 2); 68 69taylorcombine (ws**2); 70 71comment We take the inverse and convert back to REDUCE's standard 72 representation; 73 74taylorcombine (1/ws); 75taylortostandard ws; 76 77comment Some examples of Taylor kernel divsion; 78 79xx1 := taylor (sin (x), x, 0, 4); 80taylorcombine (xx/xx1); 81taylorcombine (1/xx1); 82 83tt1 := taylor (exp (x), x, 0, 3); 84tt2 := taylor (sin (x), x, 0, 3); 85tt3 := taylor (1 + tt2, x, 0, 3); 86taylorcombine(tt1/tt2); 87taylorcombine(tt1/tt3); 88taylorcombine(tt2/tt1); 89taylorcombine(tt3/tt1); 90 91 92comment Here's what I call homogeneous expansion; 93 94xx := taylor (e**(x*y), {x,y}, 0, 2); 95xx1 := taylor (sin (x+y), {x,y}, 0, 2); 96xx2 := taylor (cos (x+y), {x,y}, 0, 2); 97temp := taylorcombine (xx/xx2); 98taylorcombine (ws*xx2); 99 100comment The following shows a principal difficulty: 101 since xx1 is symmetric in x and y but has no constant term 102 it is impossible to compute 1/xx1; 103 104taylorcombine (1/xx1); 105 106comment Substitution in Taylor expressions is possible; 107 108sub (x=z, xy); 109 110comment Expression dependency in substitution is detected; 111 112sub (x=y, xy); 113 114comment It is possible to replace a Taylor variable by a constant; 115 116sub (x=4, xy); 117 118sub (x=4, xx1); 119 120sub (y=0, ws); 121 122comment This package has three switches: 123 TAYLORKEEPORIGINAL, TAYLORAUTOEXPAND, and TAYLORAUTOCOMBINE; 124 125on taylorkeeporiginal; 126 127temp := taylor (e**(x+y), x, 0, 5); 128 129taylorcombine (log (temp)); 130 131taylororiginal ws; 132 133taylorcombine (temp * e**x); 134 135on taylorautoexpand; 136 137taylorcombine ws; 138 139taylororiginal ws; 140 141taylorcombine (xx1 / x); 142 143on taylorautocombine; 144 145xx / xx2; 146 147ws * xx2; 148 149comment Another example that shows truncation if Taylor kernels 150 of different expansion order are combined; 151 152comment First we increase the number of terms to be printed; 153 154taylorprintterms := all; 155 156p := taylor (x**2 + 2, x, 0, 10); 157p - x**2; 158p - taylor (x**2, x, 0, 5); 159taylor (p - x**2, x, 0, 6); 160off taylorautocombine; 161taylorcombine(p-x**2); 162taylorcombine(p - taylor(x**2,x,0,5)); 163 164comment Switch back to finite number of terms; 165 166taylorprintterms := 6; 167 168comment Some more examples; 169 170taylor(1/(1+y^4+x^2*y^2+x^4),{x,y},0,6); 171 172taylor ((1 + x)**n, x, 0, 3); 173 174taylor (e**(-a*t) * (1 + sin(t)), t, 0, 4); 175 176operator f; 177 178taylor (1 + f(t), t, 0, 3); 179 180taylor(f(sqrt(x^2+y^2)),x,x0,4,y,y0,4); 181 182clear f; 183 184taylor (sqrt(1 + a*x + sin(x)), x, 0, 3); 185 186taylorcombine (ws**2); 187 188taylor (sqrt(1 + x), x, 0, 5); 189 190taylor ((cos(x) - sec(x))^3, x, 0, 5); 191 192taylor ((cos(x) - sec(x))^-3, x, 0, 5); 193 194taylor (sqrt(1 - k^2*sin(x)^2), x, 0, 6); 195 196taylor (sin(x + y), x, 0, 3, y, 0, 3); 197 198taylor (e^x - 1 - x,x,0,6); 199 200taylorcombine sqrt ws; 201 202taylor(sin(x)/x,x,1,2); 203 204taylor((sqrt(4+h)-2)/h,h,0,5); 205 206taylor((sqrt(x)-2)/(4-x),x,4,2); 207 208taylor((sqrt(y+4)-2)/(-y),y,0,2); 209 210taylor(x*tanh(x)/(sqrt(1-x^2)-1),x,0,3); 211 212taylor((e^(5*x)-2*x)^(1/x),x,0,2); 213 214taylor(sin x/cos x,x,pi/2,3); 215 216taylor(log x*sin(x^2)/(x*sinh x),x,0,2); 217 218taylor(1/x-1/sin x,x,0,2); 219 220taylor(tan x/log cos x,x,pi/2,2); 221 222taylor(log(x^2/(x^2-a)),x,0,3); 223 224 225comment Three more complicated examples contributed by Stan Kameny; 226 227zz2 := (z*(z-2*pi*i)*(z-pi*i/2)^2)/(sinh z-i); 228dz2 := df(zz2,z); 229z0 := pi*i/2; 230taylor(dz2,z,z0,6); 231 232zz3:=(z*(z-2*pi)*(z-pi/2)^2)/(sin z-1); 233dz3 := df(zz3,z); 234z1 := pi/2; 235taylor(dz3,z,z1,6); 236 237taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,6); 238 239taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,0); 240 241 242comment If the expansion point is not constant, it has to be taken 243 care of in differentation, as the following examples show; 244 245taylor(sin(x+a),x,a,8); 246df(ws,a); 247taylor(cos(x+a),x,a,7); 248 249 250comment A problem are non-analytical terms: rational powers and 251 logarithmic terms can be handled, but other types of essential 252 singularities cannot; 253 254taylor(sqrt(x),x,0,2); 255 256taylor(asinh(1/x),x,0,5); 257 258taylor(e**(1/x),x,0,2); 259 260comment Another example for non-integer powers; 261 262sub (y = sqrt (x), yy); 263 264comment Expansion about infinity is possible in principle...; 265 266taylor (e**(1/x), x, infinity, 5); 267xi := taylor (sin (1/x), x, infinity, 5); 268 269y1 := taylor(x/(x-1), x, infinity, 3); 270z := df(y1, x); 271 272comment ...but far from being perfect; 273 274taylor (1 / sin (x), x, infinity, 5); 275 276clear z; 277 278comment You may access the expansion with the PART operator; 279 280part(yy,0); 281part(yy,1); 282part(yy,4); 283part(yy,6); 284 285 286comment The template of a Taylor kernel can be extracted; 287 288taylortemplate yy; 289 290taylortemplate xxa; 291 292taylortemplate xi; 293 294taylortemplate xy; 295 296taylortemplate xx1; 297 298comment Here is a slightly less trivial example; 299 300exp := (sin (x) * sin (y) / (x * y))**2; 301 302taylor (exp, x, 0, 1, y, 0, 1); 303taylor (exp, x, 0, 2, y, 0, 2); 304 305tt := taylor (exp, {x,y}, 0, 2); 306 307comment An example that uses factorization; 308 309on factor; 310 311ff := y**5 - 1; 312 313zz := sub (y = taylor(e**x, x, 0, 3), ff); 314 315on exp; 316 317zz; 318 319comment A simple example of Taylor kernel differentiation; 320 321hugo := taylor(e^x,x,0,5); 322df(hugo^2,x); 323 324comment The following shows the (limited) capabilities to integrate 325 Taylor kernels. Only simple cases are supported, otherwise 326 a warning is printed and the Taylor kernels are converted to 327 standard representation; 328 329zz := taylor (sin x, x, 0, 5); 330ww := taylor (cos y, y, 0, 5); 331 332int (zz, x); 333int (ww, x); 334int (zz + ww, x); 335 336 337comment And here we present Taylor series reversion. 338 We start with the example given by Knuth for the algorithm; 339 340taylor (t - t**2, t, 0, 5); 341 342taylorrevert (ws, t, x); 343 344tan!-series := taylor (tan x, x, 0, 5); 345taylorrevert (tan!-series, x, y); 346atan!-series:=taylor (atan y, y, 0, 5); 347 348tmp := taylor (e**x, x, 0, 5); 349 350taylorrevert (tmp, x, y); 351 352taylor (log y, y, 1, 5); 353 354 355comment The following example calculates the perturbation expansion 356 of the root x = 20 of the following polynomial in terms of 357 EPS, in ROUNDED mode; 358 359poly := for r := 1 : 20 product (x - r); 360 361on rounded; 362 363tpoly := taylor (poly, x, 20, 4); 364 365taylorrevert (tpoly, x, eps); 366 367comment Some more examples using rounded mode; 368 369precision 13; 370 371taylor(sin x/x,x,0,4); 372 373taylor(sin x,x,pi/2,4); 374 375taylor(tan x,x,pi/2,4); 376 377off rounded; 378 379 380comment An example that involves computing limits of type 0/0 if 381 expansion is done via differentiation; 382 383 384taylor(sqrt((e^x - 1)/x),x,0,15); 385 386comment An example that involves intermediate non-analytical terms 387 which cancel entirely; 388 389taylor(x^(5/2)/(log(x+1)*tan(x^(3/2))),x,0,5); 390 391comment Other examples involving non-analytical terms; 392 393taylor(log(e^x-1),x,0,5); 394 395taylor(e^(1/x)*(e^x-1),x,0,5); 396 397taylor(log(x)*x^10,x,0,5); 398 399taylor(log(x)*x^10,x,0,11); 400 401taylor(log(x-a)/((a-b)*(a-c)) + log(2(x-b))/((b-c)*(b-a)) 402 + log(x-c)/((c-a)*(c-b)),x,infinity,2); 403 404ss := (sqrt(x^(2/5) +1) - x^(1/3)-1)/x^(1/3); 405taylor(exp ss,x,0,2); 406 407taylor(exp sub(x=x^15,ss),x,0,2); 408 409taylor(dilog(x),x,0,4); 410 411taylor(ei(x),x,0,4); 412 413comment In the following we demonstrate the possibiblity to compute the 414 expansion of a function which is given by a simple first order 415 differential equation: the function myexp(x) is exp(-x^2); 416 417operator myexp,myerf; 418let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 419 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 420taylor(myexp(x),x,0,5); 421taylor(myerf(x),x,0,5); 422clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 423 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 424clear myexp,erf; 425 426%%% showtime; 427 428comment There are two special operators, implicit_taylor and 429 inverse_taylor, to compute the Taylor expansion of implicit 430 or inverse functions; 431 432 433implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5); 434 435implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20); 436 437implicit_taylor(x+y^3-y,x,y,0,0,8); 438 439implicit_taylor(x+y^3-y,x,y,0,1,5); 440 441implicit_taylor(x+y^3-y,x,y,0,-1,5); 442 443implicit_taylor(y*e^y-x,x,y,0,0,5); 444 445comment This is the function exp(-1/x^2), which has an essential 446 singularity at the point 0; 447 448implicit_taylor(x^2*log y+1,x,y,0,0,3); 449 450inverse_taylor(exp(x)-1,x,y,0,8); 451 452inverse_taylor(exp(x),x,y,0,5); 453 454inverse_taylor(sqrt(x),x,y,0,5); 455 456inverse_taylor(log(1+x),x,y,0,5); 457 458inverse_taylor((e^x-e^(-x))/2,x,y,0,5); 459 460comment In the next two cases the inverse functions have a branch 461 point, therefore the computation fails; 462 463inverse_taylor((e^x+e^(-x))/2,x,y,0,5); 464 465inverse_taylor(exp(x^2-1),x,y,0,5); 466 467inverse_taylor(exp(sqrt(x))-1,x,y,0,5); 468 469inverse_taylor(x*exp(x),x,y,0,5); 470 471 472%%% showtime; 473 474 475comment An application is the problem posed by Prof. Stanley: 476 we prove that the finite difference expression below 477 corresponds to the given derivative expression; 478 479operator diff,a,f,gg; % We use gg to avoid conflict with high energy 480 % physics operator. 481 482let diff(~f,~arg) => df(f,arg); 483 484derivative_expression := 485diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) + 486diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ; 487 488finite_difference_expression := 489+a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 490+a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 491+a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 492+a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 493-f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 494-f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 495-f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 496-a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 497-gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 498-gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 499-gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 500-a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 501+f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 502+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 503+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 504+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 505-gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 506+gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2) 507-gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 508+gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 509-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 510-a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 511-a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 512+gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 513+a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 514+a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 515-gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 516+gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 517-a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 518-a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 519+gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 520+a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 521+f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 522-f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2) 523+f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 524-f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 525+a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 526+a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 527+a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 528+a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 529-f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 530-f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 531-f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 532-a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 533-gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 534-gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 535-gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 536-a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 537+f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 538+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 539+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 540+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 541-gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 542+gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2) 543-gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 544+gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 545-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 546-a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 547-a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 548+gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 549+a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 550+a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 551-gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 552+gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 553-a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 554-a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 555+gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 556+a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 557+f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 558-f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2) 559+f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 560-f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 561+f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 562+f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 563+f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 564+a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 565-f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 566-f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 567-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 568-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 569-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 570-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 571+f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 572+f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 573-f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2) 574+a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 575+a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 576+a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 577+a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 578-f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 579-f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 580-f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 581-a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 582-gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 583-gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 584-gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 585-a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 586+f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 587+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 588+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 589+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 590-gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 591+gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2) 592-gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 593+gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 594-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 595-a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 596-a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 597+gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 598+a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 599+a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 600-gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 601+gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 602-a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 603-a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 604+gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 605+a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 606+f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 607-f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2) 608+f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 609-f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 610+a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 611+a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 612+a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 613+a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 614-f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 615-f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 616-f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 617-a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 618-gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 619-gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 620-gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 621-a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 622+f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 623+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 624+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 625+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 626-gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 627+gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2) 628-gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 629+gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 630-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 631-a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 632-a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 633+gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 634+a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 635+a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 636-gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 637+gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 638-a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 639-a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 640+gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 641+a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 642+f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 643-f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2) 644+f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 645-f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 646+f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 647+f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 648+f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 649+a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 650-f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 651-f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 652-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 653-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 654-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 655-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 656+f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 657+f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 658-f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2) 659+f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2) 660+a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2) 661-f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2) 662+f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2) 663+a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2) 664-f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2) 665-a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$ 666 667comment We define abbreviations for the partial derivatives; 668 669operator ax,ay,fx,fy,gx,gy; 670operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 671operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 672operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy, 673 gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 674operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 675operator gxxxxyy,gxxxyyy,gxxyyyy; 676 677operator_diff_rules := { 678 df(a(~x,~y),~x) => ax(x,y), 679 df(a(~x,~y),~y) => ay(x,y), 680 df(f(~x,~y),~x) => fx(x,y), 681 df(f(~x,~y),~y) => fy(x,y), 682 df(gg(~x,~y),~x) => gx(x,y), 683 df(gg(~x,~y),~y) => gy(x,y), 684 685 df(ax(~x,~y),~x) => axx(x,y), 686 df(ax(~x,~y),~y) => axy(x,y), 687 df(ay(~x,~y),~x) => axy(x,y), 688 df(ay(~x,~y),~y) => ayy(x,y), 689 df(fx(~x,~y),~x) => fxx(x,y), 690 df(fx(~x,~y),~y) => fxy(x,y), 691 df(fy(~x,~y),~x) => fxy(x,y), 692 df(fy(~x,~y),~y) => fyy(x,y), 693 df(gx(~x,~y),~x) => gxx(x,y), 694 df(gx(~x,~y),~y) => gxy(x,y), 695 df(gy(~x,~y),~x) => gxy(x,y), 696 df(gy(~x,~y),~y) => gyy(x,y), 697 698 df(axx(~x,~y),~x) => axxx(x,y), 699 df(axy(~x,~y),~x) => axxy(x,y), 700 df(ayy(~x,~y),~x) => axyy(x,y), 701 df(ayy(~x,~y),~y) => ayyy(x,y), 702 df(fxx(~x,~y),~x) => fxxx(x,y), 703 df(fxy(~x,~y),~x) => fxxy(x,y), 704 df(fxy(~x,~y),~y) => fxyy(x,y), 705 df(fyy(~x,~y),~x) => fxyy(x,y), 706 df(fyy(~x,~y),~y) => fyyy(x,y), 707 df(gxx(~x,~y),~x) => gxxx(x,y), 708 df(gxx(~x,~y),~y) => gxxy(x,y), 709 df(gxy(~x,~y),~x) => gxxy(x,y), 710 df(gxy(~x,~y),~y) => gxyy(x,y), 711 df(gyy(~x,~y),~x) => gxyy(x,y), 712 df(gyy(~x,~y),~y) => gyyy(x,y), 713 714 df(axyy(~x,~y),~x) => axxyy(x,y), 715 df(axxy(~x,~y),~x) => axxxy(x,y), 716 df(ayyy(~x,~y),~x) => axyyy(x,y), 717 df(fxxy(~x,~y),~x) => fxxxy(x,y), 718 df(fxyy(~x,~y),~x) => fxxyy(x,y), 719 df(fyyy(~x,~y),~x) => fxyyy(x,y), 720 df(gxxx(~x,~y),~x) => gxxxx(x,y), 721 df(gxxy(~x,~y),~x) => gxxxy(x,y), 722 df(gxyy(~x,~y),~x) => gxxyy(x,y), 723 df(gyyy(~x,~y),~x) => gxyyy(x,y), 724 df(gyyy(~x,~y),~y) => gyyyy(x,y), 725 726 df(axxyy(~x,~y),~x) => axxxyy(x,y), 727 df(axyyy(~x,~y),~x) => axxyyy(x,y), 728 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 729 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 730 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 731 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 732 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 733 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 734 735 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 736 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 737 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y) 738}; 739 740let operator_diff_rules; 741 742texp := taylor (finite_difference_expression, dx, 0, 1, dy, 0, 1); 743 744comment You may also try to expand further but this needs a lot 745 of CPU time. Therefore the following line is commented out; 746 747%texp := taylor (finite_difference_expression, dx, 0, 2, dy, 0, 2); 748 749factor dx,dy; 750 751result := taylortostandard texp; 752 753derivative_expression - result; 754 755 756clear diff(~f,~arg); 757clearrules operator_diff_rules; 758clear diff,a,f,gg; 759clear ax,ay,fx,fy,gx,gy; 760clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 761clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 762clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 763clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 764clear gxxxxyy,gxxxyyy,gxxyyyy; 765 766taylorprintterms := 5; 767 768off taylorautoexpand,taylorkeeporiginal; 769 770%%% showtime; 771 772comment An example provided by Alan Barnes: elliptic functions; 773 774% Jacobi's elliptic functions 775% sn(x,k), cn(x,k), dn(x,k). 776% The modulus and complementary modulus are denoted by K and K!' 777% usually written mathematically as k and k' respectively 778% 779% epsilon(x,k) is the incomplete elliptic integral of the second kind 780% usually written mathematically as E(x,k) 781% 782% KK(k) is the complete elliptic integral of the first kind 783% usually written mathematically as K(k) 784% K(k) = arcsn(1,k) 785% KK!'(k) is the complementary complete integral 786% usually written mathematically as K'(k) 787% NB. K'(k) = K(k') 788% EE(k) is the complete elliptic integral of the second kind 789% usually written mathematically as E(k) 790% EE!'(k) is the complementary complete integral 791% usually written mathematically as E'(k) 792% NB. E'(k) = E(k') 793 794operator sn, cn, dn, epsilon; 795 796elliptic_rules := { 797% Differentiation rules for basic functions 798 df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 799 df(cn(~x,~k),~x) => -sn(x,k)*dn(x,k), 800 df(dn(~x,~k),~x) => -k^2*sn(x,k)*cn(x,k), 801 df(epsilon(~x,~k),~x)=> dn(x,k)^2, 802 803% k-derivatives 804% DF Lawden Elliptic Functions & Applications Springer (1989) 805 df(sn(~x,~k),~k) => (k*sn(x,k)*cn(x,k)^2 806 -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2) 807 + x*cn(x,k)*dn(x,k)/k, 808 df(cn(~x,~k),~k) => (-k*sn(x,k)^2*cn(x,k) 809 +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2) 810 - x*sn(x,k)*dn(x,k)/k, 811 df(dn(~x,~k),~k) => k*(-sn(x,k)^2*dn(x,k) 812 +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2) 813 - k*x*sn(x,k)*cn(x,k), 814 df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k) 815 -epsilon(x,k)*cn(x,k)^2)/(1-k^2) 816 -k*x*sn(x,k)^2, 817 818% parity properties 819 sn(-~x,~k) => -sn(x,k), 820 cn(-~x,~k) => cn(x,k), 821 dn(-~x,~k) => dn(x,k), 822 epsilon(-~x,~k) => -epsilon(x,k), 823 sn(~x,-~k) => sn(x,k), 824 cn(~x,-~k) => cn(x,k), 825 dn(~x,-~k) => dn(x,k), 826 epsilon(~x,-~k) => epsilon(x,k), 827 828 829% behaviour at zero 830 sn(0,~k) => 0, 831 cn(0,~k) => 1, 832 dn(0,~k) => 1, 833 epsilon(0,~k) => 0, 834 835 836% degenerate cases of modulus 837 sn(~x,0) => sin(x), 838 cn(~x,0) => cos(x), 839 dn(~x,0) => 1, 840 epsilon(~x,0) => x, 841 842 sn(~x,1) => tanh(x), 843 cn(~x,1) => 1/cosh(x), 844 dn(~x,1) => 1/cosh(x), 845 epsilon(~x,1) => tanh(x) 846}; 847 848let elliptic_rules; 849 850hugo := taylor(sn(x,k),k,0,6); 851otto := taylor(cn(x,k),k,0,6); 852taylorcombine(hugo^2 + otto^2); 853 854clearrules elliptic_rules; 855 856clear sn, cn, dn, epsilon; 857 858%%% showtime; 859 860comment That's all, folks; 861 862end; 863