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(dilog(x),x,1,4); 412 413taylor(dilog(x),x,-1,4); 414 415taylor(Ei(x),x,0,4); 416 417taylor(Ei(x),x,1,4); 418 419taylor(Ei(x),x,-1,4); 420 421COMMENT In the following we demonstrate the possibiblity to compute the 422 expansion of a function which is given by a simple first order 423 differential equation: the function myexp(x) is exp(-x^2); 424 425operator myexp,myerf; 426let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 427 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 428taylor(myexp(x),x,0,5); 429taylor(myerf(x),x,0,5); 430clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 431 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 432clear myexp,myerf; 433 434%%% showtime; 435 436COMMENT There are two special operators, implicit_taylor and 437 inverse_taylor, to compute the Taylor expansion of implicit 438 or inverse functions; 439 440 441implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5); 442 443implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20); 444 445implicit_taylor(x+y^3-y,x,y,0,0,8); 446 447implicit_taylor(x+y^3-y,x,y,0,1,5); 448 449implicit_taylor(x+y^3-y,x,y,0,-1,5); 450 451implicit_taylor(y*e^y-x,x,y,0,0,5); 452 453COMMENT This is the function exp(-1/x^2), which has an essential 454 singularity at the point 0; 455 456implicit_taylor(x^2*log y+1,x,y,0,0,3); 457 458inverse_taylor(exp(x)-1,x,y,0,8); 459 460inverse_taylor(exp(x),x,y,0,5); 461 462inverse_taylor(sqrt(x),x,y,0,5); 463 464inverse_taylor(log(1+x),x,y,0,5); 465 466inverse_taylor((e^x-e^(-x))/2,x,y,0,5); 467 468COMMENT In the next two cases the inverse functions have a branch 469 point, therefore the computation fails; 470 471inverse_taylor((e^x+e^(-x))/2,x,y,0,5); 472 473inverse_taylor(exp(x^2-1),x,y,0,5); 474 475inverse_taylor(exp(sqrt(x))-1,x,y,0,5); 476 477inverse_taylor(x*exp(x),x,y,0,5); 478 479 480%%% showtime; 481 482 483COMMENT An application is the problem posed by Prof. Stanley: 484 we prove that the finite difference expression below 485 corresponds to the given derivative expression; 486 487operator diff,a,f,gg; % We use gg to avoid conflict with high energy 488 % physics operator. 489 490let diff(~f,~arg) => df(f,arg); 491 492derivative_expression := 493diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) + 494diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ; 495 496finite_difference_expression := 497+a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 498+a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 499+a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 500+a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 501-f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 502-f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 503-f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 504-a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 505-gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 506-gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 507-gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 508-a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 509+f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 510+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 511+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 512+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 513-gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 514+gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2) 515-gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 516+gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 517-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 518-a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 519-a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 520+gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 521+a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 522+a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 523-gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 524+gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 525-a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 526-a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 527+gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 528+a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 529+f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 530-f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2) 531+f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 532-f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 533+a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 534+a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 535+a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 536+a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 537-f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 538-f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 539-f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 540-a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 541-gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 542-gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 543-gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 544-a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 545+f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 546+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 547+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 548+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 549-gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 550+gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2) 551-gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 552+gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 553-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 554-a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 555-a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 556+gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 557+a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 558+a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 559-gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 560+gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 561-a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 562-a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 563+gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 564+a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 565+f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 566-f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2) 567+f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 568-f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 569+f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 570+f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 571+f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 572+a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 573-f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 574-f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 575-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 576-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 577-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 578-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 579+f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 580+f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 581-f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2) 582+a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 583+a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 584+a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 585+a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 586-f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 587-f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 588-f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 589-a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 590-gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 591-gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 592-gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 593-a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 594+f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 595+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 596+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 597+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 598-gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 599+gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2) 600-gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 601+gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 602-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 603-a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 604-a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 605+gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 606+a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 607+a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 608-gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 609+gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 610-a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 611-a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 612+gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 613+a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 614+f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 615-f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2) 616+f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 617-f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 618+a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 619+a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 620+a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 621+a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 622-f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 623-f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 624-f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 625-a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 626-gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 627-gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 628-gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 629-a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 630+f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 631+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 632+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 633+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 634-gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 635+gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2) 636-gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 637+gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 638-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 639-a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 640-a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 641+gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 642+a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 643+a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 644-gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 645+gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 646-a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 647-a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 648+gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 649+a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 650+f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 651-f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2) 652+f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 653-f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 654+f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 655+f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 656+f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 657+a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 658-f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 659-f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 660-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 661-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 662-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 663-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 664+f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 665+f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 666-f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2) 667+f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2) 668+a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2) 669-f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2) 670+f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2) 671+a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2) 672-f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2) 673-a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$ 674 675COMMENT We define abbreviations for the partial derivatives; 676 677operator ax,ay,fx,fy,gx,gy; 678operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 679operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 680operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy, 681 gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 682operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 683operator gxxxxyy,gxxxyyy,gxxyyyy; 684 685operator_diff_rules := { 686 df(a(~x,~y),~x) => ax(x,y), 687 df(a(~x,~y),~y) => ay(x,y), 688 df(f(~x,~y),~x) => fx(x,y), 689 df(f(~x,~y),~y) => fy(x,y), 690 df(gg(~x,~y),~x) => gx(x,y), 691 df(gg(~x,~y),~y) => gy(x,y), 692 693 df(ax(~x,~y),~x) => axx(x,y), 694 df(ax(~x,~y),~y) => axy(x,y), 695 df(ay(~x,~y),~x) => axy(x,y), 696 df(ay(~x,~y),~y) => ayy(x,y), 697 df(fx(~x,~y),~x) => fxx(x,y), 698 df(fx(~x,~y),~y) => fxy(x,y), 699 df(fy(~x,~y),~x) => fxy(x,y), 700 df(fy(~x,~y),~y) => fyy(x,y), 701 df(gx(~x,~y),~x) => gxx(x,y), 702 df(gx(~x,~y),~y) => gxy(x,y), 703 df(gy(~x,~y),~x) => gxy(x,y), 704 df(gy(~x,~y),~y) => gyy(x,y), 705 706 df(axx(~x,~y),~x) => axxx(x,y), 707 df(axy(~x,~y),~x) => axxy(x,y), 708 df(ayy(~x,~y),~x) => axyy(x,y), 709 df(ayy(~x,~y),~y) => ayyy(x,y), 710 df(fxx(~x,~y),~x) => fxxx(x,y), 711 df(fxy(~x,~y),~x) => fxxy(x,y), 712 df(fxy(~x,~y),~y) => fxyy(x,y), 713 df(fyy(~x,~y),~x) => fxyy(x,y), 714 df(fyy(~x,~y),~y) => fyyy(x,y), 715 df(gxx(~x,~y),~x) => gxxx(x,y), 716 df(gxx(~x,~y),~y) => gxxy(x,y), 717 df(gxy(~x,~y),~x) => gxxy(x,y), 718 df(gxy(~x,~y),~y) => gxyy(x,y), 719 df(gyy(~x,~y),~x) => gxyy(x,y), 720 df(gyy(~x,~y),~y) => gyyy(x,y), 721 722 df(axyy(~x,~y),~x) => axxyy(x,y), 723 df(axxy(~x,~y),~x) => axxxy(x,y), 724 df(ayyy(~x,~y),~x) => axyyy(x,y), 725 df(fxxy(~x,~y),~x) => fxxxy(x,y), 726 df(fxyy(~x,~y),~x) => fxxyy(x,y), 727 df(fyyy(~x,~y),~x) => fxyyy(x,y), 728 df(gxxx(~x,~y),~x) => gxxxx(x,y), 729 df(gxxy(~x,~y),~x) => gxxxy(x,y), 730 df(gxyy(~x,~y),~x) => gxxyy(x,y), 731 df(gyyy(~x,~y),~x) => gxyyy(x,y), 732 df(gyyy(~x,~y),~y) => gyyyy(x,y), 733 734 df(axxyy(~x,~y),~x) => axxxyy(x,y), 735 df(axyyy(~x,~y),~x) => axxyyy(x,y), 736 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 737 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 738 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 739 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 740 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 741 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 742 743 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 744 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 745 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y) 746}; 747 748let operator_diff_rules; 749 750texp := taylor (finite_difference_expression, dx, 0, 1, dy, 0, 1); 751 752COMMENT You may also try to expand further but this needs a lot 753 of CPU time. Therefore the following line is commented out; 754 755%texp := taylor (finite_difference_expression, dx, 0, 2, dy, 0, 2); 756 757factor dx,dy; 758 759result := taylortostandard texp; 760 761derivative_expression - result; 762 763 764clear diff(~f,~arg); 765clearrules operator_diff_rules; 766clear diff,a,f,gg; 767clear ax,ay,fx,fy,gx,gy; 768clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 769clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 770clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 771clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 772clear gxxxxyy,gxxxyyy,gxxyyyy; 773 774taylorprintterms := 5; 775 776off taylorautoexpand,taylorkeeporiginal; 777 778%%% showtime; 779 780COMMENT An example provided by Alan Barnes: elliptic functions; 781 782% Jacobi's elliptic functions 783% sn(x,k), cn(x,k), dn(x,k). 784% The modulus and complementary modulus are denoted by K and K!' 785% usually written mathematically as k and k' respectively 786% 787% epsilon(x,k) is the incomplete elliptic integral of the second kind 788% usually written mathematically as E(x,k) 789% 790% KK(k) is the complete elliptic integral of the first kind 791% usually written mathematically as K(k) 792% K(k) = arcsn(1,k) 793% KK!'(k) is the complementary complete integral 794% usually written mathematically as K'(k) 795% NB. K'(k) = K(k') 796% EE(k) is the complete elliptic integral of the second kind 797% usually written mathematically as E(k) 798% EE!'(k) is the complementary complete integral 799% usually written mathematically as E'(k) 800% NB. E'(k) = E(k') 801 802operator sn, cn, dn, epsilon; 803 804elliptic_rules := { 805% Differentiation rules for basic functions 806 df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 807 df(cn(~x,~k),~x) => -sn(x,k)*dn(x,k), 808 df(dn(~x,~k),~x) => -k^2*sn(x,k)*cn(x,k), 809 df(epsilon(~x,~k),~x)=> dn(x,k)^2, 810 811% k-derivatives 812% DF Lawden Elliptic Functions & Applications Springer (1989) 813 df(sn(~x,~k),~k) => (k*sn(x,k)*cn(x,k)^2 814 -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2) 815 + x*cn(x,k)*dn(x,k)/k, 816 df(cn(~x,~k),~k) => (-k*sn(x,k)^2*cn(x,k) 817 +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2) 818 - x*sn(x,k)*dn(x,k)/k, 819 df(dn(~x,~k),~k) => k*(-sn(x,k)^2*dn(x,k) 820 +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2) 821 - k*x*sn(x,k)*cn(x,k), 822 df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k) 823 -epsilon(x,k)*cn(x,k)^2)/(1-k^2) 824 -k*x*sn(x,k)^2, 825 826% parity properties 827 sn(-~x,~k) => -sn(x,k), 828 cn(-~x,~k) => cn(x,k), 829 dn(-~x,~k) => dn(x,k), 830 epsilon(-~x,~k) => -epsilon(x,k), 831 sn(~x,-~k) => sn(x,k), 832 cn(~x,-~k) => cn(x,k), 833 dn(~x,-~k) => dn(x,k), 834 epsilon(~x,-~k) => epsilon(x,k), 835 836 837% behaviour at zero 838 sn(0,~k) => 0, 839 cn(0,~k) => 1, 840 dn(0,~k) => 1, 841 epsilon(0,~k) => 0, 842 843 844% degenerate cases of modulus 845 sn(~x,0) => sin(x), 846 cn(~x,0) => cos(x), 847 dn(~x,0) => 1, 848 epsilon(~x,0) => x, 849 850 sn(~x,1) => tanh(x), 851 cn(~x,1) => 1/cosh(x), 852 dn(~x,1) => 1/cosh(x), 853 epsilon(~x,1) => tanh(x) 854}; 855 856let elliptic_rules; 857 858hugo := taylor(sn(x,k),k,0,6); 859otto := taylor(cn(x,k),k,0,6); 860taylorcombine(hugo^2 + otto^2); 861 862clearrules elliptic_rules; 863 864clear sn, cn, dn, epsilon; 865 866COMMENT Examples of gamma function and its derivatives; 867 868taylor(gamma(x),x,1,3); 869 870taylor(gamma(x),x,1/2,3); 871 872taylor(gamma(x),x,a,3); 873 874taylor(gamma(x),x,-1,3); 875 876taylor(gamma(1+x),x,0, 6); 877 878COMMENT Test printing for negative expansion point; 879 880taylor(sin x,x, -1, 6); 881 882taylor(sin x,x, -1/2, 6); 883 884%%% showtime; 885 886COMMENT That's all, folks; 887 888end; 889