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 6%%% showtime; 7 8on errcont; 9 10 % disable interruption on errors 11 12COMMENT Simple Taylor expansion; 13 14 15xx := taylor (e**x, x, 0, 4); 16 17 18 1 2 1 3 1 4 5 19xx := 1 + x + ---*x + ---*x + ----*x + O(x ) 20 2 6 24 21 22yy := taylor (e**y, y, 0, 4); 23 24 25 1 2 1 3 1 4 5 26yy := 1 + y + ---*y + ---*y + ----*y + O(y ) 27 2 6 24 28 29 30COMMENT Basic operations, i.e. addition, subtraction, multiplication, 31 and division are possible: this is not done automatically if 32 the switch TAYLORAUTOCOMBINE is OFF. In this case it is 33 necessary to use taylorcombine; 34 35 36taylorcombine (xx**2); 37 38 39 2 4 3 2 4 5 401 + 2*x + 2*x + ---*x + ---*x + O(x ) 41 3 3 42 43taylorcombine (ws - xx); 44 45 46 3 2 7 3 5 4 5 47x + ---*x + ---*x + ---*x + O(x ) 48 2 6 8 49 50taylorcombine (xx**3); 51 52 53 9 2 9 3 27 4 5 541 + 3*x + ---*x + ---*x + ----*x + O(x ) 55 2 2 8 56 57 58COMMENT The result is again a Taylor kernel; 59 60 61if taylorseriesp ws then write "OK"; 62 63 64OK 65 66 67COMMENT It is not possible to combine Taylor kernels that were 68 expanded with respect to different variables; 69 70 71taylorcombine (xx**yy); 72 73 74 1 2 1 3 1 4 5 75(1 + x + ---*x + ---*x + ----*x + O(x )) 76 2 6 24 77 78 1 2 1 3 1 4 5 79**(1 + y + ---*y + ---*y + ----*y + O(y )) 80 2 6 24 81 82 83COMMENT But we can take the exponential or the logarithm 84 of a Taylor kernel; 85 86 87taylorcombine (e**xx); 88 89 90 2 5*e 3 5*e 4 5 91e + e*x + e*x + -----*x + -----*x + O(x ) 92 6 8 93 94taylorcombine log ws; 95 96 97 1 2 1 3 1 4 5 981 + x + ---*x + ---*x + ----*x + O(x ) 99 2 6 24 100 101 102COMMENT A more complicated example; 103 104 105hugo := taylor(log(1/(1-x)),x,0,5); 106 107 108 1 2 1 3 1 4 1 5 6 109hugo := x + ---*x + ---*x + ---*x + ---*x + O(x ) 110 2 3 4 5 111 112taylorcombine(exp(hugo/(1+hugo))); 113 114 115 1 4 6 1161 + x + ----*x + O(x ) 117 12 118 119 120COMMENT We may try to expand about another point; 121 122 123taylor (xx, x, 1, 2); 124 125 126 65 8 5 2 3 127---- + ---*(x - 1) + ---*(x - 1) + O((x - 1) ) 128 24 3 4 129 130 131COMMENT Arc tangent is one of the functions this package knows of; 132 133 134xxa := taylorcombine atan ws; 135 136 137 65 1536 2933040 2 3 138xxa := atan(----) + ------*(x - 1) - ----------*(x - 1) + O((x - 1) ) 139 24 4801 23049601 140 141 142COMMENT The trigonometric functions; 143 144 145taylor (tan x / x, x, 0, 2); 146 147 148 1 2 3 1491 + ---*x + O(x ) 150 3 151 152 153taylorcombine sin ws; 154 155 156 cos(1) 2 3 157sin(1) + --------*x + O(x ) 158 3 159 160 161taylor (cot x / x, x, 0, 4); 162 163 164 -2 1 1 2 2 4 5 165x - --- - ----*x - -----*x + O(x ) 166 3 45 945 167 168 169COMMENT The poles of these functions are correctly handled; 170 171 172taylor(tan x,x,pi/2,0); 173 174 175 pi -1 pi 176 - (x - ----) + O(x - ----) 177 2 2 178 179 180taylor(tan x,x,pi/2,3); 181 182 183 pi -1 1 pi 1 pi 3 pi 4 184 - (x - ----) + ---*(x - ----) + ----*(x - ----) + O((x - ----) ) 185 2 3 2 45 2 2 186 187 188COMMENT Expansion with respect to more than one kernel is possible; 189 190 191xy := taylor (e**(x+y), x, 0, 2, y, 0, 2); 192 193 194 1 2 3 3 195xy := 1 + y + ---*y + x + y*x + (4 terms) + O(x ,y ) 196 2 197 198 199taylorcombine (ws**2); 200 201 202 2 3 3 2031 + 2*y + 2*y + 2*x + 4*y*x + (4 terms) + O(x ,y ) 204 205 206COMMENT We take the inverse and convert back to REDUCE's standard 207 representation; 208 209 210taylorcombine (1/ws); 211 212 213 2 3 3 2141 - 2*x + 2*x - 2*y + 4*y*x + (4 terms) + O(x ,y ) 215 216taylortostandard ws; 217 218 219 2 2 2 2 2 2 2204*x *y - 4*x *y + 2*x - 4*x*y + 4*x*y - 2*x + 2*y - 2*y + 1 221 222 223COMMENT Some examples of Taylor kernel divsion; 224 225 226xx1 := taylor (sin (x), x, 0, 4); 227 228 229 1 3 5 230xx1 := x - ---*x + O(x ) 231 6 232 233taylorcombine (xx/xx1); 234 235 236 -1 2 1 2 3 237x + 1 + ---*x + ---*x + O(x ) 238 3 3 239 240taylorcombine (1/xx1); 241 242 243 -1 1 3 244x + ---*x + O(x ) 245 6 246 247 248tt1 := taylor (exp (x), x, 0, 3); 249 250 251 1 2 1 3 4 252tt1 := 1 + x + ---*x + ---*x + O(x ) 253 2 6 254 255tt2 := taylor (sin (x), x, 0, 3); 256 257 258 1 3 4 259tt2 := x - ---*x + O(x ) 260 6 261 262tt3 := taylor (1 + tt2, x, 0, 3); 263 264 265 1 3 4 266tt3 := 1 + x - ---*x + O(x ) 267 6 268 269taylorcombine(tt1/tt2); 270 271 272 -1 2 2 273x + 1 + ---*x + O(x ) 274 3 275 276taylorcombine(tt1/tt3); 277 278 279 1 2 1 3 4 2801 + ---*x - ---*x + O(x ) 281 2 6 282 283taylorcombine(tt2/tt1); 284 285 286 2 1 3 4 287x - x + ---*x + O(x ) 288 3 289 290taylorcombine(tt3/tt1); 291 292 293 1 2 1 3 4 2941 - ---*x + ---*x + O(x ) 295 2 6 296 297 298 299COMMENT Here's what I call homogeneous expansion; 300 301 302xx := taylor (e**(x*y), {x,y}, 0, 2); 303 304 305 3 306xx := 1 + y*x + O({x,y} ) 307 308xx1 := taylor (sin (x+y), {x,y}, 0, 2); 309 310 311 3 312xx1 := y + x + O({x,y} ) 313 314xx2 := taylor (cos (x+y), {x,y}, 0, 2); 315 316 317 1 2 1 2 3 318xx2 := 1 - ---*y - y*x - ---*x + O({x,y} ) 319 2 2 320 321temp := taylorcombine (xx/xx2); 322 323 324 1 2 1 2 3 325temp := 1 + ---*y + 2*y*x + ---*x + O({x,y} ) 326 2 2 327 328taylorcombine (ws*xx2); 329 330 331 3 3321 + y*x + O({x,y} ) 333 334 335COMMENT The following shows a principal difficulty: 336 since xx1 is symmetric in x and y but has no constant term 337 it is impossible to compute 1/xx1; 338 339 340taylorcombine (1/xx1); 341 342 343***** Not a unit in argument to invtaylor 344 345 346COMMENT Substitution in Taylor expressions is possible; 347 348 349sub (x=z, xy); 350 351 352 1 2 3 3 3531 + y + ---*y + z + y*z + (4 terms) + O(z ,y ) 354 2 355 356 357COMMENT Expression dependency in substitution is detected; 358 359 360sub (x=y, xy); 361 362 363***** Invalid substitution in Taylor kernel: dependent variables y y 364 365 366COMMENT It is possible to replace a Taylor variable by a constant; 367 368 369sub (x=4, xy); 370 371 372 13 2 3 37313 + 13*y + ----*y + O(y ) 374 2 375 376 377sub (x=4, xx1); 378 379 380 3 3814 + y + O(y ) 382 383 384sub (y=0, ws); 385 386 3874 388 389 390COMMENT This package has three switches: 391 TAYLORKEEPORIGINAL, TAYLORAUTOEXPAND, and TAYLORAUTOCOMBINE; 392 393 394on taylorkeeporiginal; 395 396 397 398temp := taylor (e**(x+y), x, 0, 5); 399 400 401 y y y 402 y y e 2 e 3 e 4 6 403temp := e + e *x + ----*x + ----*x + ----*x + (1 term) + O(x ) 404 2 6 24 405 406 407taylorcombine (log (temp)); 408 409 410 6 411y + x + O(x ) 412 413 414taylororiginal ws; 415 416 417x + y 418 419 420taylorcombine (temp * e**x); 421 422 423 y y y 424 x y y e 2 e 3 e 4 6 425e *(e + e *x + ----*x + ----*x + ----*x + (1 term) + O(x )) 426 2 6 24 427 428 429on taylorautoexpand; 430 431 432 433taylorcombine ws; 434 435 436 y y 437 y y y 2 4*e 3 2*e 4 6 438e + 2*e *x + 2*e *x + ------*x + ------*x + (1 term) + O(x ) 439 3 3 440 441 442taylororiginal ws; 443 444 445 2*x + y 446e 447 448 449taylorcombine (xx1 / x); 450 451 452 -1 2 453y*x + 1 + O({x,y} ) 454 455 456on taylorautocombine; 457 458 459 460xx / xx2; 461 462 463 1 2 1 2 3 4641 + ---*y + 2*y*x + ---*x + O({x,y} ) 465 2 2 466 467 468ws * xx2; 469 470 471 3 4721 + y*x + O({x,y} ) 473 474 475COMMENT Another example that shows truncation if Taylor kernels 476 of different expansion order are combined; 477 478 479COMMENT First we increase the number of terms to be printed; 480 481 482taylorprintterms := all; 483 484 485taylorprintterms := all 486 487 488p := taylor (x**2 + 2, x, 0, 10); 489 490 491 2 11 492p := 2 + x + O(x ) 493 494p - x**2; 495 496 497 2 11 2 498(2 + x + O(x )) - x 499 500p - taylor (x**2, x, 0, 5); 501 502 503 6 5042 + O(x ) 505 506taylor (p - x**2, x, 0, 6); 507 508 509 7 5102 + O(x ) 511 512off taylorautocombine; 513 514 515taylorcombine(p-x**2); 516 517 518 11 5192 + O(x ) 520 521taylorcombine(p - taylor(x**2,x,0,5)); 522 523 524 6 5252 + O(x ) 526 527 528COMMENT Switch back to finite number of terms; 529 530 531taylorprintterms := 6; 532 533 534taylorprintterms := 6 535 536 537COMMENT Some more examples; 538 539 540taylor(1/(1+y^4+x^2*y^2+x^4),{x,y},0,6); 541 542 543 4 2 2 4 7 5441 - y - y *x - x + O({x,y} ) 545 546 547taylor ((1 + x)**n, x, 0, 3); 548 549 550 2 551 n*(n - 1) 2 n*(n - 3*n + 2) 3 4 5521 + n*x + -----------*x + ------------------*x + O(x ) 553 2 6 554 555 556taylor (e**(-a*t) * (1 + sin(t)), t, 0, 4); 557 558 559 3 2 560 a*(a - 2) 2 - a + 3*a - 1 3 5611 + ( - a + 1)*t + -----------*t + ------------------*t 562 2 6 563 564 3 2 565 a*(a - 4*a + 4) 4 5 566 + -------------------*t + O(t ) 567 24 568 569 570operator f; 571 572 573 574taylor (1 + f(t), t, 0, 3); 575 576 577 sub(t=0,df(f(t),t,2)) 2 578f(0) + 1 + sub(t=0,df(f(t),t))*t + -----------------------*t 579 2 580 581 sub(t=0,df(f(t),t,3)) 3 4 582 + -----------------------*t + O(t ) 583 6 584 585 586taylor(f(sqrt(x^2+y^2)),x,x0,4,y,y0,4); 587 588 589 2 2 2 2 590f(sqrt(x0 + y0 )) + sub(y=y0,df(f(sqrt(x0 + y )),y))*(y - y0) 591 592 2 2 593 sub(y=y0,df(f(sqrt(x0 + y )),y,2)) 2 594 + -------------------------------------*(y - y0) 595 2 596 597 2 2 598 sub(y=y0,df(f(sqrt(x0 + y )),y,3)) 3 599 + -------------------------------------*(y - y0) 600 6 601 602 2 2 603 sub(y=y0,df(f(sqrt(x0 + y )),y,4)) 4 604 + -------------------------------------*(y - y0) 605 24 606 607 2 2 608 + sub(x=x0,df(f(sqrt(x + y0 )),x))*(x - x0) + (19 terms) 609 610 5 5 611 + O((x - x0) ,(y - y0) ) 612 613 614clear f; 615 616 617 618taylor (sqrt(1 + a*x + sin(x)), x, 0, 3); 619 620 621 2 3 2 622 a + 1 - a - 2*a - 1 2 3*a + 9*a + 9*a - 1 3 4 6231 + -------*x + -----------------*x + -----------------------*x + O(x ) 624 2 8 48 625 626 627taylorcombine (ws**2); 628 629 630 1 3 4 6311 + (a + 1)*x - ---*x + O(x ) 632 6 633 634 635taylor (sqrt(1 + x), x, 0, 5); 636 637 638 1 1 2 1 3 5 4 7 5 6 6391 + ---*x - ---*x + ----*x - -----*x + -----*x + O(x ) 640 2 8 16 128 256 641 642 643taylor ((cos(x) - sec(x))^3, x, 0, 5); 644 645 646 6 6470 + O(x ) 648 649 650taylor ((cos(x) - sec(x))^-3, x, 0, 5); 651 652 653 -6 1 -4 11 -2 347 6767 2 15377 4 6 654 - x + ---*x + -----*x - ------- - --------*x - ---------*x + O(x ) 655 2 120 15120 604800 7983360 656 657 658taylor (sqrt(1 - k^2*sin(x)^2), x, 0, 6); 659 660 661 2 2 2 2 4 2 662 k 2 k *( - 3*k + 4) 4 k *( - 45*k + 60*k - 16) 6 7 6631 - ----*x + ------------------*x + ----------------------------*x + O(x ) 664 2 24 720 665 666 667taylor (sin(x + y), x, 0, 3, y, 0, 3); 668 669 670 1 3 1 2 1 2 1 2 3 4 4 671x - ---*x + y - ---*y*x - ---*y *x + ----*y *x + (2 terms) + O(x ,y ) 672 6 2 2 12 673 674 675taylor (e^x - 1 - x,x,0,6); 676 677 678 1 2 1 3 1 4 1 5 1 6 7 679---*x + ---*x + ----*x + -----*x + -----*x + O(x ) 680 2 6 24 120 720 681 682 683taylorcombine sqrt ws; 684 685 686 1 1 2 1 3 1 4 687---------*x + -----------*x + ------------*x + -------------*x 688 sqrt(2) 6*sqrt(2) 36*sqrt(2) 270*sqrt(2) 689 690 1 5 6 691 + --------------*x + O(x ) 692 2592*sqrt(2) 693 694 695taylor(sin(x)/x,x,1,2); 696 697 698 - 2*cos(1) + sin(1) 2 699sin(1) + (cos(1) - sin(1))*(x - 1) + ----------------------*(x - 1) 700 2 701 702 3 703 + O((x - 1) ) 704 705 706taylor((sqrt(4+h)-2)/h,h,0,5); 707 708 709 1 1 1 2 5 3 7 4 21 5 6 710--- - ----*h + -----*h - -------*h + --------*h - ---------*h + O(h ) 711 4 64 512 16384 131072 2097152 712 713 714taylor((sqrt(x)-2)/(4-x),x,4,2); 715 716 717 1 1 1 2 3 718 - --- + ----*(x - 4) - -----*(x - 4) + O((x - 4) ) 719 4 64 512 720 721 722taylor((sqrt(y+4)-2)/(-y),y,0,2); 723 724 725 1 1 1 2 3 726 - --- + ----*y - -----*y + O(y ) 727 4 64 512 728 729 730taylor(x*tanh(x)/(sqrt(1-x^2)-1),x,0,3); 731 732 733 7 2 4 734 - 2 + ---*x + O(x ) 735 6 736 737 738taylor((e^(5*x)-2*x)^(1/x),x,0,2); 739 740 741 3 742 3 3 73*e 2 3 743e + 8*e *x + -------*x + O(x ) 744 3 745 746 747taylor(sin x/cos x,x,pi/2,3); 748 749 750 pi -1 1 pi 1 pi 3 pi 4 751 - (x - ----) + ---*(x - ----) + ----*(x - ----) + O((x - ----) ) 752 2 3 2 45 2 2 753 754 755taylor(log x*sin(x^2)/(x*sinh x),x,0,2); 756 757 758 1 2 3 759log(x)*(1 - ---*x + O(x )) 760 6 761 762 763taylor(1/x-1/sin x,x,0,2); 764 765 766 1 3 767 - ---*x + O(x ) 768 6 769 770 771taylor(tan x/log cos x,x,pi/2,2); 772 773 774 pi -1 pi 775 - (x - ----) + O(x - ----) 776 2 2 777------------------------------- 778 log(pi - 2*x) - log(2) 779 780 781taylor(log(x^2/(x^2-a)),x,0,3); 782 783 784 2 785 - x 786taylor(log(--------),x,0,3) 787 2 788 a - x 789 790 791 792COMMENT Three more complicated examples contributed by Stan Kameny; 793 794 795zz2 := (z*(z-2*pi*i)*(z-pi*i/2)^2)/(sinh z-i); 796 797 798 3 2 2 3 799 z*(2*i*pi - 12*i*pi*z - 9*pi *z + 4*z ) 800zz2 := ------------------------------------------- 801 4*(sinh(z) - i) 802 803dz2 := df(zz2,z); 804 805 806 3 3 2 2 807dz2 := ( - 2*cosh(z)*i*pi *z + 12*cosh(z)*i*pi*z + 9*cosh(z)*pi *z 808 809 4 3 2 810 - 4*cosh(z)*z + 2*sinh(z)*i*pi - 36*sinh(z)*i*pi*z 811 812 2 3 2 3 3 813 - 18*sinh(z)*pi *z + 16*sinh(z)*z + 18*i*pi *z - 16*i*z + 2*pi 814 815 2 2 816 - 36*pi*z )/(4*(sinh(z) - 2*sinh(z)*i - 1)) 817 818z0 := pi*i/2; 819 820 821 i*pi 822z0 := ------ 823 2 824 825taylor(dz2,z,z0,6); 826 827 828 2 829 i*(pi - 16) i*pi pi i*pi 2 830 - 2*pi + --------------*(z - ------) + ----*(z - ------) 831 4 2 2 2 832 833 2 834 i*( - 3*pi + 80) i*pi 3 pi i*pi 4 835 + -------------------*(z - ------) - ----*(z - ------) 836 120 2 24 2 837 838 2 839 i*(5*pi - 168) i*pi 5 i*pi 7 840 + -----------------*(z - ------) + (1 term) + O((z - ------) ) 841 3360 2 2 842 843 844zz3:=(z*(z-2*pi)*(z-pi/2)^2)/(sin z-1); 845 846 847 3 2 2 3 848 z*( - 2*pi + 9*pi *z - 12*pi*z + 4*z ) 849zz3 := ------------------------------------------ 850 4*(sin(z) - 1) 851 852dz3 := df(zz3,z); 853 854 855 3 2 2 3 4 856dz3 := (2*cos(z)*pi *z - 9*cos(z)*pi *z + 12*cos(z)*pi*z - 4*cos(z)*z 857 858 3 2 2 3 859 - 2*sin(z)*pi + 18*sin(z)*pi *z - 36*sin(z)*pi*z + 16*sin(z)*z 860 861 3 2 2 3 2 862 + 2*pi - 18*pi *z + 36*pi*z - 16*z )/(4*(sin(z) - 2*sin(z) + 1)) 863 864z1 := pi/2; 865 866 867 pi 868z1 := ---- 869 2 870 871taylor(dz3,z,z1,6); 872 873 874 2 2 875 pi - 16 pi pi pi 2 3*pi - 80 pi 3 8762*pi + ----------*(z - ----) + ----*(z - ----) + ------------*(z - ----) 877 4 2 2 2 120 2 878 879 2 880 pi pi 4 5*pi - 168 pi 5 pi 7 881 + ----*(z - ----) + -------------*(z - ----) + (1 term) + O((z - ----) ) 882 24 2 3360 2 2 883 884 885taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,6); 886 887 888 5 2 1313 4 2773 6 7 8891 + ---*x + ------*x - -------*x + O(x ) 890 3 1890 11907 891 892 893taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,0); 894 895 8961 + O(x) 897 898 899 900COMMENT If the expansion point is not constant, it has to be taken 901 care of in differentation, as the following examples show; 902 903 904taylor(sin(x+a),x,a,8); 905 906 907 sin(2*a) 2 cos(2*a) 3 908sin(2*a) + cos(2*a)*(x - a) - ----------*(x - a) - ----------*(x - a) 909 2 6 910 911 sin(2*a) 4 cos(2*a) 5 9 912 + ----------*(x - a) + ----------*(x - a) + (3 terms) + O((x - a) ) 913 24 120 914 915df(ws,a); 916 917 918 cos(2*a) 2 sin(2*a) 3 919cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a) + ----------*(x - a) 920 2 6 921 922 cos(2*a) 4 sin(2*a) 5 8 923 + ----------*(x - a) - ----------*(x - a) + (2 terms) + O((x - a) ) 924 24 120 925 926taylor(cos(x+a),x,a,7); 927 928 929 cos(2*a) 2 sin(2*a) 3 930cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a) + ----------*(x - a) 931 2 6 932 933 cos(2*a) 4 sin(2*a) 5 8 934 + ----------*(x - a) - ----------*(x - a) + (2 terms) + O((x - a) ) 935 24 120 936 937 938 939COMMENT A problem are non-analytical terms: rational powers and 940 logarithmic terms can be handled, but other types of essential 941 singularities cannot; 942 943 944taylor(sqrt(x),x,0,2); 945 946 947 1/2 3 948x + O(x ) 949 950 951taylor(asinh(1/x),x,0,5); 952 953 954 1 2 3 4 5 6 7 955 - log(x) + (log(2) + ---*x - ----*x + ----*x + O(x )) 956 4 32 96 957 958 959taylor(e**(1/x),x,0,2); 960 961 962 1/x 963taylor(e ,x,0,2) 964 965 966COMMENT Another example for non-integer powers; 967 968 969sub (y = sqrt (x), yy); 970 971 972 1/2 1 1 3/2 1 2 5/2 9731 + x + ---*x + ---*x + ----*x + O(x ) 974 2 6 24 975 976 977COMMENT Expansion about infinity is possible in principle...; 978 979 980taylor (e**(1/x), x, infinity, 5); 981 982 983 1 1 1 1 1 1 1 1 1 1 9841 + --- + ---*---- + ---*---- + ----*---- + -----*---- + O(----) 985 x 2 2 6 3 24 4 120 5 6 986 x x x x x 987 988xi := taylor (sin (1/x), x, infinity, 5); 989 990 991 1 1 1 1 1 1 992xi := --- - ---*---- + -----*---- + O(----) 993 x 6 3 120 5 6 994 x x x 995 996 997y1 := taylor(x/(x-1), x, infinity, 3); 998 999 1000 1 1 1 1 1001y1 := 1 + --- + ---- + ---- + O(----) 1002 x 2 3 4 1003 x x x 1004 1005z := df(y1, x); 1006 1007 1008 1 1 1 1 1009z := - ---- - 2*---- - 3*---- + O(----) 1010 2 3 4 5 1011 x x x x 1012 1013 1014COMMENT ...but far from being perfect; 1015 1016 1017taylor (1 / sin (x), x, infinity, 5); 1018 1019 1020 1 1021taylor(--------,x,infinity,5) 1022 sin(x) 1023 1024 1025clear z; 1026 1027 1028 1029COMMENT You may access the expansion with the PART operator; 1030 1031 1032part(yy,0); 1033 1034 1035plus 1036 1037part(yy,1); 1038 1039 10401 1041 1042part(yy,4); 1043 1044 1045 3 1046 y 1047---- 1048 6 1049 1050part(yy,6); 1051 1052 1053***** Expression taylor(1 + y + 1/2*y**2 + 1/6*y**3 + 1/24*y**4,y,0,4) 1054does not have part 6 1055 1056 1057 1058COMMENT The template of a Taylor kernel can be extracted; 1059 1060 1061taylortemplate yy; 1062 1063 1064{{y,0,4}} 1065 1066 1067taylortemplate xxa; 1068 1069 1070{{x,1,2}} 1071 1072 1073taylortemplate xi; 1074 1075 1076{{x,infinity,5}} 1077 1078 1079taylortemplate xy; 1080 1081 1082{{x,0,2},{y,0,2}} 1083 1084 1085taylortemplate xx1; 1086 1087 1088{{{x,y},0,2}} 1089 1090 1091COMMENT Here is a slightly less trivial example; 1092 1093 1094exp := (sin (x) * sin (y) / (x * y))**2; 1095 1096 1097 2 2 1098 sin(x) *sin(y) 1099exp := ----------------- 1100 2 2 1101 x *y 1102 1103 1104taylor (exp, x, 0, 1, y, 0, 1); 1105 1106 1107 2 2 11081 + O(x ,y ) 1109 1110taylor (exp, x, 0, 2, y, 0, 2); 1111 1112 1113 1 2 1 2 1 2 2 3 3 11141 - ---*x - ---*y + ---*y *x + O(x ,y ) 1115 3 3 9 1116 1117 1118tt := taylor (exp, {x,y}, 0, 2); 1119 1120 1121 1 2 1 2 3 1122tt := 1 - ---*y - ---*x + O({x,y} ) 1123 3 3 1124 1125 1126COMMENT An example that uses factorization; 1127 1128 1129on factor; 1130 1131 1132 1133ff := y**5 - 1; 1134 1135 1136 4 3 2 1137ff := (y + y + y + y + 1)*(y - 1) 1138 1139 1140zz := sub (y = taylor(e**x, x, 0, 3), ff); 1141 1142 1143 1 2 1 3 4 4 1 2 1 3 4 3 1144zz := ((1 + x + ---*x + ---*x + O(x )) + (1 + x + ---*x + ---*x + O(x )) 1145 2 6 2 6 1146 1147 1 2 1 3 4 2 1 2 1 3 4 1148 + (1 + x + ---*x + ---*x + O(x )) + (1 + x + ---*x + ---*x + O(x )) 1149 2 6 2 6 1150 1151 1 2 1 3 4 1152 + 1)*((1 + x + ---*x + ---*x + O(x )) - 1) 1153 2 6 1154 1155 1156on exp; 1157 1158 1159 1160zz; 1161 1162 1163 1 2 1 3 4 5 1164(1 + x + ---*x + ---*x + O(x )) - 1 1165 2 6 1166 1167 1168COMMENT A simple example of Taylor kernel differentiation; 1169 1170 1171hugo := taylor(e^x,x,0,5); 1172 1173 1174 1 2 1 3 1 4 1 5 6 1175hugo := 1 + x + ---*x + ---*x + ----*x + -----*x + O(x ) 1176 2 6 24 120 1177 1178df(hugo^2,x); 1179 1180 1181 2 8 3 4 4 5 11822 + 4*x + 4*x + ---*x + ---*x + O(x ) 1183 3 3 1184 1185 1186COMMENT The following shows the (limited) capabilities to integrate 1187 Taylor kernels. Only simple cases are supported, otherwise 1188 a warning is printed and the Taylor kernels are converted to 1189 standard representation; 1190 1191 1192zz := taylor (sin x, x, 0, 5); 1193 1194 1195 1 3 1 5 6 1196zz := x - ---*x + -----*x + O(x ) 1197 6 120 1198 1199ww := taylor (cos y, y, 0, 5); 1200 1201 1202 1 2 1 4 6 1203ww := 1 - ---*y + ----*y + O(y ) 1204 2 24 1205 1206 1207int (zz, x); 1208 1209 1210 1 2 1 4 1 6 7 1211---*x - ----*x + -----*x + O(x ) 1212 2 24 720 1213 1214int (ww, x); 1215 1216 1217 x 2 x 4 6 1218x - ---*y + ----*y + O(y ) 1219 2 24 1220 1221int (zz + ww, x); 1222 1223 1224 1 2 1 4 1 6 7 x 2 x 4 6 1225(---*x - ----*x + -----*x + O(x )) + (x - ---*y + ----*y + O(y )) 1226 2 24 720 2 24 1227 1228 1229 1230COMMENT And here we present Taylor series reversion. 1231 We start with the example given by Knuth for the algorithm; 1232 1233 1234taylor (t - t**2, t, 0, 5); 1235 1236 1237 2 6 1238t - t + O(t ) 1239 1240 1241taylorrevert (ws, t, x); 1242 1243 1244 2 3 4 5 6 1245x + x + 2*x + 5*x + 14*x + O(x ) 1246 1247 1248tan!-series := taylor (tan x, x, 0, 5); 1249 1250 1251 1 3 2 5 6 1252tan-series := x + ---*x + ----*x + O(x ) 1253 3 15 1254 1255taylorrevert (tan!-series, x, y); 1256 1257 1258 1 3 1 5 6 1259y - ---*y + ---*y + O(y ) 1260 3 5 1261 1262atan!-series:=taylor (atan y, y, 0, 5); 1263 1264 1265 1 3 1 5 6 1266atan-series := y - ---*y + ---*y + O(y ) 1267 3 5 1268 1269 1270tmp := taylor (e**x, x, 0, 5); 1271 1272 1273 1 2 1 3 1 4 1 5 6 1274tmp := 1 + x + ---*x + ---*x + ----*x + -----*x + O(x ) 1275 2 6 24 120 1276 1277 1278taylorrevert (tmp, x, y); 1279 1280 1281 1 2 1 3 1 4 1 5 6 1282y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 1283 2 3 4 5 1284 1285 1286taylor (log y, y, 1, 5); 1287 1288 1289 1 2 1 3 1 4 1 5 6 1290y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 1291 2 3 4 5 1292 1293 1294 1295COMMENT The following example calculates the perturbation expansion 1296 of the root x = 20 of the following polynomial in terms of 1297 EPS, in ROUNDED mode; 1298 1299 1300poly := for r := 1 : 20 product (x - r); 1301 1302 1303 20 19 18 17 16 15 1304poly := x - 210*x + 20615*x - 1256850*x + 53327946*x - 1672280820*x 1305 1306 14 13 12 1307 + 40171771630*x - 756111184500*x + 11310276995381*x 1308 1309 11 10 9 1310 - 135585182899530*x + 1307535010540395*x - 10142299865511450*x 1311 1312 8 7 6 1313 + 63030812099294896*x - 311333643161390640*x + 1206647803780373360*x 1314 1315 5 4 1316 - 3599979517947607200*x + 8037811822645051776*x 1317 1318 3 2 1319 - 12870931245150988800*x + 13803759753640704000*x 1320 1321 - 8752948036761600000*x + 2432902008176640000 1322 1323 1324on rounded; 1325 1326 1327 1328tpoly := taylor (poly, x, 20, 4); 1329 1330 1331 2 1332tpoly := 1.21649393692e+17*(x - 20) + 4.31564847287e+17*(x - 20) 1333 1334 3 4 1335 + 6.68609351672e+17*(x - 20) + 6.10115975015e+17*(x - 20) 1336 1337 5 1338 + O((x - 20) ) 1339 1340 1341taylorrevert (tpoly, x, eps); 1342 1343 1344 2 3 134520 + 8.22034512178e-18*eps - 2.39726594662e-34*eps + 1.09290580232e-50*eps 1346 1347 4 5 1348 - 5.97114159465e-67*eps + O(eps ) 1349 1350 1351COMMENT Some more examples using rounded mode; 1352 1353 1354precision 13; 1355 1356 135712 1358 1359 1360taylor(sin x/x,x,0,4); 1361 1362 1363 2 4 5 13641 - 0.1666666666667*x + 0.008333333333333*x + O(x ) 1365 1366 1367taylor(sin x,x,pi/2,4); 1368 1369 1370 2 4 13711 - 0.5*(x - 1.570796326795) + 0.04166666666667*(x - 1.570796326795) 1372 1373 5 1374 + O((x - 1.570796326795) ) 1375 1376 1377taylor(tan x,x,pi/2,4); 1378 1379 1380 -1 1381 - (x - 1.570796326795) + 0.3333333333333*(x - 1.570796326795) 1382 1383 3 5 1384 + 0.02222222222222*(x - 1.570796326795) + O((x - 1.570796326795) ) 1385 1386 1387off rounded; 1388 1389 1390 1391 1392COMMENT An example that involves computing limits of type 0/0 if 1393 expansion is done via differentiation; 1394 1395 1396 1397taylor(sqrt((e^x - 1)/x),x,0,15); 1398 1399 1400 1 5 2 1 3 79 4 3 5 16 14011 + ---*x + ----*x + -----*x + -------*x + -------*x + (10 terms) + O(x ) 1402 4 96 128 92160 40960 1403 1404 1405COMMENT An example that involves intermediate non-analytical terms 1406 which cancel entirely; 1407 1408 1409taylor(x^(5/2)/(log(x+1)*tan(x^(3/2))),x,0,5); 1410 1411 1412 1 1 2 7 3 139 4 67 5 6 14131 + ---*x - ----*x - ----*x - -----*x + ------*x + O(x ) 1414 2 12 24 720 1440 1415 1416 1417COMMENT Other examples involving non-analytical terms; 1418 1419 1420taylor(log(e^x-1),x,0,5); 1421 1422 1423 1 1 2 1 4 5 1424log(x) + (---*x + ----*x - ------*x + O(x )) 1425 2 24 2880 1426 1427 1428taylor(e^(1/x)*(e^x-1),x,0,5); 1429 1430 1431 1/x 1 2 1 3 1 4 1 5 6 1432e *(x + ---*x + ---*x + ----*x + -----*x + O(x )) 1433 2 6 24 120 1434 1435 1436taylor(log(x)*x^10,x,0,5); 1437 1438 1439 10 11 1440log(x)*(x + O(x )) 1441 1442 1443taylor(log(x)*x^10,x,0,11); 1444 1445 1446 10 12 1447log(x)*(x + O(x )) 1448 1449 1450taylor(log(x-a)/((a-b)*(a-c)) + log(2(x-b))/((b-c)*(b-a)) 1451 + log(x-c)/((c-a)*(c-b)),x,infinity,2); 1452 1453 1454 log(2) 1 1 1 1455 - ---------------------- - ---*---- + O(----) 1456 2 2 2 3 1457 a*b - a*c - b + b*c x x 1458 1459 1460ss := (sqrt(x^(2/5) +1) - x^(1/3)-1)/x^(1/3); 1461 1462 1463 2/5 1/3 1464 sqrt(x + 1) - x - 1 1465ss := --------------------------- 1466 1/3 1467 x 1468 1469taylor(exp ss,x,0,2); 1470 1471 1472 1 1 1/15 1 2/15 1 1/5 1 4/15 1 1/3 1473--- + -----*x + -----*x + ------*x + -------*x + --------*x 1474 e 2*e 8*e 48*e 384*e 3840*e 1475 1476 31/15 1477 + (25 terms) + O(x ) 1478 1479 1480taylor(exp sub(x=x^15,ss),x,0,2); 1481 1482 1483 1 1 1 2 3 1484--- + -----*x + -----*x + O(x ) 1485 e 2*e 8*e 1486 1487 1488taylor(dilog(x),x,0,4); 1489 1490 1491 1 2 1 3 1 4 5 1492log(x)*(x + ---*x + ---*x + ---*x + O(x )) 1493 2 3 4 1494 1495 2 1496 pi 1 2 1 3 1 4 5 1497 + (----- - x - ---*x - ---*x - ----*x + O(x )) 1498 6 4 9 16 1499 1500 1501taylor(dilog(x),x,1,4); 1502 1503 1504 1 2 1 3 1 4 5 1505 - (x - 1) + ---*(x - 1) - ---*(x - 1) + ----*(x - 1) + O((x - 1) ) 1506 4 9 16 1507 1508 1509taylor(dilog(x),x,-1,4); 1510 1511 1512 pi*( - 4*log(2)*i + pi) log(-1) log(-1) - 2 2 1513------------------------- + ---------*(x + 1) + -------------*(x + 1) 1514 4 2 8 1515 1516 log(-1) - 4 3 3*log(-1) - 20 4 5 1517 + -------------*(x + 1) + ----------------*(x + 1) + O((x + 1) ) 1518 24 192 1519 1520 1521taylor(Ei(x),x,0,4); 1522 1523 1524 1 2 1 3 1 4 5 1525log(x) + (x + ---*x + ----*x + ----*x + O(x )) + euler_gamma 1526 4 18 96 1527 1528 1529taylor(Ei(x),x,1,4); 1530 1531 1532 e 3 e 4 5 1533ei(1) + e*(x - 1) + ---*(x - 1) - ----*(x - 1) + O((x - 1) ) 1534 6 12 1535 1536 1537taylor(Ei(x),x,-1,4); 1538 1539 1540 1 1 2 5 3 2 4 1541ei(-1) - ---*(x + 1) - ---*(x + 1) - -----*(x + 1) - -----*(x + 1) 1542 e e 6*e 3*e 1543 1544 5 1545 + O((x + 1) ) 1546 1547 1548COMMENT In the following we demonstrate the possibiblity to compute the 1549 expansion of a function which is given by a simple first order 1550 differential equation: the function myexp(x) is exp(-x^2); 1551 1552 1553operator myexp,myerf; 1554 1555 1556let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 1557 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 1558 1559 1560taylor(myexp(x),x,0,5); 1561 1562 1563 2 1 4 6 15641 - x + ---*x + O(x ) 1565 2 1566 1567taylor(myerf(x),x,0,5); 1568 1569 1570 2 2 3 1 5 6 1571----------*x - ------------*x + ------------*x + O(x ) 1572 sqrt(pi) 3*sqrt(pi) 5*sqrt(pi) 1573 1574clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 1575 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 1576 1577 1578clear myexp,myerf; 1579 1580 1581 1582%%% showtime; 1583 1584COMMENT There are two special operators, implicit_taylor and 1585 inverse_taylor, to compute the Taylor expansion of implicit 1586 or inverse functions; 1587 1588 1589 1590implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5); 1591 1592 1593 1 2 1 4 6 15941 - ---*x - ---*x + O(x ) 1595 2 8 1596 1597 1598implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20); 1599 1600 1601 1 2 1 4 1 6 5 8 7 10 21 16021 - ---*x - ---*x - ----*x - -----*x - -----*x + (5 terms) + O(x ) 1603 2 8 16 128 256 1604 1605 1606implicit_taylor(x+y^3-y,x,y,0,0,8); 1607 1608 1609 3 5 7 9 1610x + x + 3*x + 12*x + O(x ) 1611 1612 1613implicit_taylor(x+y^3-y,x,y,0,1,5); 1614 1615 1616 1 3 2 1 3 105 4 3 5 6 16171 - ---*x - ---*x - ---*x - -----*x - ---*x + O(x ) 1618 2 8 2 128 2 1619 1620 1621implicit_taylor(x+y^3-y,x,y,0,-1,5); 1622 1623 1624 1 3 2 1 3 105 4 3 5 6 1625 - 1 - ---*x + ---*x - ---*x + -----*x - ---*x + O(x ) 1626 2 8 2 128 2 1627 1628 1629implicit_taylor(y*e^y-x,x,y,0,0,5); 1630 1631 1632 2 3 3 8 4 125 5 6 1633x - x + ---*x - ---*x + -----*x + O(x ) 1634 2 3 24 1635 1636 1637COMMENT This is the function exp(-1/x^2), which has an essential 1638 singularity at the point 0; 1639 1640 1641implicit_taylor(x^2*log y+1,x,y,0,0,3); 1642 1643 1644***** log 0 formed 1645 1646***** Computation of Taylor series of implicit function failed 1647 1648 1649inverse_taylor(exp(x)-1,x,y,0,8); 1650 1651 1652 1 2 1 3 1 4 1 5 1 6 9 1653y - ---*y + ---*y - ---*y + ---*y - ---*y + (2 terms) + O(y ) 1654 2 3 4 5 6 1655 1656 1657inverse_taylor(exp(x),x,y,0,5); 1658 1659 1660 1 2 1 3 1 4 1 5 6 1661y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 1662 2 3 4 5 1663 1664 1665inverse_taylor(sqrt(x),x,y,0,5); 1666 1667 1668 2 6 1669y + O(y ) 1670 1671 1672inverse_taylor(log(1+x),x,y,0,5); 1673 1674 1675 1 2 1 3 1 4 1 5 6 1676y + ---*y + ---*y + ----*y + -----*y + O(y ) 1677 2 6 24 120 1678 1679 1680inverse_taylor((e^x-e^(-x))/2,x,y,0,5); 1681 1682 1683 1 3 3 5 6 1684y - ---*y + ----*y + O(y ) 1685 6 40 1686 1687 1688COMMENT In the next two cases the inverse functions have a branch 1689 point, therefore the computation fails; 1690 1691 1692inverse_taylor((e^x+e^(-x))/2,x,y,0,5); 1693 1694 1695***** Computation of Taylor series of inverse function failed 1696 1697 1698inverse_taylor(exp(x^2-1),x,y,0,5); 1699 1700 1701***** Computation of Taylor series of inverse function failed 1702 1703 1704inverse_taylor(exp(sqrt(x))-1,x,y,0,5); 1705 1706 1707 2 3 11 4 5 5 6 1708y - y + ----*y - ---*y + O(y ) 1709 12 6 1710 1711 1712inverse_taylor(x*exp(x),x,y,0,5); 1713 1714 1715 2 3 3 8 4 125 5 6 1716y - y + ---*y - ---*y + -----*y + O(y ) 1717 2 3 24 1718 1719 1720 1721%%% showtime; 1722 1723 1724COMMENT An application is the problem posed by Prof. Stanley: 1725 we prove that the finite difference expression below 1726 corresponds to the given derivative expression; 1727 1728 1729operator diff,a,f,gg; 1730 1731 % We use gg to avoid conflict with high energy 1732 % physics operator. 1733 1734let diff(~f,~arg) => df(f,arg); 1735 1736 1737 1738derivative_expression := 1739diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) + 1740diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ; 1741 1742 1743derivative_expression := 2*a(x,y)*df(f(x,y),x,y)*df(gg(x,y),x)*df(gg(x,y),y) 1744 1745 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x,y)*df(gg(x,y),y) 1746 1747 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y,2) 1748 1749 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,y)*df(gg(x,y),x) 1750 1751 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,2)*df(gg(x,y),y) 1752 1753 + df(a(x,y),x)*df(f(x,y),y)*df(gg(x,y),x)*df(gg(x,y),y) 1754 1755 + df(a(x,y),y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y) 1756 1757 1758finite_difference_expression := 1759+a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1760+a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1761+a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1762+a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1763-f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1764-f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1765-f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1766-a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1767-gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1768-gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1769-gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1770-a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1771+f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1772+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1773+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1774+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1775-gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1776+gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1777-gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1778+gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1779-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1780-a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1781-a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1782+gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1783+a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1784+a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1785-gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1786+gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1787-a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1788-a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1789+gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1790+a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1791+f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1792-f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2) 1793+f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1794-f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1795+a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1796+a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1797+a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1798+a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1799-f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1800-f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1801-f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1802-a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1803-gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1804-gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1805-gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1806-a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1807+f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1808+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1809+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1810+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1811-gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1812+gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1813-gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1814+gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1815-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1816-a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1817-a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1818+gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1819+a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1820+a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1821-gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1822+gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1823-a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1824-a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1825+gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1826+a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1827+f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1828-f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2) 1829+f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1830-f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1831+f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 1832+f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 1833+f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 1834+a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 1835-f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 1836-f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 1837-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1838-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1839-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1840-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1841+f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 1842+f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 1843-f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2) 1844+a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1845+a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1846+a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1847+a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1848-f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1849-f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1850-f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1851-a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1852-gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1853-gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1854-gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1855-a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1856+f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1857+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1858+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1859+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1860-gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1861+gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1862-gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1863+gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1864-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1865-a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1866-a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1867+gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1868+a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1869+a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1870-gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1871+gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1872-a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1873-a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1874+gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1875+a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1876+f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1877-f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2) 1878+f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1879-f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1880+a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1881+a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1882+a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1883+a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1884-f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1885-f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1886-f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1887-a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1888-gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1889-gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1890-gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1891-a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1892+f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1893+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1894+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1895+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1896-gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1897+gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1898-gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1899+gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1900-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1901-a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1902-a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1903+gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1904+a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1905+a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1906-gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1907+gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1908-a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1909-a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1910+gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1911+a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1912+f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1913-f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2) 1914+f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1915-f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1916+f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 1917+f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 1918+f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 1919+a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 1920-f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 1921-f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 1922-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1923-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1924-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1925-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1926+f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 1927+f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 1928-f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2) 1929+f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2) 1930+a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2) 1931-f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2) 1932+f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2) 1933+a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2) 1934-f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2) 1935-a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$ 1936 1937 1938 1939COMMENT We define abbreviations for the partial derivatives; 1940 1941 1942operator ax,ay,fx,fy,gx,gy; 1943 1944 1945operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 1946 1947 1948operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 1949 1950 1951operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy, 1952 gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 1953 1954 1955operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 1956 1957 1958operator gxxxxyy,gxxxyyy,gxxyyyy; 1959 1960 1961 1962operator_diff_rules := { 1963 df(a(~x,~y),~x) => ax(x,y), 1964 df(a(~x,~y),~y) => ay(x,y), 1965 df(f(~x,~y),~x) => fx(x,y), 1966 df(f(~x,~y),~y) => fy(x,y), 1967 df(gg(~x,~y),~x) => gx(x,y), 1968 df(gg(~x,~y),~y) => gy(x,y), 1969 1970 df(ax(~x,~y),~x) => axx(x,y), 1971 df(ax(~x,~y),~y) => axy(x,y), 1972 df(ay(~x,~y),~x) => axy(x,y), 1973 df(ay(~x,~y),~y) => ayy(x,y), 1974 df(fx(~x,~y),~x) => fxx(x,y), 1975 df(fx(~x,~y),~y) => fxy(x,y), 1976 df(fy(~x,~y),~x) => fxy(x,y), 1977 df(fy(~x,~y),~y) => fyy(x,y), 1978 df(gx(~x,~y),~x) => gxx(x,y), 1979 df(gx(~x,~y),~y) => gxy(x,y), 1980 df(gy(~x,~y),~x) => gxy(x,y), 1981 df(gy(~x,~y),~y) => gyy(x,y), 1982 1983 df(axx(~x,~y),~x) => axxx(x,y), 1984 df(axy(~x,~y),~x) => axxy(x,y), 1985 df(ayy(~x,~y),~x) => axyy(x,y), 1986 df(ayy(~x,~y),~y) => ayyy(x,y), 1987 df(fxx(~x,~y),~x) => fxxx(x,y), 1988 df(fxy(~x,~y),~x) => fxxy(x,y), 1989 df(fxy(~x,~y),~y) => fxyy(x,y), 1990 df(fyy(~x,~y),~x) => fxyy(x,y), 1991 df(fyy(~x,~y),~y) => fyyy(x,y), 1992 df(gxx(~x,~y),~x) => gxxx(x,y), 1993 df(gxx(~x,~y),~y) => gxxy(x,y), 1994 df(gxy(~x,~y),~x) => gxxy(x,y), 1995 df(gxy(~x,~y),~y) => gxyy(x,y), 1996 df(gyy(~x,~y),~x) => gxyy(x,y), 1997 df(gyy(~x,~y),~y) => gyyy(x,y), 1998 1999 df(axyy(~x,~y),~x) => axxyy(x,y), 2000 df(axxy(~x,~y),~x) => axxxy(x,y), 2001 df(ayyy(~x,~y),~x) => axyyy(x,y), 2002 df(fxxy(~x,~y),~x) => fxxxy(x,y), 2003 df(fxyy(~x,~y),~x) => fxxyy(x,y), 2004 df(fyyy(~x,~y),~x) => fxyyy(x,y), 2005 df(gxxx(~x,~y),~x) => gxxxx(x,y), 2006 df(gxxy(~x,~y),~x) => gxxxy(x,y), 2007 df(gxyy(~x,~y),~x) => gxxyy(x,y), 2008 df(gyyy(~x,~y),~x) => gxyyy(x,y), 2009 df(gyyy(~x,~y),~y) => gyyyy(x,y), 2010 2011 df(axxyy(~x,~y),~x) => axxxyy(x,y), 2012 df(axyyy(~x,~y),~x) => axxyyy(x,y), 2013 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 2014 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 2015 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 2016 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 2017 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 2018 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 2019 2020 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 2021 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 2022 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y) 2023}; 2024 2025 2026operator_diff_rules := {df(a(~x,~y),~x) => ax(x,y), 2027 2028 df(a(~x,~y),~y) => ay(x,y), 2029 2030 df(f(~x,~y),~x) => fx(x,y), 2031 2032 df(f(~x,~y),~y) => fy(x,y), 2033 2034 df(gg(~x,~y),~x) => gx(x,y), 2035 2036 df(gg(~x,~y),~y) => gy(x,y), 2037 2038 df(ax(~x,~y),~x) => axx(x,y), 2039 2040 df(ax(~x,~y),~y) => axy(x,y), 2041 2042 df(ay(~x,~y),~x) => axy(x,y), 2043 2044 df(ay(~x,~y),~y) => ayy(x,y), 2045 2046 df(fx(~x,~y),~x) => fxx(x,y), 2047 2048 df(fx(~x,~y),~y) => fxy(x,y), 2049 2050 df(fy(~x,~y),~x) => fxy(x,y), 2051 2052 df(fy(~x,~y),~y) => fyy(x,y), 2053 2054 df(gx(~x,~y),~x) => gxx(x,y), 2055 2056 df(gx(~x,~y),~y) => gxy(x,y), 2057 2058 df(gy(~x,~y),~x) => gxy(x,y), 2059 2060 df(gy(~x,~y),~y) => gyy(x,y), 2061 2062 df(axx(~x,~y),~x) => axxx(x,y), 2063 2064 df(axy(~x,~y),~x) => axxy(x,y), 2065 2066 df(ayy(~x,~y),~x) => axyy(x,y), 2067 2068 df(ayy(~x,~y),~y) => ayyy(x,y), 2069 2070 df(fxx(~x,~y),~x) => fxxx(x,y), 2071 2072 df(fxy(~x,~y),~x) => fxxy(x,y), 2073 2074 df(fxy(~x,~y),~y) => fxyy(x,y), 2075 2076 df(fyy(~x,~y),~x) => fxyy(x,y), 2077 2078 df(fyy(~x,~y),~y) => fyyy(x,y), 2079 2080 df(gxx(~x,~y),~x) => gxxx(x,y), 2081 2082 df(gxx(~x,~y),~y) => gxxy(x,y), 2083 2084 df(gxy(~x,~y),~x) => gxxy(x,y), 2085 2086 df(gxy(~x,~y),~y) => gxyy(x,y), 2087 2088 df(gyy(~x,~y),~x) => gxyy(x,y), 2089 2090 df(gyy(~x,~y),~y) => gyyy(x,y), 2091 2092 df(axyy(~x,~y),~x) => axxyy(x,y), 2093 2094 df(axxy(~x,~y),~x) => axxxy(x,y), 2095 2096 df(ayyy(~x,~y),~x) => axyyy(x,y), 2097 2098 df(fxxy(~x,~y),~x) => fxxxy(x,y), 2099 2100 df(fxyy(~x,~y),~x) => fxxyy(x,y), 2101 2102 df(fyyy(~x,~y),~x) => fxyyy(x,y), 2103 2104 df(gxxx(~x,~y),~x) => gxxxx(x,y), 2105 2106 df(gxxy(~x,~y),~x) => gxxxy(x,y), 2107 2108 df(gxyy(~x,~y),~x) => gxxyy(x,y), 2109 2110 df(gyyy(~x,~y),~x) => gxyyy(x,y), 2111 2112 df(gyyy(~x,~y),~y) => gyyyy(x,y), 2113 2114 df(axxyy(~x,~y),~x) => axxxyy(x,y), 2115 2116 df(axyyy(~x,~y),~x) => axxyyy(x,y), 2117 2118 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 2119 2120 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 2121 2122 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 2123 2124 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 2125 2126 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 2127 2128 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 2129 2130 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 2131 2132 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 2133 2134 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)} 2135 2136 2137let operator_diff_rules; 2138 2139 2140 2141texp := taylor (finite_difference_expression, dx, 0, 1, dy, 0, 1); 2142 2143 2144texp := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y) 2145 2146 + 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) 2147 2148 + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y) 2149 2150 2 2 2151 + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y) + O(dx ,dy ) 2152 2153 2154COMMENT You may also try to expand further but this needs a lot 2155 of CPU time. Therefore the following line is commented out; 2156 2157 2158%texp := taylor (finite_difference_expression, dx, 0, 2, dy, 0, 2); 2159 2160factor dx,dy; 2161 2162 2163 2164result := taylortostandard texp; 2165 2166 2167result := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y) 2168 2169 + 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) 2170 2171 + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y) 2172 2173 + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y) 2174 2175 2176derivative_expression - result; 2177 2178 21790 2180 2181 2182 2183clear diff(~f,~arg); 2184 2185 2186clearrules operator_diff_rules; 2187 2188 2189clear diff,a,f,gg; 2190 2191 2192clear ax,ay,fx,fy,gx,gy; 2193 2194 2195clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 2196 2197 2198clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 2199 2200 2201clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 2202 2203 2204clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 2205 2206 2207clear gxxxxyy,gxxxyyy,gxxyyyy; 2208 2209 2210 2211taylorprintterms := 5; 2212 2213 2214taylorprintterms := 5 2215 2216 2217off taylorautoexpand,taylorkeeporiginal; 2218 2219 2220 2221%%% showtime; 2222 2223COMMENT An example provided by Alan Barnes: elliptic functions; 2224 2225 2226% Jacobi's elliptic functions 2227% sn(x,k), cn(x,k), dn(x,k). 2228% The modulus and complementary modulus are denoted by K and K!' 2229% usually written mathematically as k and k' respectively 2230% 2231% epsilon(x,k) is the incomplete elliptic integral of the second kind 2232% usually written mathematically as E(x,k) 2233% 2234% KK(k) is the complete elliptic integral of the first kind 2235% usually written mathematically as K(k) 2236% K(k) = arcsn(1,k) 2237% KK!'(k) is the complementary complete integral 2238% usually written mathematically as K'(k) 2239% NB. K'(k) = K(k') 2240% EE(k) is the complete elliptic integral of the second kind 2241% usually written mathematically as E(k) 2242% EE!'(k) is the complementary complete integral 2243% usually written mathematically as E'(k) 2244% NB. E'(k) = E(k') 2245 2246operator sn, cn, dn, epsilon; 2247 2248 2249 2250elliptic_rules := { 2251% Differentiation rules for basic functions 2252 df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 2253 df(cn(~x,~k),~x) => -sn(x,k)*dn(x,k), 2254 df(dn(~x,~k),~x) => -k^2*sn(x,k)*cn(x,k), 2255 df(epsilon(~x,~k),~x)=> dn(x,k)^2, 2256 2257% k-derivatives 2258% DF Lawden Elliptic Functions & Applications Springer (1989) 2259 df(sn(~x,~k),~k) => (k*sn(x,k)*cn(x,k)^2 2260 -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2) 2261 + x*cn(x,k)*dn(x,k)/k, 2262 df(cn(~x,~k),~k) => (-k*sn(x,k)^2*cn(x,k) 2263 +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2) 2264 - x*sn(x,k)*dn(x,k)/k, 2265 df(dn(~x,~k),~k) => k*(-sn(x,k)^2*dn(x,k) 2266 +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2) 2267 - k*x*sn(x,k)*cn(x,k), 2268 df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k) 2269 -epsilon(x,k)*cn(x,k)^2)/(1-k^2) 2270 -k*x*sn(x,k)^2, 2271 2272% parity properties 2273 sn(-~x,~k) => -sn(x,k), 2274 cn(-~x,~k) => cn(x,k), 2275 dn(-~x,~k) => dn(x,k), 2276 epsilon(-~x,~k) => -epsilon(x,k), 2277 sn(~x,-~k) => sn(x,k), 2278 cn(~x,-~k) => cn(x,k), 2279 dn(~x,-~k) => dn(x,k), 2280 epsilon(~x,-~k) => epsilon(x,k), 2281 2282 2283% behaviour at zero 2284 sn(0,~k) => 0, 2285 cn(0,~k) => 1, 2286 dn(0,~k) => 1, 2287 epsilon(0,~k) => 0, 2288 2289 2290% degenerate cases of modulus 2291 sn(~x,0) => sin(x), 2292 cn(~x,0) => cos(x), 2293 dn(~x,0) => 1, 2294 epsilon(~x,0) => x, 2295 2296 sn(~x,1) => tanh(x), 2297 cn(~x,1) => 1/cosh(x), 2298 dn(~x,1) => 1/cosh(x), 2299 epsilon(~x,1) => tanh(x) 2300}; 2301 2302 2303elliptic_rules := {df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 2304 2305 df(cn(~x,~k),~x) => - sn(x,k)*dn(x,k), 2306 2307 2 2308 df(dn(~x,~k),~x) => - k *sn(x,k)*cn(x,k), 2309 2310 2 2311 df(epsilon(~x,~k),~x) => dn(x,k) , 2312 2313 2 dn(x,k) 2314 k*sn(x,k)*cn(x,k) - epsilon(x,k)*cn(x,k)*--------- 2315 k 2316 df(sn(~x,~k),~k) => ----------------------------------------------------- 2317 2 2318 1 - k 2319 2320 dn(x,k) 2321 + x*cn(x,k)*---------, 2322 k 2323 2324 2 dn(x,k) 2325 - k*sn(x,k) *cn(x,k) + epsilon(x,k)*sn(x,k)*--------- 2326 k 2327 df(cn(~x,~k),~k) => -------------------------------------------------------- 2328 2 2329 1 - k 2330 2331 dn(x,k) 2332 - x*sn(x,k)*---------, 2333 k 2334 2335 2 2336 - sn(x,k) *dn(x,k) + epsilon(x,k)*sn(x,k)*cn(x,k) 2337 df(dn(~x,~k),~k) => k*---------------------------------------------------- 2338 2 2339 1 - k 2340 2341 - k*x*sn(x,k)*cn(x,k), 2342 2343 df(epsilon(~x,~k),~k) 2344 2345 2 2346 sn(x,k)*cn(x,k)*dn(x,k) - epsilon(x,k)*cn(x,k) 2 2347 => k*------------------------------------------------- - k*x*sn(x,k) , 2348 2 2349 1 - k 2350 2351 sn( - ~x,~k) => - sn(x,k), 2352 2353 cn( - ~x,~k) => cn(x,k), 2354 2355 dn( - ~x,~k) => dn(x,k), 2356 2357 epsilon( - ~x,~k) => - epsilon(x,k), 2358 2359 sn(~x, - ~k) => sn(x,k), 2360 2361 cn(~x, - ~k) => cn(x,k), 2362 2363 dn(~x, - ~k) => dn(x,k), 2364 2365 epsilon(~x, - ~k) => epsilon(x,k), 2366 2367 sn(0,~k) => 0, 2368 2369 cn(0,~k) => 1, 2370 2371 dn(0,~k) => 1, 2372 2373 epsilon(0,~k) => 0, 2374 2375 sn(~x,0) => sin(x), 2376 2377 cn(~x,0) => cos(x), 2378 2379 dn(~x,0) => 1, 2380 2381 epsilon(~x,0) => x, 2382 2383 sn(~x,1) => tanh(x), 2384 2385 1 2386 cn(~x,1) => ---------, 2387 cosh(x) 2388 2389 1 2390 dn(~x,1) => ---------, 2391 cosh(x) 2392 2393 epsilon(~x,1) => tanh(x)} 2394 2395 2396let elliptic_rules; 2397 2398 2399 2400hugo := taylor(sn(x,k),k,0,6); 2401 2402 2403 2 2 2404 cos(x)*(cos(x) *x + cos(x)*sin(x) + sin(x) *x - 2*x) 2 2405hugo := sin(x) + ------------------------------------------------------*k + ( 2406 4 2407 2408 5 4 2 4 2409 cos(x) *x - 2*cos(x) *sin(x)*x + 5*cos(x) *sin(x) 2410 2411 3 2 3 2 3 2 2412 - 10*cos(x) *sin(x) *x + 6*cos(x) *x - 4*cos(x) *sin(x) *x 2413 2414 2 3 2 2 2 2415 + cos(x) *sin(x) + 8*cos(x) *sin(x)*x + 4*cos(x) *sin(x) 2416 2417 4 2 2418 - 11*cos(x)*sin(x) *x + 30*cos(x)*sin(x) *x - 16*cos(x)*x 2419 2420 5 2 3 2 2 4 2421 - 2*sin(x) *x + 8*sin(x) *x - 8*sin(x)*x )/64*k + ( 2422 2423 7 3 7 6 2 2424 - 6*cos(x) *x + 17*cos(x) *x - 99*cos(x) *sin(x)*x 2425 2426 6 5 2 3 5 2 2427 + 21*cos(x) *sin(x) - 18*cos(x) *sin(x) *x - 71*cos(x) *sin(x) *x 2428 2429 5 3 5 4 3 2 2430 + 36*cos(x) *x - 18*cos(x) *x - 135*cos(x) *sin(x) *x 2431 2432 4 3 4 2 4 2433 - 133*cos(x) *sin(x) + 324*cos(x) *sin(x)*x + 172*cos(x) *sin(x) 2434 2435 3 4 3 3 4 2436 - 18*cos(x) *sin(x) *x - 13*cos(x) *sin(x) *x 2437 2438 3 2 3 3 2 3 3 2439 + 72*cos(x) *sin(x) *x - 156*cos(x) *sin(x) *x - 72*cos(x) *x 2440 2441 3 2 5 2 2 5 2442 + 160*cos(x) *x + 27*cos(x) *sin(x) *x - 118*cos(x) *sin(x) 2443 2444 2 3 2 2 2 2445 + 176*cos(x) *sin(x) - 108*cos(x) *sin(x)*x + 32*cos(x) *sin(x) 2446 2447 6 3 6 4 3 2448 - 6*cos(x)*sin(x) *x + 75*cos(x)*sin(x) *x + 36*cos(x)*sin(x) *x 2449 2450 4 2 3 2 2451 - 498*cos(x)*sin(x) *x - 72*cos(x)*sin(x) *x + 888*cos(x)*sin(x) *x 2452 2453 3 7 2 5 2 2454 + 48*cos(x)*x - 384*cos(x)*x + 63*sin(x) *x - 324*sin(x) *x 2455 2456 3 2 2 6 7 2457 + 540*sin(x) *x - 288*sin(x)*x )/2304*k + O(k ) 2458 2459otto := taylor(cn(x,k),k,0,6); 2460 2461 2462 2 2 2463 sin(x)*( - cos(x) *x - cos(x)*sin(x) - sin(x) *x + 2*x) 2 2464otto := cos(x) + ---------------------------------------------------------*k + 2465 4 2466 2467 5 2 4 3 2 2 2468 ( - 2*cos(x) *x - 5*cos(x) *sin(x)*x - 4*cos(x) *sin(x) *x 2469 2470 3 2 3 2 2 3 2471 - 7*cos(x) *sin(x) + 8*cos(x) *x + 2*cos(x) *sin(x) *x 2472 2473 2 4 2 4 2474 + 2*cos(x) *sin(x)*x - 2*cos(x)*sin(x) *x - 3*cos(x)*sin(x) 2475 2476 2 2 2 2 5 2477 + 8*cos(x)*sin(x) *x - 4*cos(x)*sin(x) - 8*cos(x)*x + 7*sin(x) *x 2478 2479 3 4 7 2 2480 - 22*sin(x) *x + 16*sin(x)*x)/64*k + ( - 9*cos(x) *x 2481 2482 6 3 6 5 2 2 2483 + 6*cos(x) *sin(x)*x - 71*cos(x) *sin(x)*x + 135*cos(x) *sin(x) *x 2484 2485 5 2 5 2 4 3 3 2486 - 66*cos(x) *sin(x) - 36*cos(x) *x + 18*cos(x) *sin(x) *x 2487 2488 4 3 4 3 4 2489 - cos(x) *sin(x) *x - 36*cos(x) *sin(x)*x + 18*cos(x) *sin(x)*x 2490 2491 3 4 2 3 4 2492 + 297*cos(x) *sin(x) *x + 61*cos(x) *sin(x) 2493 2494 3 2 2 3 2 3 2 2495 - 720*cos(x) *sin(x) *x - 208*cos(x) *sin(x) + 252*cos(x) *x 2496 2497 2 5 3 2 5 2498 + 18*cos(x) *sin(x) *x + 31*cos(x) *sin(x) *x 2499 2500 2 3 3 2 3 2501 - 72*cos(x) *sin(x) *x - 24*cos(x) *sin(x) *x 2502 2503 2 3 2 6 2 2504 + 72*cos(x) *sin(x)*x + 56*cos(x) *sin(x)*x + 153*cos(x)*sin(x) *x 2505 2506 6 4 2 4 2507 + 91*cos(x)*sin(x) - 684*cos(x)*sin(x) *x - 212*cos(x)*sin(x) 2508 2509 2 2 2 2 2510 + 900*cos(x)*sin(x) *x - 32*cos(x)*sin(x) - 288*cos(x)*x 2511 2512 7 3 7 5 3 5 2513 + 6*sin(x) *x - 39*sin(x) *x - 36*sin(x) *x + 318*sin(x) *x 2514 2515 3 3 3 3 2516 + 72*sin(x) *x - 672*sin(x) *x - 48*sin(x)*x + 384*sin(x)*x)/2304 2517 2518 6 7 2519 *k + O(k ) 2520 2521taylorcombine(hugo^2 + otto^2); 2522 2523 2524 2 2 7 2525cos(x) + sin(x) + O(k ) 2526 2527 2528clearrules elliptic_rules; 2529 2530 2531 2532clear sn, cn, dn, epsilon; 2533 2534 2535 2536COMMENT Examples of gamma function and its derivatives; 2537 2538 2539taylor(gamma(x),x,1,3); 2540 2541 2542 2 2 2543 6*euler_gamma + pi 2 25441 - euler_gamma*(x - 1) + ----------------------*(x - 1) 2545 12 2546 2547 3 2 2548 - 4*zeta(3) - 2*euler_gamma - euler_gamma*pi 3 4 2549 + -------------------------------------------------*(x - 1) + O((x - 1) ) 2550 12 2551 2552 2553taylor(gamma(x),x,1/2,3); 2554 2555 2556 1 2557sqrt(pi) + sqrt(pi)*( - 2*log(2) - euler_gamma)*(x - ---) + 2558 2 2559 2560 2 1 2 2561 sqrt(pi)*(4*log(2) + 4*log(2)*euler_gamma + polygamma(1,---) + euler_gamma ) 2562 2 2563------------------------------------------------------------------------------- 2564 2 2565 2566 1 2 3 2 2567*(x - ---) + (sqrt(pi)*( - 8*log(2) - 12*log(2) *euler_gamma 2568 2 2569 2570 1 2 2571 - 6*log(2)*polygamma(1,---) - 6*log(2)*euler_gamma 2572 2 2573 2574 1 1 2575 + polygamma(2,---) - 3*polygamma(1,---)*euler_gamma 2576 2 2 2577 2578 3 1 3 1 4 2579 - euler_gamma ))/6*(x - ---) + O((x - ---) ) 2580 2 2 2581 2582 2583taylor(gamma(x),x,a,3); 2584 2585 2586gamma(a) + gamma(a)*psi(a)*(x - a) 2587 2588 2 2589 gamma(a)*(polygamma(1,a) + psi(a) ) 2 2590 + -------------------------------------*(x - a) 2591 2 2592 2593 3 2594 gamma(a)*(polygamma(2,a) + 3*polygamma(1,a)*psi(a) + psi(a) ) 3 2595 + ---------------------------------------------------------------*(x - a) 2596 6 2597 2598 4 2599 + O((x - a) ) 2600 2601 2602taylor(gamma(x),x,-1,3); 2603 2604 2605***** Psi undefined for nonpositive integer argument -1 2606 2607***** Zero divisor in computation of constant term 2608 2609***** Psi undefined for nonpositive integer argument -1 2610 2611***** Error during expansion (possible singularity!) 2612 2613***** Psi undefined for nonpositive integer argument -1 2614 2615***** Error during expansion (possible singularity!) 2616 2617***** Psi undefined for nonpositive integer argument -1 2618 2619***** Error during expansion (possible singularity!) 2620 2621***** Psi undefined for nonpositive integer argument -1 2622 2623***** Error during expansion (possible singularity!) 2624 2625 2626taylor(gamma(1+x),x,0, 6); 2627 2628 2629 2 2 2630 6*euler_gamma + pi 2 26311 - euler_gamma*x + ----------------------*x 2632 12 2633 2634 3 2 2635 - 4*zeta(3) - 2*euler_gamma - euler_gamma*pi 3 2636 + -------------------------------------------------*x 2637 12 2638 2639 4 2 2 4 2640 160*zeta(3)*euler_gamma + 20*euler_gamma + 20*euler_gamma *pi + 3*pi 4 2641 + -------------------------------------------------------------------------*x 2642 480 2643 2644 7 2645 + (2 terms) + O(x ) 2646 2647 2648COMMENT Test printing for negative expansion point; 2649 2650 2651taylor(sin x,x, -1, 6); 2652 2653 2654 sin(1) 2 cos(1) 3 2655 - sin(1) + cos(1)*(x + 1) + --------*(x + 1) - --------*(x + 1) 2656 2 6 2657 2658 sin(1) 4 7 2659 - --------*(x + 1) + (2 terms) + O((x + 1) ) 2660 24 2661 2662 2663taylor(sin x,x, -1/2, 6); 2664 2665 2666 1 1 2667 sin(---) cos(---) 2668 1 1 1 2 1 2 2 1 3 2669 - sin(---) + cos(---)*(x + ---) + ----------*(x + ---) - ----------*(x + ---) 2670 2 2 2 2 2 6 2 2671 2672 1 2673 sin(---) 2674 2 1 4 1 7 2675 - ----------*(x + ---) + (2 terms) + O((x + ---) ) 2676 24 2 2 2677 2678 2679%%% showtime; 2680 2681COMMENT That's all, folks; 2682 2683 2684end; 2685 2686