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(ei(x),x,0,4); 1502 1503 1504 1 2 1 3 1 4 5 1505log(x) - psi(1) + (x + ---*x + ----*x + ----*x + O(x )) 1506 4 18 96 1507 1508 1509comment In the following we demonstrate the possibiblity to compute the 1510 expansion of a function which is given by a simple first order 1511 differential equation: the function myexp(x) is exp(-x^2); 1512 1513 1514operator myexp,myerf; 1515 1516 1517let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 1518 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 1519 1520 1521taylor(myexp(x),x,0,5); 1522 1523 1524 2 1 4 6 15251 - x + ---*x + O(x ) 1526 2 1527 1528taylor(myerf(x),x,0,5); 1529 1530 1531 2 2 3 1 5 6 1532----------*x - ------------*x + ------------*x + O(x ) 1533 sqrt(pi) 3*sqrt(pi) 5*sqrt(pi) 1534 1535clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1, 1536 df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0}; 1537 1538 1539clear myexp,erf; 1540 1541 1542 1543%%% showtime; 1544 1545comment There are two special operators, implicit_taylor and 1546 inverse_taylor, to compute the Taylor expansion of implicit 1547 or inverse functions; 1548 1549 1550 1551implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5); 1552 1553 1554 1 2 1 4 6 15551 - ---*x - ---*x + O(x ) 1556 2 8 1557 1558 1559implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20); 1560 1561 1562 1 2 1 4 1 6 5 8 7 10 21 15631 - ---*x - ---*x - ----*x - -----*x - -----*x + (5 terms) + O(x ) 1564 2 8 16 128 256 1565 1566 1567implicit_taylor(x+y^3-y,x,y,0,0,8); 1568 1569 1570 3 5 7 9 1571x + x + 3*x + 12*x + O(x ) 1572 1573 1574implicit_taylor(x+y^3-y,x,y,0,1,5); 1575 1576 1577 1 3 2 1 3 105 4 3 5 6 15781 - ---*x - ---*x - ---*x - -----*x - ---*x + O(x ) 1579 2 8 2 128 2 1580 1581 1582implicit_taylor(x+y^3-y,x,y,0,-1,5); 1583 1584 1585 1 3 2 1 3 105 4 3 5 6 1586 - 1 - ---*x + ---*x - ---*x + -----*x - ---*x + O(x ) 1587 2 8 2 128 2 1588 1589 1590implicit_taylor(y*e^y-x,x,y,0,0,5); 1591 1592 1593 2 3 3 8 4 125 5 6 1594x - x + ---*x - ---*x + -----*x + O(x ) 1595 2 3 24 1596 1597 1598comment This is the function exp(-1/x^2), which has an essential 1599 singularity at the point 0; 1600 1601 1602implicit_taylor(x^2*log y+1,x,y,0,0,3); 1603 1604 1605***** Computation of Taylor series of implicit function failed 1606 Input expression non-zero at given point 1607 1608 1609inverse_taylor(exp(x)-1,x,y,0,8); 1610 1611 1612 1 2 1 3 1 4 1 5 1 6 9 1613y - ---*y + ---*y - ---*y + ---*y - ---*y + (2 terms) + O(y ) 1614 2 3 4 5 6 1615 1616 1617inverse_taylor(exp(x),x,y,0,5); 1618 1619 1620 1 2 1 3 1 4 1 5 6 1621y - 1 - ---*(y - 1) + ---*(y - 1) - ---*(y - 1) + ---*(y - 1) + O((y - 1) ) 1622 2 3 4 5 1623 1624 1625inverse_taylor(sqrt(x),x,y,0,5); 1626 1627 1628 2 6 1629y + O(y ) 1630 1631 1632inverse_taylor(log(1+x),x,y,0,5); 1633 1634 1635 1 2 1 3 1 4 1 5 6 1636y + ---*y + ---*y + ----*y + -----*y + O(y ) 1637 2 6 24 120 1638 1639 1640inverse_taylor((e^x-e^(-x))/2,x,y,0,5); 1641 1642 1643 1 3 3 5 6 1644y - ---*y + ----*y + O(y ) 1645 6 40 1646 1647 1648comment In the next two cases the inverse functions have a branch 1649 point, therefore the computation fails; 1650 1651 1652inverse_taylor((e^x+e^(-x))/2,x,y,0,5); 1653 1654 1655***** Computation of Taylor series of inverse function failed 1656 1657 1658inverse_taylor(exp(x^2-1),x,y,0,5); 1659 1660 1661***** Computation of Taylor series of inverse function failed 1662 1663 1664inverse_taylor(exp(sqrt(x))-1,x,y,0,5); 1665 1666 1667 2 3 11 4 5 5 6 1668y - y + ----*y - ---*y + O(y ) 1669 12 6 1670 1671 1672inverse_taylor(x*exp(x),x,y,0,5); 1673 1674 1675 2 3 3 8 4 125 5 6 1676y - y + ---*y - ---*y + -----*y + O(y ) 1677 2 3 24 1678 1679 1680 1681%%% showtime; 1682 1683 1684comment An application is the problem posed by Prof. Stanley: 1685 we prove that the finite difference expression below 1686 corresponds to the given derivative expression; 1687 1688 1689operator diff,a,f,gg; 1690 1691 % We use gg to avoid conflict with high energy 1692 % physics operator. 1693 1694let diff(~f,~arg) => df(f,arg); 1695 1696 1697 1698derivative_expression := 1699diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) + 1700diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ; 1701 1702 1703derivative_expression := 2*a(x,y)*df(f(x,y),x,y)*df(gg(x,y),x)*df(gg(x,y),y) 1704 1705 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x,y)*df(gg(x,y),y) 1706 1707 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y,2) 1708 1709 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,y)*df(gg(x,y),x) 1710 1711 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,2)*df(gg(x,y),y) 1712 1713 + df(a(x,y),x)*df(f(x,y),y)*df(gg(x,y),x)*df(gg(x,y),y) 1714 1715 + df(a(x,y),y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y) 1716 1717 1718finite_difference_expression := 1719+a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1720+a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1721+a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1722+a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1723-f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1724-f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1725-f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1726-a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2) 1727-gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1728-gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1729-gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1730-a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1731+f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1732+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1733+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1734+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2) 1735-gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1736+gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1737-gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1738+gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1739-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1740-a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1741-a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1742+gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1743+a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1744+a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2) 1745-gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1746+gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1747-a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1748-a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1749+gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2) 1750+a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2) 1751+f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1752-f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2) 1753+f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1754-f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2) 1755+a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1756+a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1757+a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1758+a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1759-f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1760-f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1761-f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1762-a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2) 1763-gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1764-gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1765-gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1766-a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1767+f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1768+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1769+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1770+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2) 1771-gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1772+gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1773-gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1774+gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1775-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1776-a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1777-a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1778+gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1779+a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1780+a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2) 1781-gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1782+gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1783-a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1784-a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1785+gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2) 1786+a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2) 1787+f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1788-f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2) 1789+f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1790-f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2) 1791+f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 1792+f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 1793+f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2) 1794+a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2) 1795-f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 1796-f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2) 1797-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1798-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1799-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1800-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2) 1801+f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 1802+f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2) 1803-f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2) 1804+a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1805+a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1806+a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1807+a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1808-f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1809-f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1810-f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1811-a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2) 1812-gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1813-gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1814-gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1815-a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1816+f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1817+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1818+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1819+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2) 1820-gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1821+gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1822-gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1823+gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1824-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1825-a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1826-a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1827+gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1828+a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1829+a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2) 1830-gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1831+gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1832-a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1833-a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1834+gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2) 1835+a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2) 1836+f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1837-f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2) 1838+f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1839-f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2) 1840+a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1841+a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1842+a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1843+a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1844-f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1845-f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1846-f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1847-a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2) 1848-gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1849-gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1850-gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1851-a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1852+f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1853+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1854+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1855+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2) 1856-gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1857+gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1858-gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1859+gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1860-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1861-a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1862-a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1863+gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1864+a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1865+a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2) 1866-gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1867+gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1868-a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1869-a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1870+gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2) 1871+a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2) 1872+f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1873-f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2) 1874+f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1875-f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2) 1876+f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 1877+f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 1878+f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2) 1879+a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2) 1880-f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 1881-f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2) 1882-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1883-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1884-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1885-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2) 1886+f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 1887+f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2) 1888-f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2) 1889+f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2) 1890+a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2) 1891-f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2) 1892+f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2) 1893+a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2) 1894-f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2) 1895-a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$ 1896 1897 1898 1899comment We define abbreviations for the partial derivatives; 1900 1901 1902operator ax,ay,fx,fy,gx,gy; 1903 1904 1905operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 1906 1907 1908operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 1909 1910 1911operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy, 1912 gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 1913 1914 1915operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 1916 1917 1918operator gxxxxyy,gxxxyyy,gxxyyyy; 1919 1920 1921 1922operator_diff_rules := { 1923 df(a(~x,~y),~x) => ax(x,y), 1924 df(a(~x,~y),~y) => ay(x,y), 1925 df(f(~x,~y),~x) => fx(x,y), 1926 df(f(~x,~y),~y) => fy(x,y), 1927 df(gg(~x,~y),~x) => gx(x,y), 1928 df(gg(~x,~y),~y) => gy(x,y), 1929 1930 df(ax(~x,~y),~x) => axx(x,y), 1931 df(ax(~x,~y),~y) => axy(x,y), 1932 df(ay(~x,~y),~x) => axy(x,y), 1933 df(ay(~x,~y),~y) => ayy(x,y), 1934 df(fx(~x,~y),~x) => fxx(x,y), 1935 df(fx(~x,~y),~y) => fxy(x,y), 1936 df(fy(~x,~y),~x) => fxy(x,y), 1937 df(fy(~x,~y),~y) => fyy(x,y), 1938 df(gx(~x,~y),~x) => gxx(x,y), 1939 df(gx(~x,~y),~y) => gxy(x,y), 1940 df(gy(~x,~y),~x) => gxy(x,y), 1941 df(gy(~x,~y),~y) => gyy(x,y), 1942 1943 df(axx(~x,~y),~x) => axxx(x,y), 1944 df(axy(~x,~y),~x) => axxy(x,y), 1945 df(ayy(~x,~y),~x) => axyy(x,y), 1946 df(ayy(~x,~y),~y) => ayyy(x,y), 1947 df(fxx(~x,~y),~x) => fxxx(x,y), 1948 df(fxy(~x,~y),~x) => fxxy(x,y), 1949 df(fxy(~x,~y),~y) => fxyy(x,y), 1950 df(fyy(~x,~y),~x) => fxyy(x,y), 1951 df(fyy(~x,~y),~y) => fyyy(x,y), 1952 df(gxx(~x,~y),~x) => gxxx(x,y), 1953 df(gxx(~x,~y),~y) => gxxy(x,y), 1954 df(gxy(~x,~y),~x) => gxxy(x,y), 1955 df(gxy(~x,~y),~y) => gxyy(x,y), 1956 df(gyy(~x,~y),~x) => gxyy(x,y), 1957 df(gyy(~x,~y),~y) => gyyy(x,y), 1958 1959 df(axyy(~x,~y),~x) => axxyy(x,y), 1960 df(axxy(~x,~y),~x) => axxxy(x,y), 1961 df(ayyy(~x,~y),~x) => axyyy(x,y), 1962 df(fxxy(~x,~y),~x) => fxxxy(x,y), 1963 df(fxyy(~x,~y),~x) => fxxyy(x,y), 1964 df(fyyy(~x,~y),~x) => fxyyy(x,y), 1965 df(gxxx(~x,~y),~x) => gxxxx(x,y), 1966 df(gxxy(~x,~y),~x) => gxxxy(x,y), 1967 df(gxyy(~x,~y),~x) => gxxyy(x,y), 1968 df(gyyy(~x,~y),~x) => gxyyy(x,y), 1969 df(gyyy(~x,~y),~y) => gyyyy(x,y), 1970 1971 df(axxyy(~x,~y),~x) => axxxyy(x,y), 1972 df(axyyy(~x,~y),~x) => axxyyy(x,y), 1973 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 1974 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 1975 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 1976 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 1977 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 1978 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 1979 1980 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 1981 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 1982 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y) 1983}; 1984 1985 1986operator_diff_rules := {df(a(~x,~y),~x) => ax(x,y), 1987 1988 df(a(~x,~y),~y) => ay(x,y), 1989 1990 df(f(~x,~y),~x) => fx(x,y), 1991 1992 df(f(~x,~y),~y) => fy(x,y), 1993 1994 df(gg(~x,~y),~x) => gx(x,y), 1995 1996 df(gg(~x,~y),~y) => gy(x,y), 1997 1998 df(ax(~x,~y),~x) => axx(x,y), 1999 2000 df(ax(~x,~y),~y) => axy(x,y), 2001 2002 df(ay(~x,~y),~x) => axy(x,y), 2003 2004 df(ay(~x,~y),~y) => ayy(x,y), 2005 2006 df(fx(~x,~y),~x) => fxx(x,y), 2007 2008 df(fx(~x,~y),~y) => fxy(x,y), 2009 2010 df(fy(~x,~y),~x) => fxy(x,y), 2011 2012 df(fy(~x,~y),~y) => fyy(x,y), 2013 2014 df(gx(~x,~y),~x) => gxx(x,y), 2015 2016 df(gx(~x,~y),~y) => gxy(x,y), 2017 2018 df(gy(~x,~y),~x) => gxy(x,y), 2019 2020 df(gy(~x,~y),~y) => gyy(x,y), 2021 2022 df(axx(~x,~y),~x) => axxx(x,y), 2023 2024 df(axy(~x,~y),~x) => axxy(x,y), 2025 2026 df(ayy(~x,~y),~x) => axyy(x,y), 2027 2028 df(ayy(~x,~y),~y) => ayyy(x,y), 2029 2030 df(fxx(~x,~y),~x) => fxxx(x,y), 2031 2032 df(fxy(~x,~y),~x) => fxxy(x,y), 2033 2034 df(fxy(~x,~y),~y) => fxyy(x,y), 2035 2036 df(fyy(~x,~y),~x) => fxyy(x,y), 2037 2038 df(fyy(~x,~y),~y) => fyyy(x,y), 2039 2040 df(gxx(~x,~y),~x) => gxxx(x,y), 2041 2042 df(gxx(~x,~y),~y) => gxxy(x,y), 2043 2044 df(gxy(~x,~y),~x) => gxxy(x,y), 2045 2046 df(gxy(~x,~y),~y) => gxyy(x,y), 2047 2048 df(gyy(~x,~y),~x) => gxyy(x,y), 2049 2050 df(gyy(~x,~y),~y) => gyyy(x,y), 2051 2052 df(axyy(~x,~y),~x) => axxyy(x,y), 2053 2054 df(axxy(~x,~y),~x) => axxxy(x,y), 2055 2056 df(ayyy(~x,~y),~x) => axyyy(x,y), 2057 2058 df(fxxy(~x,~y),~x) => fxxxy(x,y), 2059 2060 df(fxyy(~x,~y),~x) => fxxyy(x,y), 2061 2062 df(fyyy(~x,~y),~x) => fxyyy(x,y), 2063 2064 df(gxxx(~x,~y),~x) => gxxxx(x,y), 2065 2066 df(gxxy(~x,~y),~x) => gxxxy(x,y), 2067 2068 df(gxyy(~x,~y),~x) => gxxyy(x,y), 2069 2070 df(gyyy(~x,~y),~x) => gxyyy(x,y), 2071 2072 df(gyyy(~x,~y),~y) => gyyyy(x,y), 2073 2074 df(axxyy(~x,~y),~x) => axxxyy(x,y), 2075 2076 df(axyyy(~x,~y),~x) => axxyyy(x,y), 2077 2078 df(fxxyy(~x,~y),~x) => fxxxyy(x,y), 2079 2080 df(fxyyy(~x,~y),~x) => fxxyyy(x,y), 2081 2082 df(gxxxy(~x,~y),~x) => gxxxxy(x,y), 2083 2084 df(gxxyy(~x,~y),~x) => gxxxyy(x,y), 2085 2086 df(gxyyy(~x,~y),~x) => gxxyyy(x,y), 2087 2088 df(gyyyy(~x,~y),~x) => gxyyyy(x,y), 2089 2090 df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y), 2091 2092 df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y), 2093 2094 df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)} 2095 2096 2097let operator_diff_rules; 2098 2099 2100 2101texp := taylor (finite_difference_expression, dx, 0, 1, dy, 0, 1); 2102 2103 2104texp := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y) 2105 2106 + 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) 2107 2108 + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y) 2109 2110 2 2 2111 + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y) + O(dx ,dy ) 2112 2113 2114comment You may also try to expand further but this needs a lot 2115 of CPU time. Therefore the following line is commented out; 2116 2117 2118%texp := taylor (finite_difference_expression, dx, 0, 2, dy, 0, 2); 2119 2120factor dx,dy; 2121 2122 2123 2124result := taylortostandard texp; 2125 2126 2127result := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y) 2128 2129 + 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) 2130 2131 + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y) 2132 2133 + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y) 2134 2135 2136derivative_expression - result; 2137 2138 21390 2140 2141 2142 2143clear diff(~f,~arg); 2144 2145 2146clearrules operator_diff_rules; 2147 2148 2149clear diff,a,f,gg; 2150 2151 2152clear ax,ay,fx,fy,gx,gy; 2153 2154 2155clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy; 2156 2157 2158clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy; 2159 2160 2161clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy; 2162 2163 2164clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy; 2165 2166 2167clear gxxxxyy,gxxxyyy,gxxyyyy; 2168 2169 2170 2171taylorprintterms := 5; 2172 2173 2174taylorprintterms := 5 2175 2176 2177off taylorautoexpand,taylorkeeporiginal; 2178 2179 2180 2181%%% showtime; 2182 2183comment An example provided by Alan Barnes: elliptic functions; 2184 2185 2186% Jacobi's elliptic functions 2187% sn(x,k), cn(x,k), dn(x,k). 2188% The modulus and complementary modulus are denoted by K and K!' 2189% usually written mathematically as k and k' respectively 2190% 2191% epsilon(x,k) is the incomplete elliptic integral of the second kind 2192% usually written mathematically as E(x,k) 2193% 2194% KK(k) is the complete elliptic integral of the first kind 2195% usually written mathematically as K(k) 2196% K(k) = arcsn(1,k) 2197% KK!'(k) is the complementary complete integral 2198% usually written mathematically as K'(k) 2199% NB. K'(k) = K(k') 2200% EE(k) is the complete elliptic integral of the second kind 2201% usually written mathematically as E(k) 2202% EE!'(k) is the complementary complete integral 2203% usually written mathematically as E'(k) 2204% NB. E'(k) = E(k') 2205 2206operator sn, cn, dn, epsilon; 2207 2208 2209 2210elliptic_rules := { 2211% Differentiation rules for basic functions 2212 df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 2213 df(cn(~x,~k),~x) => -sn(x,k)*dn(x,k), 2214 df(dn(~x,~k),~x) => -k^2*sn(x,k)*cn(x,k), 2215 df(epsilon(~x,~k),~x)=> dn(x,k)^2, 2216 2217% k-derivatives 2218% DF Lawden Elliptic Functions & Applications Springer (1989) 2219 df(sn(~x,~k),~k) => (k*sn(x,k)*cn(x,k)^2 2220 -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2) 2221 + x*cn(x,k)*dn(x,k)/k, 2222 df(cn(~x,~k),~k) => (-k*sn(x,k)^2*cn(x,k) 2223 +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2) 2224 - x*sn(x,k)*dn(x,k)/k, 2225 df(dn(~x,~k),~k) => k*(-sn(x,k)^2*dn(x,k) 2226 +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2) 2227 - k*x*sn(x,k)*cn(x,k), 2228 df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k) 2229 -epsilon(x,k)*cn(x,k)^2)/(1-k^2) 2230 -k*x*sn(x,k)^2, 2231 2232% parity properties 2233 sn(-~x,~k) => -sn(x,k), 2234 cn(-~x,~k) => cn(x,k), 2235 dn(-~x,~k) => dn(x,k), 2236 epsilon(-~x,~k) => -epsilon(x,k), 2237 sn(~x,-~k) => sn(x,k), 2238 cn(~x,-~k) => cn(x,k), 2239 dn(~x,-~k) => dn(x,k), 2240 epsilon(~x,-~k) => epsilon(x,k), 2241 2242 2243% behaviour at zero 2244 sn(0,~k) => 0, 2245 cn(0,~k) => 1, 2246 dn(0,~k) => 1, 2247 epsilon(0,~k) => 0, 2248 2249 2250% degenerate cases of modulus 2251 sn(~x,0) => sin(x), 2252 cn(~x,0) => cos(x), 2253 dn(~x,0) => 1, 2254 epsilon(~x,0) => x, 2255 2256 sn(~x,1) => tanh(x), 2257 cn(~x,1) => 1/cosh(x), 2258 dn(~x,1) => 1/cosh(x), 2259 epsilon(~x,1) => tanh(x) 2260}; 2261 2262 2263elliptic_rules := {df(sn(~x,~k),~x) => cn(x,k)*dn(x,k), 2264 2265 df(cn(~x,~k),~x) => - sn(x,k)*dn(x,k), 2266 2267 2 2268 df(dn(~x,~k),~x) => - k *sn(x,k)*cn(x,k), 2269 2270 2 2271 df(epsilon(~x,~k),~x) => dn(x,k) , 2272 2273 2 dn(x,k) 2274 k*sn(x,k)*cn(x,k) - epsilon(x,k)*cn(x,k)*--------- 2275 k 2276 df(sn(~x,~k),~k) => ----------------------------------------------------- 2277 2 2278 1 - k 2279 2280 dn(x,k) 2281 + x*cn(x,k)*---------, 2282 k 2283 2284 2 dn(x,k) 2285 - k*sn(x,k) *cn(x,k) + epsilon(x,k)*sn(x,k)*--------- 2286 k 2287 df(cn(~x,~k),~k) => -------------------------------------------------------- 2288 2 2289 1 - k 2290 2291 dn(x,k) 2292 - x*sn(x,k)*---------, 2293 k 2294 2295 2 2296 - sn(x,k) *dn(x,k) + epsilon(x,k)*sn(x,k)*cn(x,k) 2297 df(dn(~x,~k),~k) => k*---------------------------------------------------- 2298 2 2299 1 - k 2300 2301 - k*x*sn(x,k)*cn(x,k), 2302 2303 df(epsilon(~x,~k),~k) 2304 2305 2 2306 sn(x,k)*cn(x,k)*dn(x,k) - epsilon(x,k)*cn(x,k) 2 2307 => k*------------------------------------------------- - k*x*sn(x,k) , 2308 2 2309 1 - k 2310 2311 sn( - ~x,~k) => - sn(x,k), 2312 2313 cn( - ~x,~k) => cn(x,k), 2314 2315 dn( - ~x,~k) => dn(x,k), 2316 2317 epsilon( - ~x,~k) => - epsilon(x,k), 2318 2319 sn(~x, - ~k) => sn(x,k), 2320 2321 cn(~x, - ~k) => cn(x,k), 2322 2323 dn(~x, - ~k) => dn(x,k), 2324 2325 epsilon(~x, - ~k) => epsilon(x,k), 2326 2327 sn(0,~k) => 0, 2328 2329 cn(0,~k) => 1, 2330 2331 dn(0,~k) => 1, 2332 2333 epsilon(0,~k) => 0, 2334 2335 sn(~x,0) => sin(x), 2336 2337 cn(~x,0) => cos(x), 2338 2339 dn(~x,0) => 1, 2340 2341 epsilon(~x,0) => x, 2342 2343 sn(~x,1) => tanh(x), 2344 2345 1 2346 cn(~x,1) => ---------, 2347 cosh(x) 2348 2349 1 2350 dn(~x,1) => ---------, 2351 cosh(x) 2352 2353 epsilon(~x,1) => tanh(x)} 2354 2355 2356let elliptic_rules; 2357 2358 2359 2360hugo := taylor(sn(x,k),k,0,6); 2361 2362 2363 2 2 2364 cos(x)*(cos(x) *x + cos(x)*sin(x) + sin(x) *x - 2*x) 2 2365hugo := sin(x) + ------------------------------------------------------*k + ( 2366 4 2367 2368 5 4 2 4 2369 cos(x) *x - 2*cos(x) *sin(x)*x + 5*cos(x) *sin(x) 2370 2371 3 2 3 2 3 2 2372 - 10*cos(x) *sin(x) *x + 6*cos(x) *x - 4*cos(x) *sin(x) *x 2373 2374 2 3 2 2 2 2375 + cos(x) *sin(x) + 8*cos(x) *sin(x)*x + 4*cos(x) *sin(x) 2376 2377 4 2 2378 - 11*cos(x)*sin(x) *x + 30*cos(x)*sin(x) *x - 16*cos(x)*x 2379 2380 5 2 3 2 2 4 2381 - 2*sin(x) *x + 8*sin(x) *x - 8*sin(x)*x )/64*k + ( 2382 2383 7 3 7 6 2 2384 - 6*cos(x) *x + 17*cos(x) *x - 99*cos(x) *sin(x)*x 2385 2386 6 5 2 3 5 2 2387 + 21*cos(x) *sin(x) - 18*cos(x) *sin(x) *x - 71*cos(x) *sin(x) *x 2388 2389 5 3 5 4 3 2 2390 + 36*cos(x) *x - 18*cos(x) *x - 135*cos(x) *sin(x) *x 2391 2392 4 3 4 2 4 2393 - 133*cos(x) *sin(x) + 324*cos(x) *sin(x)*x + 172*cos(x) *sin(x) 2394 2395 3 4 3 3 4 2396 - 18*cos(x) *sin(x) *x - 13*cos(x) *sin(x) *x 2397 2398 3 2 3 3 2 3 3 2399 + 72*cos(x) *sin(x) *x - 156*cos(x) *sin(x) *x - 72*cos(x) *x 2400 2401 3 2 5 2 2 5 2402 + 160*cos(x) *x + 27*cos(x) *sin(x) *x - 118*cos(x) *sin(x) 2403 2404 2 3 2 2 2 2405 + 176*cos(x) *sin(x) - 108*cos(x) *sin(x)*x + 32*cos(x) *sin(x) 2406 2407 6 3 6 4 3 2408 - 6*cos(x)*sin(x) *x + 75*cos(x)*sin(x) *x + 36*cos(x)*sin(x) *x 2409 2410 4 2 3 2 2411 - 498*cos(x)*sin(x) *x - 72*cos(x)*sin(x) *x + 888*cos(x)*sin(x) *x 2412 2413 3 7 2 5 2 2414 + 48*cos(x)*x - 384*cos(x)*x + 63*sin(x) *x - 324*sin(x) *x 2415 2416 3 2 2 6 7 2417 + 540*sin(x) *x - 288*sin(x)*x )/2304*k + O(k ) 2418 2419otto := taylor(cn(x,k),k,0,6); 2420 2421 2422 2 2 2423 sin(x)*( - cos(x) *x - cos(x)*sin(x) - sin(x) *x + 2*x) 2 2424otto := cos(x) + ---------------------------------------------------------*k + 2425 4 2426 2427 5 2 4 3 2 2 2428 ( - 2*cos(x) *x - 5*cos(x) *sin(x)*x - 4*cos(x) *sin(x) *x 2429 2430 3 2 3 2 2 3 2431 - 7*cos(x) *sin(x) + 8*cos(x) *x + 2*cos(x) *sin(x) *x 2432 2433 2 4 2 4 2434 + 2*cos(x) *sin(x)*x - 2*cos(x)*sin(x) *x - 3*cos(x)*sin(x) 2435 2436 2 2 2 2 5 2437 + 8*cos(x)*sin(x) *x - 4*cos(x)*sin(x) - 8*cos(x)*x + 7*sin(x) *x 2438 2439 3 4 7 2 2440 - 22*sin(x) *x + 16*sin(x)*x)/64*k + ( - 9*cos(x) *x 2441 2442 6 3 6 5 2 2 2443 + 6*cos(x) *sin(x)*x - 71*cos(x) *sin(x)*x + 135*cos(x) *sin(x) *x 2444 2445 5 2 5 2 4 3 3 2446 - 66*cos(x) *sin(x) - 36*cos(x) *x + 18*cos(x) *sin(x) *x 2447 2448 4 3 4 3 4 2449 - cos(x) *sin(x) *x - 36*cos(x) *sin(x)*x + 18*cos(x) *sin(x)*x 2450 2451 3 4 2 3 4 2452 + 297*cos(x) *sin(x) *x + 61*cos(x) *sin(x) 2453 2454 3 2 2 3 2 3 2 2455 - 720*cos(x) *sin(x) *x - 208*cos(x) *sin(x) + 252*cos(x) *x 2456 2457 2 5 3 2 5 2458 + 18*cos(x) *sin(x) *x + 31*cos(x) *sin(x) *x 2459 2460 2 3 3 2 3 2461 - 72*cos(x) *sin(x) *x - 24*cos(x) *sin(x) *x 2462 2463 2 3 2 6 2 2464 + 72*cos(x) *sin(x)*x + 56*cos(x) *sin(x)*x + 153*cos(x)*sin(x) *x 2465 2466 6 4 2 4 2467 + 91*cos(x)*sin(x) - 684*cos(x)*sin(x) *x - 212*cos(x)*sin(x) 2468 2469 2 2 2 2 2470 + 900*cos(x)*sin(x) *x - 32*cos(x)*sin(x) - 288*cos(x)*x 2471 2472 7 3 7 5 3 5 2473 + 6*sin(x) *x - 39*sin(x) *x - 36*sin(x) *x + 318*sin(x) *x 2474 2475 3 3 3 3 2476 + 72*sin(x) *x - 672*sin(x) *x - 48*sin(x)*x + 384*sin(x)*x)/2304 2477 2478 6 7 2479 *k + O(k ) 2480 2481taylorcombine(hugo^2 + otto^2); 2482 2483 2484 2 2 7 2485cos(x) + sin(x) + O(k ) 2486 2487 2488clearrules elliptic_rules; 2489 2490 2491 2492clear sn, cn, dn, epsilon; 2493 2494 2495 2496%%% showtime; 2497 2498comment That's all, folks; 2499 2500 2501end; 2502 2503Tested on x86_64-pc-windows CSL 2504Time (counter 1): 126 ms 2505 2506End of Lisp run after 0.12+0.04 seconds 2507real 0.33 2508user 0.03 2509sys 0.06 2510