1on errcont; 2 3 4bounds (x,x=(1 .. 2)); 5 6 71..2 8 9bounds (2*x,x=(1 .. 2)); 10 11 122..4 13 14bounds (x**3,x=(1 .. 2)); 15 16 171..8 18 19bounds (x*y,x=(1 .. 2),y=(-1 .. 0)); 20 21 22 - 2..0 23 24bounds (x**3+y,x=(1 .. 2),y=(-1 .. 0)); 25 26 270..8 28 29bounds (x**3/y,{x=(1 .. 2),y=(-1 .. -0.5)}); 30 31 32 - 16..-1 33 34bounds (x**3/y,x=(1 .. 2),y=(-1 .. -0.5)); 35 36 37 - 16..-1 38 39 % unbounded expression (pole at y=0) 40bounds (x**3/y,x=(1 .. 2),y=(-1 .. 0.5)); 41 42 43***** unbounded in range 44 45 46on rounded; 47 48 49bounds(e**x,x=(1 .. 2)); 50 51 522.71828182846..7.38905609893 53 54bounds((1/2)**x,x=(1 .. 2)); 55 56 570.25..0.5 58 59off rounded; 60 61 62 63bounds(abs x,x=(1 .. 2)); 64 65 661..2 67 68bounds(abs x,x=(-3 .. 2)); 69 70 710..3 72 73bounds(abs x,x=(-3 .. -2)); 74 75 762..3 77 78 79bounds(sin x,x=(1 .. 2)); 80 81 82 - 1..1 83 84 85on rounded; 86 87 88 89bounds(sin x,x=(1 .. 2)); 90 91 920.841470984808..1 93 94bounds(sin x,x=(1 .. 10)); 95 96 97 - 1..1 98 99bounds(sin x,x=(1001 .. 1002)); 100 101 1020.167266541974..0.919990597586 103 104 105bounds(log x,x=(1 .. 10)); 106 107 1080..2.30258509299 109 110 111bounds(tan x,x=(1 .. 1.1)); 112 113 1141.55740772465..1.96475965725 115 116 117bounds(cot x,x=(1 .. 1.1)); 118 119 1200.508968105239..0.642092615934 121 122bounds(asin x,x=(-0.6 .. 0.6)); 123 124 125 - 0.643501108793..0.643501108793 126 127bounds(acos x,x=(-0.6 .. 0.6)); 128 129 1300.927295218002..2.21429743559 131 132 133bounds(sqrt(x),x=(1 .. 1.1)); 134 135 1361..1.04880884817 137 138bounds(x**(7/3),x=(1 .. 1.1)); 139 140 1411..1.2490589397 142 143bounds(x**y,x=(1 .. 1.1),y=(2 .. 4)); 144 145 1461..1.4641 147 148 149off rounded; 150 151 152 153 154% MINIMA (steepest descent) 155 156% Rosenbrock function (minimum extremely hard to find). 157fktn := 100*(x1^2-x2)^2 + (1-x1)^2; 158 159 160 4 2 2 2 161fktn := 100*x1 - 200*x1 *x2 + x1 - 2*x1 + 100*x2 + 1 162 163num_min(fktn, x1=-1.2, x2=1, accuracy=6); 164 165 166{0.00000021870243927,{x1=0.999532844813,x2=0.999068072135}} 167 168 169% infinitely many local minima 170num_min(sin(x)+x/5, x=1); 171 172 173{ - 1.33422674663,{x= - 1.77215430279}} 174 175 176% bivariate polynomial 177num_min(x^4 + 3 x^2 * y + 5 y^2 + x + y, x=0.1, y=0.2); 178 179 180{ - 0.832523282274,{x= - 0.889601609042,y= - 0.33989805551}} 181 182 183 184% ROOTS (non polynomial: damped Newton) 185 186 num_solve (cos x -x, x=0,accuracy=6); 187 188 189{x=0.739085133215} 190 191 192 % automatically randomized starting point 193num_solve (cos x -x,x, accuracy=6); 194 195 196{x=0.739085133215} 197 198 199 % syntactical errors: forms do not evaluate to purely 200 % numerical values 201num_solve (cos x -x, x=a); 202 203 204***** a invalid as number 205 206num_solve (cos x -a, x=0); 207 208 209***** error during function evaluation (e.g. singularity) 210 211 212num_solve (sin x = 0, x=3); 213 214 215{x=3.14159265359} 216 217 218 % blows up: no real solution exists 219num_solve(sin x = 2, x=1); 220 221 222***** Newton method does not converge 223 224 225 % solution in complex plane(only fond with complex starting point): 226on complex; 227 228 229*** Domain mode rounded changed to complex-rounded 230 231num_solve(sin x = 2, x=1+i); 232 233 234{x=1.57079632542 + 1.31695789681*i} 235 236off complex; 237 238 239*** Domain mode complex-rounded changed to rounded 240 241 242 % blows up for derivative 0 in starting point 243num_solve(x^2-1, x=0); 244 245 246***** division by zero 247 248 249 % succeeds because of perturbed starting point 250num_solve(x^2-1, x=0.1); 251 252 253{x=1.00000000033} 254 255 256 % bivariate equation system 257num_solve({sin x=cos y, x + y = 1},{x=1,y=2}); 258 259 260{x= - 4.99778714379,y=5.99778714379} 261 262on rounded,evallhseqp; 263 264 265sub(ws,{sin x=cos y, x + y = 1}); 266 267 268{0.959549629982=0.959549629988,1=1} 269 270off rounded,evallhseqp; 271 272 273 274 % temporal member of the Barry Simon test sequence 275sys :={sin (x) + y^2 + log(z) = 7, 276 3*x + 2^y - z^3 = -90, 277 x^2 + y^2 + z^(1/2) = 16}; 278 279 280 2 281sys := {sin(x) + y + log(z)=7, 282 283 y 3 284 3*x + 2 - z =-90, 285 286 2 2 1/2 287 x + y + z =16} 288 289sol:=num_solve(sys,{x=1,y=1,z=1}); 290 291 292sol := {x=2.93087675819,y= - 2.29328251176,z=4.62601269017} 293 294on rounded; 295 296 297for each s in sys collect sub(sol,lhs s-rhs s); 298 299 300{0,0,0} 301 302off rounded; 303 304 305clear sys,sol; 306 307 308 309 % 2 examples taken from Nowak/Weimann (Tech.Rep TR91-10, ZIB Berlin) 310 311 % #1: exp/sin combination 312 313on rounded; 314 315 316sys := {e**(x1**2 + x2**2)-3, x1 + x2 - sin(3(x1 + x2))}; 317 318 319 2 2 320 x1 + x2 321sys := {e - 3, 322 323 - sin(3*x1 + 3*x2) + x1 + x2} 324 325num_solve(sys,x1=0.81, x2=0.82); 326 327 328*** precision increased to 14 329 330*** precision increased to 18 331 332*** precision increased to 20 333 334{x1= - 0.256625076922,x2=1.01624596361} 335 336sub(ws,sys); 337 338 339{0,0} 340 341 342 % 2nd example (semiconductor simulation), here computed with 343 % intermediate steps printed 344 345alpha := 38.683; 346 347 348alpha := 38.683 349 350ni := 1.22e10; 351 352 353ni := 1.22e+10 354 355v := 100; 356 357 358v := 100 359 360d := 1e17; 361 362 363d := 1.0e+17 364 365sys := { e**(alpha*(x3-x1)) - e**(alpha*(x1-x2)) - d/ni, 366 x2, 367 x3, 368 e**(alpha*(x6-x4)) - e**(alpha*(x4-x5)) + d/ni, 369 x5 - v, 370 x6 - v}; 371 372 373 77.366*x1 38.683*x1 + 38.683*x2 374sys := {( - e - 8.19672131148e+6*e 375 376 38.683*x2 + 38.683*x3 38.683*x1 + 38.683*x2 377 + e )/e , 378 379 x2, 380 381 x3, 382 383 77.366*x4 38.683*x4 + 38.683*x5 384 ( - e + 8.19672131148e+6*e 385 386 38.683*x5 + 38.683*x6 38.683*x4 + 38.683*x5 387 + e )/e , 388 389 x5 - 100, 390 391 x6 - 100} 392 393on trnumeric; 394 395 396num_solve(sys,x1=1,x2=2,x3=3,x4=4,x5=5,x6=6,iterations=100); 397 398 399*** computing symbolic Jacobian 400 401*** starting Newton iteration 402 403*** Newton iteration 404 405*** evalute function in the initial point 406 407*** evaluate Jacobian (or its inverse) 408 409*** compute the next point 410 411*** evalute function in the next point 412 413*** compute the point difference 414 415*** compute the size of the point difference 416 4171. residue=(1.46329673989e+33 , 0 , 0 , 1.46329673988e+33 , 0.0 , 0.0) 418 419, step length=95.0 420 at ( - 1.97414885092 , 0 , 0 , 98.0258511491 , 100.0 , 100.0) 421 422*** evaluate Jacobian (or its inverse) 423 424*** compute the next point 425 426*** evalute function in the next point 427 428*** compute the point difference 429 430*** compute the size of the point difference 431 4322. residue=(5.38316786938e+32 , 0.0 , 0.0 , 5.38316786935e+32 , 0.0 , 0.0) 433 434, step length=0.0258511490836 435 at ( - 1.94829770183 , 0.0 , 0.0 , 98.0517022982 , 100.0 , 100.0) 436 437*** evaluate Jacobian (or its inverse) 438 439*** compute the next point 440 441*** evalute function in the next point 442 443*** compute the point difference 444 445*** compute the size of the point difference 446 4473. residue=(1.98035678752e+32 , 0.0 , 0.0 , 1.98035678751e+32 , 0.0 , 0.0) 448 449, step length=0.0258511490836 450 at ( - 1.92244655275 , 0.0 , 0.0 , 98.0775534473 , 100.0 , 100.0) 451 452*** evaluate Jacobian (or its inverse) 453 454*** compute the next point 455 456*** evalute function in the next point 457 458*** compute the point difference 459 460*** compute the size of the point difference 461 4624. residue=(7.28532548313e+31 , 0.0 , 0.0 , 7.28532548309e+31 , 0.0 , 0.0) 463 464, step length=0.0258511490836 465 at ( - 1.89659540367 , 0.0 , 0.0 , 98.1034045963 , 100.0 , 100.0) 466 467*** evaluate Jacobian (or its inverse) 468 469*** compute the next point 470 471*** evalute function in the next point 472 473*** compute the point difference 474 475*** compute the size of the point difference 476 4775. residue=(2.68012146749e+31 , 0.0 , 0.0 , 2.68012146747e+31 , 0.0 , 0.0) 478 479, step length=0.0258511490836 480 at ( - 1.87074425458 , 0.0 , 0.0 , 98.1292557454 , 100.0 , 100.0) 481 482*** evaluate Jacobian (or its inverse) 483 484*** compute the next point 485 486*** evalute function in the next point 487 488*** compute the point difference 489 490*** compute the size of the point difference 491 4926. residue=(9.8596158773e+30 , 0.0 , 0.0 , 9.85961587725e+30 , 0.0 , 0.0) 493 494, step length=0.0258511490836 495 at ( - 1.8448931055 , 0.0 , 0.0 , 98.1551068945 , 100.0 , 100.0) 496 497*** evaluate Jacobian (or its inverse) 498 499*** compute the next point 500 501*** evalute function in the next point 502 503*** compute the point difference 504 505*** compute the size of the point difference 506 5077. residue=(3.62714997911e+30 , 0.0 , 0.0 , 3.62714997909e+30 , 0.0 , 0.0) 508 509, step length=0.0258511490836 510 at ( - 1.81904195641 , 0.0 , 0.0 , 98.1809580436 , 100.0 , 100.0) 511 512*** evaluate Jacobian (or its inverse) 513 514*** compute the next point 515 516*** evalute function in the next point 517 518*** compute the point difference 519 520*** compute the size of the point difference 521 5228. residue=(1.33435390736e+30 , 0.0 , 0.0 , 1.33435390735e+30 , 0.0 , 0.0) 523 524, step length=0.0258511490836 525 at ( - 1.79319080733 , 0.0 , 0.0 , 98.2068091927 , 100.0 , 100.0) 526 527*** evaluate Jacobian (or its inverse) 528 529*** compute the next point 530 531*** evalute function in the next point 532 533*** compute the point difference 534 535*** compute the size of the point difference 536 5379. residue=(4.90881369764e+29 , 0.0 , 0.0 , 4.90881369762e+29 , 0.0 , 0.0) 538 539, step length=0.0258511490836 540 at ( - 1.76733965825 , 0.0 , 0.0 , 98.2326603418 , 100.0 , 100.0) 541 542*** evaluate Jacobian (or its inverse) 543 544*** compute the next point 545 546*** evalute function in the next point 547 548*** compute the point difference 549 550*** compute the size of the point difference 551 55210. residue=(1.8058516399e+29 , 0.0 , 0.0 , 1.80585163991e+29 , 0.0 , 0.0) 553 554, step length=0.0258511490836 555 at ( - 1.74148850916 , 0.0 , 0.0 , 98.2585114908 , 100.0 , 100.0) 556 557*** evaluate Jacobian (or its inverse) 558 559*** compute the next point 560 561*** evalute function in the next point 562 563*** compute the point difference 564 565*** compute the size of the point difference 566 56711. residue=(6.64335692126e+28 , 0.0 , 0.0 , 6.64335692127e+28 , 0.0 , 0.0) 568 569, step length=0.0258511490836 570 at ( - 1.71563736008 , 0.0 , 0.0 , 98.2843626399 , 100.0 , 100.0) 571 572*** evaluate Jacobian (or its inverse) 573 574*** compute the next point 575 576*** evalute function in the next point 577 578*** compute the point difference 579 580*** compute the size of the point difference 581 58212. residue=(2.4439544317e+28 , 0.0 , 0.0 , 2.4439544317e+28 , 0.0 , 0.0) 583 584, step length=0.0258511490836 585 at ( - 1.689786211 , 0.0 , 0.0 , 98.310213789 , 100.0 , 100.0) 586 587*** evaluate Jacobian (or its inverse) 588 589*** compute the next point 590 591*** evalute function in the next point 592 593*** compute the point difference 594 595*** compute the size of the point difference 596 59713. residue=(8.99080590581e+27 , 0.0 , 0.0 , 8.99080590582e+27 , 0.0 , 0.0) 598 599, step length=0.0258511490836 600 at ( - 1.66393506191 , 0.0 , 0.0 , 98.3360649381 , 100.0 , 100.0) 601 602*** evaluate Jacobian (or its inverse) 603 604*** compute the next point 605 606*** evalute function in the next point 607 608*** compute the point difference 609 610*** compute the size of the point difference 611 61214. residue=(3.30753265231e+27 , 0.0 , 0.0 , 3.30753265232e+27 , 0.0 , 0.0) 613 614, step length=0.0258511490836 615 at ( - 1.63808391283 , 0.0 , 0.0 , 98.3619160872 , 100.0 , 100.0) 616 617*** evaluate Jacobian (or its inverse) 618 619*** compute the next point 620 621*** evalute function in the next point 622 623*** compute the point difference 624 625*** compute the size of the point difference 626 62715. residue=(1.21677326379e+27 , 0.0 , 0.0 , 1.21677326379e+27 , 0.0 , 0.0) 628 629, step length=0.0258511490836 630 at ( - 1.61223276375 , 0.0 , 0.0 , 98.3877672363 , 100.0 , 100.0) 631 632*** evaluate Jacobian (or its inverse) 633 634*** compute the next point 635 636*** evalute function in the next point 637 638*** compute the point difference 639 640*** compute the size of the point difference 641 64216. residue=(4.47625868315e+26 , 0.0 , 0.0 , 4.47625868316e+26 , 0.0 , 0.0) 643 644, step length=0.0258511490836 645 at ( - 1.58638161466 , 0.0 , 0.0 , 98.4136183853 , 100.0 , 100.0) 646 647*** evaluate Jacobian (or its inverse) 648 649*** compute the next point 650 651*** evalute function in the next point 652 653*** compute the point difference 654 655*** compute the size of the point difference 656 65717. residue=(1.64672354289e+26 , 0.0 , 0.0 , 1.6467235429e+26 , 0.0 , 0.0) 658 659, step length=0.0258511490836 660 at ( - 1.56053046558 , 0.0 , 0.0 , 98.4394695344 , 100.0 , 100.0) 661 662*** evaluate Jacobian (or its inverse) 663 664*** compute the next point 665 666*** evalute function in the next point 667 668*** compute the point difference 669 670*** compute the size of the point difference 671 67218. residue=(6.05795736724e+25 , 0.0 , 0.0 , 6.05795736725e+25 , 0.0 , 0.0) 673 674, step length=0.0258511490836 675 at ( - 1.5346793165 , 0.0 , 0.0 , 98.4653206835 , 100.0 , 100.0) 676 677*** evaluate Jacobian (or its inverse) 678 679*** compute the next point 680 681*** evalute function in the next point 682 683*** compute the point difference 684 685*** compute the size of the point difference 686 68719. residue=(2.2285979709e+25 , 0.0 , 0.0 , 2.2285979709e+25 , 0.0 , 0.0) 688 689, step length=0.0258511490836 690 at ( - 1.50882816741 , 0.0 , 0.0 , 98.4911718326 , 100.0 , 100.0) 691 692*** evaluate Jacobian (or its inverse) 693 694*** compute the next point 695 696*** evalute function in the next point 697 698*** compute the point difference 699 700*** compute the size of the point difference 701 70220. residue=(8.19855376131e+24 , 0.0 , 0.0 , 8.19855376132e+24 , 0.0 , 0.0) 703 704, step length=0.0258511490836 705 at ( - 1.48297701833 , 0.0 , 0.0 , 98.5170229817 , 100.0 , 100.0) 706 707*** evaluate Jacobian (or its inverse) 708 709*** compute the next point 710 711*** evalute function in the next point 712 713*** compute the point difference 714 715*** compute the size of the point difference 716 71721. residue=(3.01607937612e+24 , 0.0 , 0.0 , 3.01607937613e+24 , 0.0 , 0.0) 718 719, step length=0.0258511490836 720 at ( - 1.45712586924 , 0.0 , 0.0 , 98.5428741308 , 100.0 , 100.0) 721 722*** evaluate Jacobian (or its inverse) 723 724*** compute the next point 725 726*** evalute function in the next point 727 728*** compute the point difference 729 730*** compute the size of the point difference 731 73222. residue=(1.10955359542e+24 , 0.0 , 0.0 , 1.10955359542e+24 , 0.0 , 0.0) 733 734, step length=0.0258511490836 735 at ( - 1.43127472016 , 0.0 , 0.0 , 98.5687252798 , 100.0 , 100.0) 736 737*** evaluate Jacobian (or its inverse) 738 739*** compute the next point 740 741*** evalute function in the next point 742 743*** compute the point difference 744 745*** compute the size of the point difference 746 74723. residue=(4.08181956632e+23 , 0.0 , 0.0 , 4.08181956633e+23 , 0.0 , 0.0) 748 749, step length=0.0258511490836 750 at ( - 1.40542357108 , 0.0 , 0.0 , 98.5945764289 , 100.0 , 100.0) 751 752*** evaluate Jacobian (or its inverse) 753 754*** compute the next point 755 756*** evalute function in the next point 757 758*** compute the point difference 759 760*** compute the size of the point difference 761 76224. residue=(1.50161750102e+23 , 0.0 , 0.0 , 1.50161750102e+23 , 0.0 , 0.0) 763 764, step length=0.0258511490836 765 at ( - 1.37957242199 , 0.0 , 0.0 , 98.620427578 , 100.0 , 100.0) 766 767*** evaluate Jacobian (or its inverse) 768 769*** compute the next point 770 771*** evalute function in the next point 772 773*** compute the point difference 774 775*** compute the size of the point difference 776 77725. residue=(5.52414207128e+22 , 0.0 , 0.0 , 5.52414207133e+22 , 0.0 , 0.0) 778 779, step length=0.0258511490836 780 at ( - 1.35372127291 , 0.0 , 0.0 , 98.6462787271 , 100.0 , 100.0) 781 782*** evaluate Jacobian (or its inverse) 783 784*** compute the next point 785 786*** evalute function in the next point 787 788*** compute the point difference 789 790*** compute the size of the point difference 791 79226. residue=(2.03221829814e+22 , 0.0 , 0.0 , 2.03221829815e+22 , 0.0 , 0.0) 793 794, step length=0.0258511490836 795 at ( - 1.32787012383 , 0.0 , 0.0 , 98.6721298762 , 100.0 , 100.0) 796 797*** evaluate Jacobian (or its inverse) 798 799*** compute the next point 800 801*** evalute function in the next point 802 803*** compute the point difference 804 805*** compute the size of the point difference 806 80727. residue=(7.47611331856e+21 , 0.0 , 0.0 , 7.47611331863e+21 , 0.0 , 0.0) 808 809, step length=0.0258511490836 810 at ( - 1.30201897474 , 0.0 , 0.0 , 98.6979810253 , 100.0 , 100.0) 811 812*** evaluate Jacobian (or its inverse) 813 814*** compute the next point 815 816*** evalute function in the next point 817 818*** compute the point difference 819 820*** compute the size of the point difference 821 82228. residue=(2.75030838977e+21 , 0.0 , 0.0 , 2.75030838979e+21 , 0.0 , 0.0) 823 824, step length=0.0258511490836 825 at ( - 1.27616782566 , 0.0 , 0.0 , 98.7238321743 , 100.0 , 100.0) 826 827*** evaluate Jacobian (or its inverse) 828 829*** compute the next point 830 831*** evalute function in the next point 832 833*** compute the point difference 834 835*** compute the size of the point difference 836 83729. residue=(1.01178191348e+21 , 0.0 , 0.0 , 1.01178191349e+21 , 0.0 , 0.0) 838 839, step length=0.0258511490836 840 at ( - 1.25031667658 , 0.0 , 0.0 , 98.7496833234 , 100.0 , 100.0) 841 842*** evaluate Jacobian (or its inverse) 843 844*** compute the next point 845 846*** evalute function in the next point 847 848*** compute the point difference 849 850*** compute the size of the point difference 851 85230. residue=(3.72213764917e+20 , 0.0 , 0.0 , 3.72213764921e+20 , 0.0 , 0.0) 853 854, step length=0.0258511490836 855 at ( - 1.22446552749 , 0.0 , 0.0 , 98.7755344725 , 100.0 , 100.0) 856 857*** evaluate Jacobian (or its inverse) 858 859*** compute the next point 860 861*** evalute function in the next point 862 863*** compute the point difference 864 865*** compute the size of the point difference 866 86731. residue=(1.36929791834e+20 , 0.0 , 0.0 , 1.36929791835e+20 , 0.0 , 0.0) 868 869, step length=0.0258511490836 870 at ( - 1.19861437841 , 0.0 , 0.0 , 98.8013856216 , 100.0 , 100.0) 871 872*** evaluate Jacobian (or its inverse) 873 874*** compute the next point 875 876*** evalute function in the next point 877 878*** compute the point difference 879 880*** compute the size of the point difference 881 88232. residue=(5.03736552996e+19 , 0.0 , 0.0 , 5.03736553001e+19 , 0.0 , 0.0) 883 884, step length=0.0258511490836 885 at ( - 1.17276322933 , 0.0 , 0.0 , 98.8272367707 , 100.0 , 100.0) 886 887*** evaluate Jacobian (or its inverse) 888 889*** compute the next point 890 891*** evalute function in the next point 892 893*** compute the point difference 894 895*** compute the size of the point difference 896 89733. residue=(1.85314321614e+19 , 0.0 , 0.0 , 1.85314321616e+19 , 0.0 , 0.0) 898 899, step length=0.0258511490836 900 at ( - 1.14691208024 , 0.0 , 0.0 , 98.8530879198 , 100.0 , 100.0) 901 902*** evaluate Jacobian (or its inverse) 903 904*** compute the next point 905 906*** evalute function in the next point 907 908*** compute the point difference 909 910*** compute the size of the point difference 911 91234. residue=(6.81733290764e+18 , 0.0 , 0.0 , 6.81733290771e+18 , 0.0 , 0.0) 913 914, step length=0.0258511490836 915 at ( - 1.12106093116 , 0.0 , 0.0 , 98.8789390688 , 100.0 , 100.0) 916 917*** evaluate Jacobian (or its inverse) 918 919*** compute the next point 920 921*** evalute function in the next point 922 923*** compute the point difference 924 925*** compute the size of the point difference 926 92735. residue=(2.50795662034e+18 , 0.0 , 0.0 , 2.50795662037e+18 , 0.0 , 0.0) 928 929, step length=0.0258511490836 930 at ( - 1.09520978207 , 0.0 , 0.0 , 98.9047902179 , 100.0 , 100.0) 931 932*** evaluate Jacobian (or its inverse) 933 934*** compute the next point 935 936*** evalute function in the next point 937 938*** compute the point difference 939 940*** compute the size of the point difference 941 94236. residue=(9.2262567997e+17 , 0.0 , 0.0 , 9.22625679991e+17 , 0.0 , 0.0) 943 944, step length=0.0258511490837 945 at ( - 1.06935863299 , 0.0 , 0.0 , 98.930641367 , 100.0 , 100.0) 946 947*** evaluate Jacobian (or its inverse) 948 949*** compute the next point 950 951*** evalute function in the next point 952 953*** compute the point difference 954 955*** compute the size of the point difference 956 95737. residue=(3.39415019556e+17 , 0.0 , 0.0 , 3.39415019568e+17 , 0.0 , 0.0) 958 959, step length=0.0258511490838 960 at ( - 1.04350748391 , 0.0 , 0.0 , 98.9564925161 , 100.0 , 100.0) 961 962*** evaluate Jacobian (or its inverse) 963 964*** compute the next point 965 966*** evalute function in the next point 967 968*** compute the point difference 969 970*** compute the size of the point difference 971 97238. residue=(1.24863807717e+17 , 0.0 , 0.0 , 1.24863807725e+17 , 0.0 , 0.0) 973 974, step length=0.0258511490842 975 at ( - 1.01765633483 , 0.0 , 0.0 , 98.9823436652 , 100.0 , 100.0) 976 977*** evaluate Jacobian (or its inverse) 978 979*** compute the next point 980 981*** evalute function in the next point 982 983*** compute the point difference 984 985*** compute the size of the point difference 986 98739. residue=(4.59348278034e+16 , 0.0 , 0.0 , 4.59348278107e+16 , 0.0 , 0.0) 988 989, step length=0.0258511490853 990 at ( - 0.991805185743 , 0.0 , 0.0 , 99.0081948143 , 100.0 , 100.0) 991 992*** evaluate Jacobian (or its inverse) 993 994*** compute the next point 995 996*** evalute function in the next point 997 998*** compute the point difference 999 1000*** compute the size of the point difference 1001 100240. residue=(1.68984787804e+16 , 0.0 , 0.0 , 1.68984787874e+16 , 0.0 , 0.0) 1003 1004, step length=0.0258511490882 1005 at ( - 0.965954036664 , 0.0 , 0.0 , 99.0340459633 , 100.0 , 100.0) 1006 1007*** evaluate Jacobian (or its inverse) 1008 1009*** compute the next point 1010 1011*** evalute function in the next point 1012 1013*** compute the point difference 1014 1015*** compute the size of the point difference 1016 101741. residue=(6.21660292823e+15 , 0.0 , 0.0 , 6.21660293516e+15 , 0.0 , 0.0) 1018 1019, step length=0.0258511490961 1020 at ( - 0.940102887593 , 0.0 , 0.0 , 99.0598971124 , 100.0 , 100.0) 1021 1022*** evaluate Jacobian (or its inverse) 1023 1024*** compute the next point 1025 1026*** evalute function in the next point 1027 1028*** compute the point difference 1029 1030*** compute the size of the point difference 1031 103242. residue=(2.28696040906e+15 , 0.0 , 0.0 , 2.28696041594e+15 , 0.0 , 0.0) 1033 1034, step length=0.0258511491177 1035 at ( - 0.914251738544 , 0.0 , 0.0 , 99.0857482616 , 100.0 , 100.0) 1036 1037*** evaluate Jacobian (or its inverse) 1038 1039*** compute the next point 1040 1041*** evalute function in the next point 1042 1043*** compute the point difference 1044 1045*** compute the size of the point difference 1046 104743. residue=(8.41325715099e+14 , 0.0 , 0.0 , 8.41325721962e+14 , 0.0 , 0.0) 1048 1049, step length=0.0258511491762 1050 at ( - 0.888400589553 , 0.0 , 0.0 , 99.1115994107 , 100.0 , 100.0) 1051 1052*** evaluate Jacobian (or its inverse) 1053 1054*** compute the next point 1055 1056*** evalute function in the next point 1057 1058*** compute the point difference 1059 1060*** compute the size of the point difference 1061 106244. residue=(3.09506431748e+14 , 0.0 , 0.0 , 3.09506438604e+14 , 0.0 , 0.0) 1063 1064, step length=0.0258511493354 1065 at ( - 0.862549440721 , 0.0 , 0.0 , 99.1374505601 , 100.0 , 100.0) 1066 1067*** evaluate Jacobian (or its inverse) 1068 1069*** compute the next point 1070 1071*** evalute function in the next point 1072 1073*** compute the point difference 1074 1075*** compute the size of the point difference 1076 107745. residue=(1.13861050984e+14 , 0.0 , 0.0 , 1.13861057839e+14 , 0.0 , 0.0) 1078 1079, step length=0.0258511497682 1080 at ( - 0.836698292322 , 0.0 , 0.0 , 99.1633017098 , 100.0 , 100.0) 1081 1082*** evaluate Jacobian (or its inverse) 1083 1084*** compute the next point 1085 1086*** evalute function in the next point 1087 1088*** compute the point difference 1089 1090*** compute the size of the point difference 1091 109246. residue=(4.18871376415e+13 , 0.0 , 0.0 , 4.18871444947e+13 , 0.0 , 0.0) 1093 1094, step length=0.0258511509446 1095 at ( - 0.8108471451 , 0.0 , 0.0 , 99.1891528608 , 100.0 , 100.0) 1096 1097*** evaluate Jacobian (or its inverse) 1098 1099*** compute the next point 1100 1101*** evalute function in the next point 1102 1103*** compute the point difference 1104 1105*** compute the size of the point difference 1106 110747. residue=(1.54094146219e+13 , 0.0 , 0.0 , 1.54094214749e+13 , 0.0 , 0.0) 1108 1109, step length=0.0258511541423 1110 at ( - 0.784996001075 , 0.0 , 0.0 , 99.2150040149 , 100.0 , 100.0) 1111 1112*** evaluate Jacobian (or its inverse) 1113 1114*** compute the next point 1115 1116*** evalute function in the next point 1117 1118*** compute the point difference 1119 1120*** compute the size of the point difference 1121 112248. residue=(5.66880467397e+12 , 0.0 , 0.0 , 5.6688115269e+12 , 0.0 , 0.0) 1123 1124, step length=0.0258511628346 1125 at ( - 0.759144865742 , 0.0 , 0.0 , 99.2408551778 , 100.0 , 100.0) 1126 1127*** evaluate Jacobian (or its inverse) 1128 1129*** compute the next point 1130 1131*** evalute function in the next point 1132 1133*** compute the point difference 1134 1135*** compute the size of the point difference 1136 113749. residue=(2.08543452966e+12 , 0.0 , 0.0 , 2.08544138253e+12 , 0.0 , 0.0) 1138 1139, step length=0.0258511864627 1140 at ( - 0.733293754037 , 0.0 , 0.0 , 99.2667063642 , 100.0 , 100.0) 1141 1142*** evaluate Jacobian (or its inverse) 1143 1144*** compute the next point 1145 1146*** evalute function in the next point 1147 1148*** compute the point difference 1149 1150*** compute the size of the point difference 1151 115250. residue=(767186323467.0 , 0.0 , 0.0 , 767193176319.0 , 0.0 , 0.0) 1153 1154, step length=0.0258512506906 1155 at ( - 0.70744270656 , 0.0 , 0.0 , 99.2925576149 , 100.0 , 100.0) 1156 1157*** evaluate Jacobian (or its inverse) 1158 1159*** compute the next point 1160 1161*** evalute function in the next point 1162 1163*** compute the point difference 1164 1165*** compute the size of the point difference 1166 116751. residue=(282229910057.0 , 0.0 , 0.0 , 282236762901.0 , 0.0 , 0.0) 1168 1169, step length=0.0258514252812 1170 at ( - 0.681591833671 , 0.0 , 0.0 , 99.3184090402 , 100.0 , 100.0) 1171 1172*** evaluate Jacobian (or its inverse) 1173 1174*** compute the next point 1175 1176*** evalute function in the next point 1177 1178*** compute the point difference 1179 1180*** compute the size of the point difference 1181 118252. residue=(103824415727.0 , 0.0 , 0.0 , 103831268569.0 , 0.0 , 0.0) 1183 1184, step length=0.0258518998746 1185 at ( - 0.655741435353 , 0.0 , 0.0 , 99.3442609401 , 100.0 , 100.0) 1186 1187*** evaluate Jacobian (or its inverse) 1188 1189*** compute the next point 1190 1191*** evalute function in the next point 1192 1193*** compute the point difference 1194 1195*** compute the size of the point difference 1196 119753. residue=(3.81927022457e+10 , 0.0 , 0.0 , 3.8199555087e+10 , 0.0 , 0.0) 1198 1199, step length=0.0258531900044 1200 at ( - 0.629892327003 , 0.0 , 0.0 , 99.3701141301 , 100.0 , 100.0) 1201 1202*** evaluate Jacobian (or its inverse) 1203 1204*** compute the next point 1205 1206*** evalute function in the next point 1207 1208*** compute the point difference 1209 1210*** compute the size of the point difference 1211 121254. residue=(1.40481443717e+10 , 0.0 , 0.0 , 1.40549972128e+10 , 0.0 , 0.0) 1213 1214, step length=0.0258566973195 1215 at ( - 0.604046724769 , 0.0 , 0.0 , 99.3959708274 , 100.0 , 100.0) 1216 1217*** evaluate Jacobian (or its inverse) 1218 1219*** compute the next point 1220 1221*** evalute function in the next point 1222 1223*** compute the point difference 1224 1225*** compute the size of the point difference 1226 122755. residue=(5.16585846951e+9 , 0.0 , 0.0 , 5.17271131074e+9 , 0.0 , 0.0) 1228 1229, step length=0.0258662339895 1230 at ( - 0.578210650353 , 0.0 , 0.0 , 99.4218370614 , 100.0 , 100.0) 1231 1232*** evaluate Jacobian (or its inverse) 1233 1234*** compute the next point 1235 1236*** evalute function in the next point 1237 1238*** compute the point difference 1239 1240*** compute the size of the point difference 1241 124256. residue=(1.89824960589e+9 , 0.0 , 0.0 , 1.90510244878e+9 , 0.0 , 0.0) 1243 1244, step length=0.025892178044 1245 at ( - 0.552400454575 , 0.0 , 0.0 , 99.4477292394 , 100.0 , 100.0) 1246 1247*** evaluate Jacobian (or its inverse) 1248 1249*** compute the next point 1250 1251*** evalute function in the next point 1252 1253*** compute the point difference 1254 1255*** compute the size of the point difference 1256 125757. residue=(6.96167585052e+8 , 0.0 , 0.0 , 7.03020440594e+8 , 0.0 , 0.0) 1258 1259, step length=0.0259628545108 1260 at ( - 0.526660451901 , 0.0 , 0.0 , 99.4736920939 , 100.0 , 100.0) 1261 1262*** evaluate Jacobian (or its inverse) 1263 1264*** compute the next point 1265 1266*** evalute function in the next point 1267 1268*** compute the point difference 1269 1270*** compute the size of the point difference 1271 127258. residue=(2.53957444815e+8 , 0.0 , 0.0 , 2.60810394004e+8 , 0.0 , 0.0) 1273 1274, step length=0.026156110844 1275 at ( - 0.501110133883 , 0.0 , 0.0 , 99.4998482048 , 100.0 , 100.0) 1276 1277*** evaluate Jacobian (or its inverse) 1278 1279*** compute the next point 1280 1281*** evalute function in the next point 1282 1283*** compute the point difference 1284 1285*** compute the size of the point difference 1286 128759. residue=(9.13074482936e+7 , 0.0 , 0.0 , 9.81610893493e+7 , 0.0 , 0.0) 1288 1289, step length=0.0266899582516 1290 at ( - 0.476067267452 , 0.0 , 0.0 , 99.526538163 , 100.0 , 100.0) 1291 1292*** evaluate Jacobian (or its inverse) 1293 1294*** compute the next point 1295 1296*** evalute function in the next point 1297 1298*** compute the point difference 1299 1300*** compute the size of the point difference 1301 130260. residue=(3.15519019482e+7 , 0.0 , 0.0 , 3.8410646839e+7 , 0.0 , 0.0) 1303 1304, step length=0.0282064667415 1305 at ( - 0.452345623749 , 0.0 , 0.0 , 99.5547446298 , 100.0 , 100.0) 1306 1307*** evaluate Jacobian (or its inverse) 1308 1309*** compute the next point 1310 1311*** evalute function in the next point 1312 1313*** compute the point difference 1314 1315*** compute the size of the point difference 1316 131761. residue=(9.77481469301e+6 , 0.0 , 0.0 , 1.66708124709e+7 , 0.0 , 0.0) 1318 1319, step length=0.0328642948738 1320 at ( - 0.431825342687 , 0.0 , 0.0 , 99.5876089247 , 100.0 , 100.0) 1321 1322*** evaluate Jacobian (or its inverse) 1323 1324*** compute the next point 1325 1326*** evalute function in the next point 1327 1328*** compute the point difference 1329 1330*** compute the size of the point difference 1331 133262. residue=(2.23533931681e+6 , 0.0 , 0.0 , 9.38172386018e+6 , 0.0 , 0.0) 1333 1334, step length=0.0508561508748 1335 at ( - 0.417764764235 , 0.0 , 0.0 , 99.6384650755 , 100.0 , 100.0) 1336 1337*** evaluate Jacobian (or its inverse) 1338 1339*** compute the next point 1340 1341*** evalute function in the next point 1342 1343*** compute the point difference 1344 1345*** compute the size of the point difference 1346 134763. residue=(2.23262484734e+5 , 0.0 , 0.0 , 8.19715321429e+6 , 0.0 , 0.0) 1348 1349, step length=0.20466482746 1350 at ( - 0.41222548566 , 0.0 , 0.0 , 99.843129903 , 100.0 , 100.0) 1351 1352*** evaluate Jacobian (or its inverse) 1353 1354*** compute the next point 1355 1356*** evalute function in the next point 1357 1358*** compute the point difference 1359 1360*** compute the size of the point difference 1361 1362*** evalute function in the next point 1363 1364*** compute the point difference 1365 1366*** compute the size of the point difference 1367 1368*** reduce the difference limititeration to its half 1369 1370*** evalute function in the next point 1371 1372*** compute the point difference 1373 1374*** compute the size of the point difference 1375 1376*** reduce the difference limititeration to its half 1377 1378*** evalute function in the next point 1379 1380*** compute the point difference 1381 1382*** compute the size of the point difference 1383 1384*** reduce the difference limititeration to its half 1385 1386*** evalute function in the next point 1387 1388*** compute the point difference 1389 1390*** compute the size of the point difference 1391 1392*** reduce the difference limititeration to its half 1393 1394*** evalute function in the next point 1395 1396*** compute the point difference 1397 1398*** compute the size of the point difference 1399 1400*** reduce the difference limititeration to its half 1401 1402*** evalute function in the next point 1403 1404*** compute the point difference 1405 1406*** compute the size of the point difference 1407 1408*** reduce the difference limititeration to its half 1409 1410*** evalute function in the next point 1411 1412*** compute the point difference 1413 1414*** compute the size of the point difference 1415 1416*** reduce the difference limititeration to its half 1417 1418*** evalute function in the next point 1419 1420*** compute the point difference 1421 1422*** compute the size of the point difference 1423 142464. residue=(2.23088062724e+5 , 0.0 , 0.0 , 8.19035290355e+6 , 0.0 , 0.0) 1425 1426, step length=0.383303021247 1427 at ( - 0.412224950141 , 0.0 , 0.0 , 100.226432924 , 100.0 , 100.0) 1428 1429*** evaluate Jacobian (or its inverse) 1430 1431*** compute the next point 1432 1433*** evalute function in the next point 1434 1435*** compute the point difference 1436 1437*** compute the size of the point difference 1438 1439*** reduce the difference limititeration to its half 1440 1441*** evalute function in the next point 1442 1443*** compute the point difference 1444 1445*** compute the size of the point difference 1446 1447*** reduce the difference limititeration to its half 1448 1449*** evalute function in the next point 1450 1451*** compute the point difference 1452 1453*** compute the size of the point difference 1454 1455*** reduce the difference limititeration to its half 1456 1457*** evalute function in the next point 1458 1459*** compute the point difference 1460 1461*** compute the size of the point difference 1462 1463*** reduce the difference limititeration to its half 1464 1465*** evalute function in the next point 1466 1467*** compute the point difference 1468 1469*** compute the size of the point difference 1470 1471*** reduce the difference limititeration to its half 1472 1473*** evalute function in the next point 1474 1475*** compute the point difference 1476 1477*** compute the size of the point difference 1478 1479*** reduce the difference limititeration to its half 1480 1481*** evalute function in the next point 1482 1483*** compute the point difference 1484 1485*** compute the size of the point difference 1486 1487*** reduce the difference limititeration to its half 1488 1489*** evalute function in the next point 1490 1491*** compute the point difference 1492 1493*** compute the size of the point difference 1494 1495*** reduce the difference limititeration to its half 1496 1497*** evalute function in the next point 1498 1499*** compute the point difference 1500 1501*** compute the size of the point difference 1502 150365. residue=(2.22216670074e+5 , 0.0 , 0.0 , 7.2288078e+6 , 0.0 , 0.0) 1504 1505, step length=0.129870827146 1506 at ( - 0.412222274586 , 0.0 , 0.0 , 100.356303751 , 100.0 , 100.0) 1507 1508*** evaluate Jacobian (or its inverse) 1509 1510*** compute the next point 1511 1512*** evalute function in the next point 1513 1514*** compute the point difference 1515 1516*** compute the size of the point difference 1517 1518*** reduce the difference limititeration to its half 1519 1520*** evalute function in the next point 1521 1522*** compute the point difference 1523 1524*** compute the size of the point difference 1525 1526*** reduce the difference limititeration to its half 1527 1528*** evalute function in the next point 1529 1530*** compute the point difference 1531 1532*** compute the size of the point difference 1533 153466. residue=(1.66845393097e+5 , 0.0 , 0.0 , 1.93472870727e+6 , 0.0 , 0.0) 1535 1536, step length=0.0482669644337 1537 at ( - 0.412051690235 , 0.0 , 0.0 , 100.404570716 , 100.0 , 100.0) 1538 1539*** evaluate Jacobian (or its inverse) 1540 1541*** compute the next point 1542 1543*** evalute function in the next point 1544 1545*** compute the point difference 1546 1547*** compute the size of the point difference 1548 154967. residue=(1653.19388834 , 0.0 , 0.0 , - 3.3219398641e+5 , 0.0 , 0.0) 1550 1551, step length=0.00798706792058 1552 at ( - 0.411535983805 , 0.0 , 0.0 , 100.412557784 , 100.0 , 100.0) 1553 1554*** evaluate Jacobian (or its inverse) 1555 1556*** compute the next point 1557 1558*** evalute function in the next point 1559 1560*** compute the point difference 1561 1562*** compute the size of the point difference 1563 156468. residue=(0.166671264823 , 0.0 , 0.0 , - 6386.15622619 , 0.0 , 0.0) 1565 1566, step length=0.00100688023827 1567 at ( - 0.411530770947 , 0.0 , 0.0 , 100.411550903 , 100.0 , 100.0) 1568 1569*** evaluate Jacobian (or its inverse) 1570 1571*** compute the next point 1572 1573*** evalute function in the next point 1574 1575*** compute the point difference 1576 1577*** compute the size of the point difference 1578 157969. residue=( - 0.000000122 , 0.0 , 0.0 , - 2.48519530414 , 0.0 , 0.0) 1580 1581, step length=0.00002012523636 1582 at ( - 0.411530770421 , 0.0 , 0.0 , 100.411530778 , 100.0 , 100.0) 1583 1584{x1= - 0.411530770421,x2=0.0,x3=0.0,x4=100.41153077,x5=100.0,x6=100.0} 1585 1586off trnumeric; 1587 1588 1589clear alpha,ni,v,d,sys; 1590 1591 1592off rounded; 1593 1594 1595 1596% INTEGRALS 1597 1598num_int( x**2,x=(1 .. 2),accuracy=3); 1599 1600 16012.33333333333 1602 1603 1604 % 1st case: using formal integral 1605needle := 1/(10**-4 + x**2); 1606 1607 1608 10000 1609needle := -------------- 1610 2 1611 10000*x + 1 1612 1613num_int(needle,x=(-1 .. 1),accuracy=3); 1614 1615 1616312.159332022 1617 % 312.16 1618 1619 % no formal integral, but easy Chebyshev fit 1620num_int(sin x/x,x=(1 .. 10)); 1621 1622 16230.712264523852 1624 1625 1626 % using a Chebyshev fit of order 60 1627num_int(exp(-x**2),x=(-10 .. 10),accuracy=3); 1628 1629 16301.77245387654 1631 % 1.772 1632 1633 % cases with singularities 1634 1635num_int(1/sqrt x ,x=(0 .. 1),accuracy=2); 1636 1637 16382 1639 % 1.999 1640 1641num_int(1/sqrt abs x ,x=(-1 .. 1),iterations=50); 1642 1643 16443.99999231465 1645 % 3.999 1646 1647 % simple multidimensional integrals 1648num_int(x+y,x=(0 .. 1),y=(2 .. 3)); 1649 1650 16513.0 1652 1653 1654num_int(sin(x+y),x=(0 .. 1),y=(0 .. 1)); 1655 1656 16570.773135425204 1658 1659 1660% some integrals with infinite bounds 1661 1662on rounded; 1663 1664 % for the error function 1665 1666num_int(e^(-x) ,x=(0 .. infinity)); 1667 1668 16691.00000034605 1670 % 1.000 1671 16722/sqrt(pi)* num_int(e^(-x^2) ,x=(0 .. infinity)); 1673 1674 16751.00000003784 1676 % 1.00 1677 16782/sqrt(pi)* num_int(e^(-x^2), x=(-infinity .. infinity)); 1679 1680 16812.00000007569 1682 % 2.00 1683 1684num_int(sin(x) * e^(-x), x=(0 .. infinity)); 1685 1686 16870.500000522701 1688 % 0.500 1689 1690off rounded; 1691 1692 1693 1694% APPROXIMATION 1695 1696 %approximate sin x by a cubic polynomial 1697num_fit(sin x,{1,x,x**2,x**3},x=for i:=0:20 collect 0.1*i); 1698 1699 1700 3 2 1701{ - 0.0847539694989*x - 0.134641944765*x + 1.06263064633*x - 0.00519313406462, 1702 1703 { - 0.00519313406462,1.06263064633, - 0.134641944765, - 0.0847539694989}} 1704 1705 1706 % approximate x**2 by a harmonic series in the interval [0,1] 1707num_fit(x**2,1 . for i:=1:5 join {sin(i*x)}, 1708 x=for i:=0:10 collect i/10); 1709 1710 1711{ - 1.30957809531*sin(5*x) + 7.16375571726*sin(4*x) - 18.5490183807*sin(3*x) 1712 1713 + 26.5601711119*sin(2*x) - 19.4492187125*sin(x) - 0.00197199703147, 1714 1715 { - 0.00197199703147, - 19.4492187125,26.5601711119, - 18.5490183807, 1716 1717 7.16375571726, - 1.30957809531}} 1718 1719 1720 % approximate a set of points by a polynomial 1721pts:=for i:=1 step 0.1 until 3 collect i$ 1722 1723 1724vals:=for each p in pts collect (p+2)**3$ 1725 1726 1727num_fit(vals,{1,x,x**2,x**3},x=pts); 1728 1729 1730 3 2 1731{1.00000000002*x + 5.99999999991*x + 12.0000000002*x + 7.99999999989,{ 1732 1733 7.99999999989,12.0000000002,5.99999999991,1.00000000002}} 1734 1735 % compute the approximation error 1736on rounded; 1737 1738 1739first ws - (x+2)**3; 1740 1741 1742 3 2 17430.0000000000155546686642*x - 0.0000000000946407396896*x 1744 1745 + 0.00000000018181722794*x - 0.000000000109046105479 1746 1747off rounded; 1748 1749 1750 1751 1752 1753% ODE SOLUTION (Runge-Kutta) 1754 1755depend(y,x); 1756 1757 1758 1759 % approximate y=y(x) with df(y,x)=2y in interval [0 : 5] 1760num_odesolve(df(y,x)=y,y=2,x=(0 .. 5),iterations=20); 1761 1762 1763{{x,y}, 1764 1765 {0.0,2.0}, 1766 1767 {0.25,2.56805083337}, 1768 1769 {0.5,3.2974425414}, 1770 1771 {0.75,4.23400003322}, 1772 1773 {1.0,5.43656365691}, 1774 1775 {1.25,6.98068591491}, 1776 1777 {1.5,8.96337814065}, 1778 1779 {1.75,11.509205352}, 1780 1781 {2.0,14.7781121978}, 1782 1783 {2.25,18.9754716726}, 1784 1785 {2.5,24.3649879213}, 1786 1787 {2.75,31.2852637682}, 1788 1789 {3.0,40.1710738461}, 1790 1791 {3.25,51.5806798341}, 1792 1793 {3.5,66.2309039169}, 1794 1795 {3.75,85.0421639995}, 1796 1797 {4.0,109.196300065}, 1798 1799 {4.25,140.210824692}, 1800 1801 {4.5,180.034262599}, 1802 1803 {4.75,231.168569052}, 1804 1805 {5.0,296.826318202}} 1806 1807 1808 % same with negative direction 1809num_odesolve(df(y,x)=y,y=2,x=(0 .. -5),iterations=20); 1810 1811 1812{{x,y}, 1813 1814 {0.0,2.0}, 1815 1816 {-0.25,1.55760156614}, 1817 1818 {-0.5,1.21306131943}, 1819 1820 {-0.75,0.944733105483}, 1821 1822 {-1.0,0.735758882344}, 1823 1824 {-1.25,0.573009593722}, 1825 1826 {-1.5,0.446260320298}, 1827 1828 {-1.75,0.347547886902}, 1829 1830 {-2.0,0.270670566474}, 1831 1832 {-2.25,0.210798449125}, 1833 1834 {-2.5,0.164169997249}, 1835 1836 {-2.75,0.127855722414}, 1837 1838 {-3.0,0.0995741367363}, 1839 1840 {-3.25,0.0775484156639}, 1841 1842 {-3.5,0.060394766845}, 1843 1844 {-3.75,0.0470354917124}, 1845 1846 {-4.0,0.0366312777778}, 1847 1848 {-4.25,0.0285284678182}, 1849 1850 {-4.5,0.0222179930767}, 1851 1852 {-4.75,0.0173033904064}, 1853 1854 {-5.0,0.0134758939983}} 1855 1856 1857 % giving a nice picture when plotted 1858num_odesolve(df(y,x)=1- x*y**2 ,y=0,x=(0 .. 4),iterations=20); 1859 1860 1861{{x,y}, 1862 1863 {0.0,0.0}, 1864 1865 {0.2,0.199600912188}, 1866 1867 {0.4,0.393714914166}, 1868 1869 {0.6,0.569482634406}, 1870 1871 {0.8,0.710657687564}, 1872 1873 {1.0,0.805480022354}, 1874 1875 {1.2,0.852604291055}, 1876 1877 {1.4,0.860563377356}, 1878 1879 {1.6,0.842333334456}, 1880 1881 {1.8,0.80999200878}, 1882 1883 {2.0,0.772211952811}, 1884 1885 {2.2,0.734163640068}, 1886 1887 {2.4,0.698433235122}, 1888 1889 {2.6,0.666019196492}, 1890 1891 {2.8,0.637070046905}, 1892 1893 {3.0,0.611341375657}, 1894 1895 {3.2,0.588447372601}, 1896 1897 {3.4,0.567985133759}, 1898 1899 {3.6,0.549587947292}, 1900 1901 {3.8,0.532942255624}, 1902 1903 {4.0,0.517787833735}} 1904 1905 1906 % system of ordinary differential equations 1907depend(y,x); 1908 1909 1910depend(z,x); 1911 1912 1913num_odesolve( 1914 {df(z,x) = y, df(y,x)= y+x}, 1915 {z=2, y=4}, 1916 x=(0 .. 5),iterations=20); 1917 1918 1919{{x,z,y}, 1920 1921 {0.0,2.0,4.0}, 1922 1923 {0.25,3.13887708344,5.17012708344}, 1924 1925 {0.5,4.61860635349,6.74360635349}, 1926 1927 {0.75,6.55375008305,8.83500008305}, 1928 1929 {1.0,9.09140914227,11.5914091423}, 1930 1931 {1.25,12.4204647873,15.2017147873}, 1932 1933 {1.5,16.7834453516,19.9084453516}, 1934 1935 {1.75,22.4917633799,26.0230133799}, 1936 1937 {2.0,29.9452804945,33.9452804945}, 1938 1939 {2.25,39.6574291816,44.1886791816}, 1940 1941 {2.5,52.2874698032,57.4124698032}, 1942 1943 {2.75,68.6819094205,74.4631594205}, 1944 1945 {3.0,89.9276846154,96.4276846154}, 1946 1947 {3.25,117.420449585,124.701699585}, 1948 1949 {3.5,152.952259792,161.077259792}, 1950 1951 {3.75,198.824159999,207.855409999}, 1952 1953 {4.0,257.990750164,267.990750164}, 1954 1955 {4.25,334.245811731,345.277061731}, 1956 1957 {4.5,432.460656499,444.585656499}, 1958 1959 {4.75,558.890172631,572.171422631}, 1960 1961 {5.0,721.565795506,736.065795506}} 1962 1963 1964 1965%----------------- Chebyshev fit ------------------------- 1966 1967on rounded; 1968 1969 1970 1971func := x**2 * (x**2 - 2) * sin x; 1972 1973 1974 2 2 1975func := sin(x)*x *(x - 2) 1976 1977ord := 15; 1978 1979 1980ord := 15 1981 1982 1983cx:=chebyshev_fit(func,x=(0 .. 2),ord)$ 1984 1985 1986cp:=first cx; 1987 1988 1989 13 12 11 1990cp := 0.000000620105096185*x + 0.0000168737305191*x - 0.000269014332288*x 1991 1992 10 9 8 1993 + 0.000155646029006*x + 0.00848163760265*x + 0.000272748604876*x 1994 1995 7 6 5 1996 - 0.183540091904*x + 0.00010680840443*x + 1.33329694616*x 1997 1998 4 3 2 1999 + 0.00000770692780683*x - 2.00000091554*x + 0.0000000501515695639*x 2000 2001 - 0.000000000784989184766*x - 4.86055640181e-13 2002 2003cc:=second cx; 2004 2005 2006cc := {2.69320512829, 2007 2008 2.76751928466, 2009 2010 2.25642507569, 2011 2012 0.955452569949, 2013 2014 0.0509075944268, 2015 2016 - 0.0868248678183, 2017 2018 - 0.0170919216091, 2019 2020 0.00104527137626, 2021 2022 0.000349190502034, 2023 2024 - 0.00000253521592323, 2025 2026 - 0.00000280798840641, 2027 2028 - 0.0000000157676044858, 2029 2030 0.0000000121753402195, 2031 2032 0.000000000118269801846, 2033 2034 - 0.0000000000331230439026} 2035 2036 2037for u:=0 step 0.2 until 2 do write 2038 "x:",u," true value:",sub(x=u,func), 2039 " Chebyshev eval:", chebyshev_eval(cc,x=(0 .. 2),x=u), 2040 " Chebyshev polynomial:",sub(x=u,cp); 2041 2042 2043x:0 true value:0 Chebyshev eval: - 4.85167461761e-13 Chebyshev polynomial: 2044 2045 - 4.86055640181e-13 2046 2047x:0.2 true value: - 0.0155756755343 Chebyshev eval: - 0.0155756755339 2048 2049 Chebyshev polynomial: - 0.015575675548 2050 2051x:0.4 true value: - 0.114644759976 Chebyshev eval: - 0.114644759976 2052 2053 Chebyshev polynomial: - 0.114644759974 2054 2055x:0.6 true value: - 0.333364916292 Chebyshev eval: - 0.333364916292 2056 2057 Chebyshev polynomial: - 0.333364916295 2058 2059x:0.8 true value: - 0.624386741519 Chebyshev eval: - 0.624386741519 2060 2061 Chebyshev polynomial: - 0.624386741504 2062 2063x:1 true value: - 0.841470984808 Chebyshev eval: - 0.841470984808 2064 2065 Chebyshev polynomial: - 0.841470984841 2066 2067x:1.2 true value: - 0.751596318924 Chebyshev eval: - 0.751596318924 2068 2069 Chebyshev polynomial: - 0.751596318876 2070 2071x:1.4 true value: - 0.0772592588311 Chebyshev eval: - 0.0772592588311 2072 2073 Chebyshev polynomial: - 0.0772592588864 2074 2075x:1.6 true value:1.43298871732 Chebyshev eval:1.43298871732 2076 2077 Chebyshev polynomial:1.43298871738 2078 2079x:1.8 true value:3.91253024182 Chebyshev eval:3.91253024182 2080 2081 Chebyshev polynomial:3.91253024177 2082 2083x:2.0 true value:7.27437941461 Chebyshev eval:7.27437941461 2084 2085 Chebyshev polynomial:7.27437941467 2086 2087 2088% integral 2089 2090 % integrate coefficients 2091ci := chebyshev_int(cc,x=(0 .. 2)); 2092 2093 2094ci := {0.0310113015322, 2095 2096 0.2183900263, 2097 2098 0.453016678678, 2099 2100 0.367586246877, 2101 2102 0.130284679721, 2103 2104 0.00679995160359, 2105 2106 - 0.00732251159954, 2107 2108 - 0.00124579372222, 2109 2110 0.0000654879120115, 2111 2112 0.0000195554716911, 2113 2114 - 0.000000125972415937, 2115 2116 - 0.000000128189261211, 2117 2118 - 0.000000000661911428653, 2119 2120 0.000000000469556279362, 2121 2122 4.22392149449e-12} 2123 2124 2125 % compare with true values (normalized absolute term) 2126ci0:=chebyshev_eval(ci,x=(0 .. 2),x=0)$ 2127 2128 2129ifunc := int(func,x)$ 2130 2131 if0 := sub(x=0,ifunc); 2132 2133 2134if0 := - 28.0 2135 2136 2137for u:=0 step 0.2 until 2 do write 2138 {u,sub(x=u,ifunc) - if0, 2139 chebyshev_eval(ci,x=(0 .. 2),x=u) - ci0}; 2140 2141 2142{0,0,0} 2143 2144{0.2, - 0.000785836355117, - 0.00078583635293} 2145 2146{0.4, - 0.0119047051867, - 0.0119047051858} 2147 2148{0.6, - 0.0548116700418, - 0.0548116700408} 2149 2150{0.8, - 0.150297976106, - 0.150297976105} 2151 2152{1, - 0.299838223412, - 0.29983822341} 2153 2154{1.2, - 0.466528961073, - 0.466528961072} 2155 2156{1.4, - 0.561460555384, - 0.561460555383} 2157 2158{1.6, - 0.441445769516, - 0.441445769514} 2159 2160{1.8,0.0768452822437,0.0768452822437} 2161 2162{2.0,1.18309971762,1.18309971762} 2163 2164 2165% derivative 2166 2167 % differentiate coefficients 2168cd := chebyshev_df(cc,x=(0 .. 2))$ 2169 2170 2171 % compute coefficients of derivative 2172cds := second chebyshev_fit(df(func,x),x=(0 .. 2),ord)$ 2173 2174 2175 % compare coefficients 2176for i:=1:ord do write {part(cd,i),part(cds,i)}; 2177 2178 2179{10.4140931324,10.4140931324} 2180 2181{9.23338917839,9.2333891784} 2182 2183{4.87905456308,4.87905456307} 2184 2185{0.207688875651,0.207688875654} 2186 2187{ - 0.853660856614, - 0.853660856625} 2188 2189{ - 0.199571879764, - 0.19957187976} 2190 2191{0.0145878215687,0.0145878215579} 2192 2193{0.00553117954514,0.00553117954883} 2194 2195{ - 0.000045977698902, - 0.0000459777097756} 2196 2197{ - 0.0000558684874082, - 0.0000558684837245} 2198 2199{ - 0.00000034381228384, - 0.00000034382315144} 2200 2201{0.000000291280720039,0.00000029128440905} 2202 2203{0.00000000307501484799,0.00000000306414359587} 2204 2205{ - 0.000000000927445229273, - 0.000000000923750392908} 2206 2207{0, - 0.0000000000109242472928} 2208 2209 2210clear func,ord,cc,cx,cd,cds,ci,ci0; 2211 2212 2213 2214% One from ISSAC '97 -- should be ~ 1.10*10^36300 2215 2216f := x^(x^x); 2217 2218 2219 x 2220 x 2221f := x 2222 num_int(f,x= (1 .. 6),iterations=40); 2223 2224 2225*** ROUNDBF turned on to increase accuracy 2226 22271.10267709584e+36300 2228 2229 2230off rounded; 2231 2232 2233 2234end; 2235 2236Tested on x86_64-pc-windows CSL 2237Time (counter 1): 186 ms plus GC time: 32 ms 2238 2239End of Lisp run after 0.18+0.09 seconds 2240real 0.45 2241user 0.00 2242sys 0.06 2243