reduce-historical/r38/log/numeric.rlg

Tue Apr 15 00:34:16 2008 run on win32
on errcont;


bounds (x,x=(1 .. 2));


1 .. 2

bounds (2*x,x=(1 .. 2));


2 .. 4

bounds (x**3,x=(1 .. 2));


1 .. 8

bounds (x*y,x=(1 .. 2),y=(-1 .. 0));


 - 2 .. 0

bounds (x**3+y,x=(1 .. 2),y=(-1 .. 0));


0 .. 8

bounds (x**3/y,{x=(1 .. 2),y=(-1 .. -0.5)});


 - 16 .. -1

bounds (x**3/y,x=(1 .. 2),y=(-1 .. -0.5));


 - 16 .. -1

   % unbounded expression (pole at y=0)
bounds (x**3/y,x=(1 .. 2),y=(-1 .. 0.5));


***** unbounded in range 


on rounded;


bounds(e**x,x=(1 .. 2));


2.71828182846 .. 7.38905609893

bounds((1/2)**x,x=(1 .. 2));


0.25 .. 0.5

off rounded;



bounds(abs x,x=(1 .. 2));


1 .. 2

bounds(abs x,x=(-3 .. 2));


0 .. 3

bounds(abs x,x=(-3 .. -2));


2 .. 3


bounds(sin x,x=(1 .. 2));


 - 1 .. 1

 
on rounded;



bounds(sin x,x=(1 .. 2));


0.841470984808 .. 1

bounds(sin x,x=(1 .. 10));


 - 1 .. 1

bounds(sin x,x=(1001 .. 1002));


0.167266541974 .. 0.919990597586


bounds(log x,x=(1 .. 10));


0 .. 2.30258509299


bounds(tan x,x=(1 .. 1.1));


1.55740772465 .. 1.96475965725


bounds(cot x,x=(1 .. 1.1));


0.508968105239 .. 0.642092615934

bounds(asin x,x=(-0.6 .. 0.6));


 - 0.643501108793 .. 0.643501108793

bounds(acos x,x=(-0.6 .. 0.6));


0.927295218002 .. 2.21429743559


bounds(sqrt(x),x=(1 .. 1.1));


1 .. 1.04880884817

bounds(x**(7/3),x=(1 .. 1.1));


1 .. 1.2490589397

bounds(x**y,x=(1 .. 1.1),y=(2 .. 4));


1 .. 1.4641

 
off rounded;




% MINIMA  (steepest descent)

% Rosenbrock function (minimum extremely hard to find).
fktn := 100*(x1^2-x2)^2 + (1-x1)^2;


              4         2        2                2
fktn := 100*x1  - 200*x1 *x2 + x1  - 2*x1 + 100*x2  + 1

num_min(fktn, x1=-1.2, x2=1, accuracy=6);


{0.00000021870243927,{x1=0.999532844813,x2=0.999068072135}}


% infinitely many local minima
num_min(sin(x)+x/5, x=1);


{ - 1.33422674663,{x= - 1.77215430279}}


% bivariate polynomial
num_min(x^4 + 3 x^2 * y + 5 y^2 + x + y, x=0.1, y=0.2);


{ - 0.832523282274,{x= - 0.889601609042,y= - 0.33989805551}}



% ROOTS (non polynomial: damped Newton)

 num_solve (cos x -x, x=0,accuracy=6);


{x=0.739085133215}


   % automatically randomized starting point
num_solve (cos x -x,x, accuracy=6);


{x=0.739085133215}
  
 
   % syntactical errors: forms do not evaluate to purely 
   % numerical values
num_solve (cos x -x, x=a);


***** a invalid as number

num_solve (cos x -a, x=0);


***** error during function evaluation (e.g. singularity) 


num_solve (sin x = 0, x=3);


{x=3.14159265359}


  % blows up: no real solution exists
num_solve(sin x = 2, x=1);


***** Newton method does not converge 


  % solution in complex plane(only fond with complex starting point):
on complex;


*** Domain mode rounded changed to complex-rounded 

num_solve(sin x = 2, x=1+i);


{x=1.57079632542 + 1.31695789681*i}

off complex;


*** Domain mode complex-rounded changed to rounded 


  % blows up for derivative 0 in starting point 
num_solve(x^2-1, x=0);


***** division by zero 


  % succeeds because of perturbed starting point
num_solve(x^2-1, x=0.1);


{x=1.00000000033}


  % bivariate equation system
num_solve({sin x=cos y, x + y = 1},{x=1,y=2});


{x= - 52.1216769476,y=53.1216769476}

on rounded,evallhseqp;


sub(ws,{sin x=cos y, x + y = 1});


{ - 0.959549629985= - 0.959549629985,1=1}

off rounded,evallhseqp;


 
  % temporal member of the Barry Simon test sequence
sys :={sin (x) + y^2 + log(z) = 7,
       3*x + 2^y - z^3 = -90,
       x^2 + y^2 + z^(1/2) = 16};


                  2
sys := {sin(x) + y  + log(z)=7,

               y    3
        3*x + 2  - z =-90,

         2    2    1/2
        x  + y  + z   =16}

sol:=num_solve(sys,{x=1,y=1,z=1});


sol := {x=2.93087675819,y= - 2.29328251176,z=4.62601269017}

on rounded;


for each s in sys collect sub(sol,lhs s-rhs s);


{0,0,0}
  
off rounded;


clear sys,sol;


 
  % 2 examples taken from Nowak/Weimann (Tech.Rep TR91-10, ZIB Berlin)
 
  % #1: exp/sin combination

on rounded;


sys := {e**(x1**2 + x2**2)-3, x1 + x2 - sin(3(x1 + x2))};


           2     2
         x1  + x2
sys := {e          - 3,

         - sin(3*x1 + 3*x2) + x1 + x2}

num_solve(sys,x1=0.81, x2=0.82);


*** precision increased to 14 

*** precision increased to 18 

*** precision increased to 20 

{x1= - 0.256625076922,x2=1.01624596361}

sub(ws,sys);


{0,0}


  % 2nd example (semiconductor simulation), here computed with 
  % intermediate steps printed

alpha := 38.683;


alpha := 38.683

ni := 1.22e10;


ni := 1.22e+10

v := 100;


v := 100

d := 1e17;


d := 1.0e+17

sys := { e**(alpha*(x3-x1)) - e**(alpha*(x1-x2)) - d/ni,
         x2,
         x3,
         e**(alpha*(x6-x4)) - e**(alpha*(x4-x5)) + d/ni,
         x5 - v,
         x6 - v};


             77.366*x1                     38.683*x1 + 38.683*x2
sys := {( - e          - 8.19672131148e+6*e

             38.683*x2 + 38.683*x3   38.683*x1 + 38.683*x2
          + e                     )/e                     ,

        x2,

        x3,

             77.366*x4                     38.683*x4 + 38.683*x5
        ( - e          + 8.19672131148e+6*e

             38.683*x5 + 38.683*x6   38.683*x4 + 38.683*x5
          + e                     )/e                     ,

        x5 - 100,

        x6 - 100}

on trnumeric;


num_solve(sys,x1=1,x2=2,x3=3,x4=4,x5=5,x6=6,iterations=100);


*** computing symbolic Jacobian 

*** starting Newton iteration 

1. residue=(1.46329673989e+33 , 0 , 0 , 1.46329673988e+33 , 0.0 , 0.0)

, step length=163.473870223
 at ( - 1.97414885092 , 0 , 0 , 98.0258511491 , 100.0 , 100.0) 

2. residue=(5.38316786938e+32 , 0.0 , 0.0 , 5.38316786935e+32 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.94829770183 , 0.0 , 0.0 , 98.0517022982 , 100.0 , 100.0) 

3. residue=(1.98035678752e+32 , 0.0 , 0.0 , 1.98035678751e+32 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.92244655275 , 0.0 , 0.0 , 98.0775534473 , 100.0 , 100.0) 

4. residue=(7.28532548313e+31 , 0.0 , 0.0 , 7.28532548309e+31 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.89659540367 , 0.0 , 0.0 , 98.1034045963 , 100.0 , 100.0) 

5. residue=(2.68012146749e+31 , 0.0 , 0.0 , 2.68012146747e+31 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.87074425458 , 0.0 , 0.0 , 98.1292557454 , 100.0 , 100.0) 

6. residue=(9.8596158773e+30 , 0.0 , 0.0 , 9.85961587725e+30 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.8448931055 , 0.0 , 0.0 , 98.1551068945 , 100.0 , 100.0) 

7. residue=(3.62714997911e+30 , 0.0 , 0.0 , 3.62714997909e+30 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.81904195641 , 0.0 , 0.0 , 98.1809580436 , 100.0 , 100.0) 

8. residue=(1.33435390736e+30 , 0.0 , 0.0 , 1.33435390735e+30 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.79319080733 , 0.0 , 0.0 , 98.2068091927 , 100.0 , 100.0) 

9. residue=(4.90881369764e+29 , 0.0 , 0.0 , 4.90881369762e+29 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.76733965825 , 0.0 , 0.0 , 98.2326603418 , 100.0 , 100.0) 

10. residue=(1.8058516399e+29 , 0.0 , 0.0 , 1.80585163991e+29 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.74148850916 , 0.0 , 0.0 , 98.2585114908 , 100.0 , 100.0) 

11. residue=(6.64335692126e+28 , 0.0 , 0.0 , 6.64335692127e+28 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.71563736008 , 0.0 , 0.0 , 98.2843626399 , 100.0 , 100.0) 

12. residue=(2.4439544317e+28 , 0.0 , 0.0 , 2.4439544317e+28 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.689786211 , 0.0 , 0.0 , 98.310213789 , 100.0 , 100.0) 

13. residue=(8.99080590581e+27 , 0.0 , 0.0 , 8.99080590582e+27 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.66393506191 , 0.0 , 0.0 , 98.3360649381 , 100.0 , 100.0) 

14. residue=(3.30753265231e+27 , 0.0 , 0.0 , 3.30753265232e+27 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.63808391283 , 0.0 , 0.0 , 98.3619160872 , 100.0 , 100.0) 

15. residue=(1.21677326379e+27 , 0.0 , 0.0 , 1.21677326379e+27 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.61223276375 , 0.0 , 0.0 , 98.3877672363 , 100.0 , 100.0) 

16. residue=(4.47625868315e+26 , 0.0 , 0.0 , 4.47625868316e+26 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.58638161466 , 0.0 , 0.0 , 98.4136183853 , 100.0 , 100.0) 

17. residue=(1.64672354289e+26 , 0.0 , 0.0 , 1.6467235429e+26 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.56053046558 , 0.0 , 0.0 , 98.4394695344 , 100.0 , 100.0) 

18. residue=(6.05795736724e+25 , 0.0 , 0.0 , 6.05795736725e+25 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.5346793165 , 0.0 , 0.0 , 98.4653206835 , 100.0 , 100.0) 

19. residue=(2.2285979709e+25 , 0.0 , 0.0 , 2.2285979709e+25 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.50882816741 , 0.0 , 0.0 , 98.4911718326 , 100.0 , 100.0) 

20. residue=(8.19855376131e+24 , 0.0 , 0.0 , 8.19855376132e+24 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.48297701833 , 0.0 , 0.0 , 98.5170229817 , 100.0 , 100.0) 

21. residue=(3.01607937612e+24 , 0.0 , 0.0 , 3.01607937613e+24 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.45712586924 , 0.0 , 0.0 , 98.5428741308 , 100.0 , 100.0) 

22. residue=(1.10955359542e+24 , 0.0 , 0.0 , 1.10955359542e+24 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.43127472016 , 0.0 , 0.0 , 98.5687252798 , 100.0 , 100.0) 

23. residue=(4.08181956632e+23 , 0.0 , 0.0 , 4.08181956633e+23 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.40542357108 , 0.0 , 0.0 , 98.5945764289 , 100.0 , 100.0) 

24. residue=(1.50161750102e+23 , 0.0 , 0.0 , 1.50161750102e+23 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.37957242199 , 0.0 , 0.0 , 98.620427578 , 100.0 , 100.0) 

25. residue=(5.52414207128e+22 , 0.0 , 0.0 , 5.52414207133e+22 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.35372127291 , 0.0 , 0.0 , 98.6462787271 , 100.0 , 100.0) 

26. residue=(2.03221829814e+22 , 0.0 , 0.0 , 2.03221829815e+22 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.32787012383 , 0.0 , 0.0 , 98.6721298762 , 100.0 , 100.0) 

27. residue=(7.47611331856e+21 , 0.0 , 0.0 , 7.47611331863e+21 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.30201897474 , 0.0 , 0.0 , 98.6979810253 , 100.0 , 100.0) 

28. residue=(2.75030838977e+21 , 0.0 , 0.0 , 2.75030838979e+21 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.27616782566 , 0.0 , 0.0 , 98.7238321743 , 100.0 , 100.0) 

29. residue=(1.01178191348e+21 , 0.0 , 0.0 , 1.01178191349e+21 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.25031667658 , 0.0 , 0.0 , 98.7496833234 , 100.0 , 100.0) 

30. residue=(3.72213764917e+20 , 0.0 , 0.0 , 3.72213764921e+20 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.22446552749 , 0.0 , 0.0 , 98.7755344725 , 100.0 , 100.0) 

31. residue=(1.36929791834e+20 , 0.0 , 0.0 , 1.36929791835e+20 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.19861437841 , 0.0 , 0.0 , 98.8013856216 , 100.0 , 100.0) 

32. residue=(5.03736552996e+19 , 0.0 , 0.0 , 5.03736553001e+19 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.17276322933 , 0.0 , 0.0 , 98.8272367707 , 100.0 , 100.0) 

33. residue=(1.85314321614e+19 , 0.0 , 0.0 , 1.85314321616e+19 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.14691208024 , 0.0 , 0.0 , 98.8530879198 , 100.0 , 100.0) 

34. residue=(6.81733290764e+18 , 0.0 , 0.0 , 6.81733290771e+18 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.12106093116 , 0.0 , 0.0 , 98.8789390688 , 100.0 , 100.0) 

35. residue=(2.50795662034e+18 , 0.0 , 0.0 , 2.50795662037e+18 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.09520978207 , 0.0 , 0.0 , 98.9047902179 , 100.0 , 100.0) 

36. residue=(9.2262567997e+17 , 0.0 , 0.0 , 9.22625679991e+17 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.06935863299 , 0.0 , 0.0 , 98.930641367 , 100.0 , 100.0) 

37. residue=(3.39415019556e+17 , 0.0 , 0.0 , 3.39415019568e+17 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.04350748391 , 0.0 , 0.0 , 98.9564925161 , 100.0 , 100.0) 

38. residue=(1.24863807717e+17 , 0.0 , 0.0 , 1.24863807725e+17 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 1.01765633483 , 0.0 , 0.0 , 98.9823436652 , 100.0 , 100.0) 

39. residue=(4.59348278034e+16 , 0.0 , 0.0 , 4.59348278107e+16 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.991805185743 , 0.0 , 0.0 , 99.0081948143 , 100.0 , 100.0) 

40. residue=(1.68984787804e+16 , 0.0 , 0.0 , 1.68984787874e+16 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.965954036664 , 0.0 , 0.0 , 99.0340459633 , 100.0 , 100.0) 

41. residue=(6.21660292823e+15 , 0.0 , 0.0 , 6.21660293516e+15 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.940102887593 , 0.0 , 0.0 , 99.0598971124 , 100.0 , 100.0) 

42. residue=(2.28696040906e+15 , 0.0 , 0.0 , 2.28696041594e+15 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.914251738544 , 0.0 , 0.0 , 99.0857482616 , 100.0 , 100.0) 

43. residue=(8.41325715099e+14 , 0.0 , 0.0 , 8.41325721962e+14 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.888400589553 , 0.0 , 0.0 , 99.1115994107 , 100.0 , 100.0) 

44. residue=(3.09506431748e+14 , 0.0 , 0.0 , 3.09506438604e+14 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.862549440721 , 0.0 , 0.0 , 99.1374505601 , 100.0 , 100.0) 

45. residue=(1.13861050984e+14 , 0.0 , 0.0 , 1.13861057839e+14 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.836698292322 , 0.0 , 0.0 , 99.1633017098 , 100.0 , 100.0) 

46. residue=(4.18871376415e+13 , 0.0 , 0.0 , 4.18871444947e+13 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.8108471451 , 0.0 , 0.0 , 99.1891528608 , 100.0 , 100.0) 

47. residue=(1.54094146219e+13 , 0.0 , 0.0 , 1.54094214749e+13 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.784996001075 , 0.0 , 0.0 , 99.2150040149 , 100.0 , 100.0) 

48. residue=(5.66880467397e+12 , 0.0 , 0.0 , 5.6688115269e+12 , 0.0 , 0.0)

, step length=0.0365590456369
 at ( - 0.759144865742 , 0.0 , 0.0 , 99.2408551778 , 100.0 , 100.0) 

49. residue=(2.08543452966e+12 , 0.0 , 0.0 , 2.08544138253e+12 , 0.0 , 0.0)

, step length=0.036559045637
 at ( - 0.733293754037 , 0.0 , 0.0 , 99.2667063642 , 100.0 , 100.0) 

50. residue=(767186323467.0 , 0.0 , 0.0 , 767193176319.0 , 0.0 , 0.0)

, step length=0.0365590456375
 at ( - 0.70744270656 , 0.0 , 0.0 , 99.2925576149 , 100.0 , 100.0) 

51. residue=(282229910057.0 , 0.0 , 0.0 , 282236762901.0 , 0.0 , 0.0)

, step length=0.0365590456414
 at ( - 0.681591833671 , 0.0 , 0.0 , 99.3184090402 , 100.0 , 100.0) 

52. residue=(103824415727.0 , 0.0 , 0.0 , 103831268569.0 , 0.0 , 0.0)

, step length=0.0365590456703
 at ( - 0.655741435353 , 0.0 , 0.0 , 99.3442609401 , 100.0 , 100.0) 

53. residue=(3.81927022457e+10 , 0.0 , 0.0 , 3.8199555087e+10 , 0.0 , 0.0)

, step length=0.0365590458835
 at ( - 0.629892327003 , 0.0 , 0.0 , 99.3701141301 , 100.0 , 100.0) 

54. residue=(1.40481443717e+10 , 0.0 , 0.0 , 1.40549972128e+10 , 0.0 , 0.0)

, step length=0.0365590474585
 at ( - 0.604046724769 , 0.0 , 0.0 , 99.3959708274 , 100.0 , 100.0) 

55. residue=(5.16585846951e+9 , 0.0 , 0.0 , 5.17271131074e+9 , 0.0 , 0.0)

, step length=0.0365590590969
 at ( - 0.578210650353 , 0.0 , 0.0 , 99.4218370614 , 100.0 , 100.0) 

56. residue=(1.89824960589e+9 , 0.0 , 0.0 , 1.90510244878e+9 , 0.0 , 0.0)

, step length=0.0365591450932
 at ( - 0.552400454575 , 0.0 , 0.0 , 99.4477292394 , 100.0 , 100.0) 

57. residue=(6.96167585052e+8 , 0.0 , 0.0 , 7.03020440594e+8 , 0.0 , 0.0)

, step length=0.0365597805245
 at ( - 0.526660451901 , 0.0 , 0.0 , 99.4736920939 , 100.0 , 100.0) 

58. residue=(2.53957444815e+8 , 0.0 , 0.0 , 2.60810394004e+8 , 0.0 , 0.0)

, step length=0.0365644757288
 at ( - 0.501110133883 , 0.0 , 0.0 , 99.4998482048 , 100.0 , 100.0) 

59. residue=(9.13074482936e+7 , 0.0 , 0.0 , 9.81610893493e+7 , 0.0 , 0.0)

, step length=0.036599167075
 at ( - 0.476067267452 , 0.0 , 0.0 , 99.526538163 , 100.0 , 100.0) 

60. residue=(3.15519019482e+7 , 0.0 , 0.0 , 3.8410646839e+7 , 0.0 , 0.0)

, step length=0.0368554086395
 at ( - 0.452345623749 , 0.0 , 0.0 , 99.5547446298 , 100.0 , 100.0) 

61. residue=(9.77481469301e+6 , 0.0 , 0.0 , 1.66708124709e+7 , 0.0 , 0.0)

, step length=0.0387445972031
 at ( - 0.431825342687 , 0.0 , 0.0 , 99.5876089247 , 100.0 , 100.0) 

62. residue=(2.23533931681e+6 , 0.0 , 0.0 , 9.38172386018e+6 , 0.0 , 0.0)

, step length=0.0527640781991
 at ( - 0.417764764235 , 0.0 , 0.0 , 99.6384650755 , 100.0 , 100.0) 

63. residue=(2.23262484734e+5 , 0.0 , 0.0 , 8.19715321429e+6 , 0.0 , 0.0)

, step length=0.204739774363
 at ( - 0.41222548566 , 0.0 , 0.0 , 99.843129903 , 100.0 , 100.0) 

64. residue=(2.23088062724e+5 , 0.0 , 0.0 , 8.19035290355e+6 , 0.0 , 0.0)

, step length=0.383303021247
 at ( - 0.412224950141 , 0.0 , 0.0 , 100.226432924 , 100.0 , 100.0) 

65. residue=(2.22216670074e+5 , 0.0 , 0.0 , 7.2288078e+6 , 0.0 , 0.0)

, step length=0.129870827173
 at ( - 0.412222274586 , 0.0 , 0.0 , 100.356303751 , 100.0 , 100.0) 

66. residue=(1.66845393097e+5 , 0.0 , 0.0 , 1.93472870727e+6 , 0.0 , 0.0)

, step length=0.048267265871
 at ( - 0.412051690235 , 0.0 , 0.0 , 100.404570716 , 100.0 , 100.0) 

67. residue=(1653.19388834 , 0.0 , 0.0 ,  - 3.3219398641e+5 , 0.0 , 0.0)

, step length=0.00800369958771
 at ( - 0.411535983805 , 0.0 , 0.0 , 100.412557784 , 100.0 , 100.0) 

68. residue=(0.166671264823 , 0.0 , 0.0 ,  - 6386.15622619 , 0.0 , 0.0)

, step length=0.00100689373229
 at ( - 0.411530770947 , 0.0 , 0.0 , 100.411550903 , 100.0 , 100.0) 

69. residue=( - 0.000000122 , 0.0 , 0.0 ,  - 2.48519530414 , 0.0 , 0.0)

, step length=0.0000201252363669
 at ( - 0.411530770421 , 0.0 , 0.0 , 100.411530778 , 100.0 , 100.0) 

{x1= - 0.411530770421,x2=0.0,x3=0.0,x4=100.41153077,x5=100.0,x6=100.0}

off trnumeric;


clear alpha,ni,v,d,sys;


off rounded;



% INTEGRALS
 
num_int( x**2,x=(1 .. 2),accuracy=3);


 7
---
 3


  % 1st case: using formal integral
needle := 1/(10**-4 + x**2);


              10000
needle := --------------
                  2
           10000*x  + 1

num_int(needle,x=(-1 .. 1),accuracy=3);


312.159332022
           % 312.16

  % no formal integral, but easy Chebyshev fit
num_int(sin x/x,x=(1 .. 10));


0.712264523852


  % using a Chebyshev fit of order 60
num_int(exp(-x**2),x=(-10 .. 10),accuracy=3);


1.77245387654
     % 1.772

   % cases with singularities

num_int(1/sqrt x ,x=(0 .. 1),accuracy=2);


2
          % 1.999

num_int(1/sqrt abs x ,x=(-1 .. 1),iterations=50);


3.99999231465
     % 3.999

   % simple multidimensional integrals
num_int(x+y,x=(0 .. 1),y=(2 .. 3));


3.0


num_int(sin(x+y),x=(0 .. 1),y=(0 .. 1));


0.773135425204


% some integrals with infinite bounds
 
on rounded;

 % for the error function

num_int(e^(-x) ,x=(0 .. infinity));


1.00000034605
  % 1.000

2/sqrt(pi)* num_int(e^(-x^2) ,x=(0 .. infinity));


1.00000003784
 % 1.00

2/sqrt(pi)* num_int(e^(-x^2), x=(-infinity .. infinity));


2.00000007569
 % 2.00

num_int(sin(x) * e^(-x), x=(0 .. infinity));


0.500000522701
 % 0.500
 
off rounded;


 
% APPROXIMATION
 
  %approximate sin x by a cubic polynomial 
num_fit(sin x,{1,x,x**2,x**3},x=for i:=0:20 collect 0.1*i);


                     3                   2
{ - 0.0847539694989*x  - 0.134641944765*x  + 1.06263064633*x - 0.00519313406463,

 { - 0.00519313406463,1.06263064633, - 0.134641944765, - 0.0847539694989}}

 
  % approximate x**2 by a harmonic series in the interval [0,1]
num_fit(x**2,1 . for i:=1:5 join {sin(i*x)},
               x=for i:=0:10 collect i/10);


{ - 1.3095780871*sin(5*x) + 7.1637556683*sin(4*x) - 18.549018248*sin(3*x)

  + 26.5601709095*sin(2*x) - 19.4492185507*sin(x) - 0.00197199704297,

 { - 0.00197199704297, - 19.4492185507,26.5601709095, - 18.549018248,

  7.1637556683, - 1.3095780871}}

 
  % approximate a set of points by a polynomial
pts:=for i:=1 step 0.1 until 3 collect i$


vals:=for each p in pts collect (p+2)**3$


num_fit(vals,{1,x,x**2,x**3},x=pts);


      3                  2
{1.0*x  + 5.99999999998*x  + 12.0*x + 7.99999999998,{7.99999999998,12.0,

    5.99999999998,1.0}}

  % compute the approximation error
on rounded;


first ws - (x+2)**3;


                   3                             2
2.50954812486e-12*x  - 0.0000000000155884194442*x  + 0.0000000000306474845502*x

 - 0.0000000000188205007134

off rounded;




 
% ODE SOLUTION (Runge-Kutta)
 
depend(y,x);



   % approximate y=y(x) with df(y,x)=2y in interval [0 : 5]
num_odesolve(df(y,x)=y,y=2,x=(0 .. 5),iterations=20);


{{x,y},

 {0.0,2.0},

 {0.25,2.56805083337},

 {0.5,3.2974425414},

 {0.75,4.23400003322},

 {1.0,5.43656365691},

 {1.25,6.98068591491},

 {1.5,8.96337814065},

 {1.75,11.509205352},

 {2.0,14.7781121978},

 {2.25,18.9754716726},

 {2.5,24.3649879213},

 {2.75,31.2852637682},

 {3.0,40.1710738461},

 {3.25,51.5806798341},

 {3.5,66.2309039169},

 {3.75,85.0421639995},

 {4.0,109.196300065},

 {4.25,140.210824692},

 {4.5,180.034262599},

 {4.75,231.168569052},

 {5.0,296.826318202}}

 
   % same with negative direction
num_odesolve(df(y,x)=y,y=2,x=(0 .. -5),iterations=20);


{{x,y},

 {0.0,2.0},

 {-0.25,1.55760156614},

 {-0.5,1.21306131943},

 {-0.75,0.944733105483},

 {-1.0,0.735758882344},

 {-1.25,0.573009593722},

 {-1.5,0.446260320298},

 {-1.75,0.347547886902},

 {-2.0,0.270670566474},

 {-2.25,0.210798449125},

 {-2.5,0.164169997249},

 {-2.75,0.127855722414},

 {-3.0,0.0995741367363},

 {-3.25,0.0775484156639},

 {-3.5,0.060394766845},

 {-3.75,0.0470354917124},

 {-4.0,0.0366312777778},

 {-4.25,0.0285284678182},

 {-4.5,0.0222179930767},

 {-4.75,0.0173033904064},

 {-5.0,0.0134758939983}}


   % giving a nice picture when plotted
num_odesolve(df(y,x)=1- x*y**2 ,y=0,x=(0 .. 4),iterations=20);


{{x,y},

 {0.0,0.0},

 {0.2,0.199600912188},

 {0.4,0.393714914166},

 {0.6,0.569482634406},

 {0.8,0.710657687564},

 {1.0,0.805480022354},

 {1.2,0.852604291055},

 {1.4,0.860563377356},

 {1.6,0.842333334456},

 {1.8,0.80999200878},

 {2.0,0.772211952811},

 {2.2,0.734163640068},

 {2.4,0.698433235122},

 {2.6,0.666019196492},

 {2.8,0.637070046905},

 {3.0,0.611341375657},

 {3.2,0.588447372601},

 {3.4,0.567985133759},

 {3.6,0.549587947292},

 {3.8,0.532942255624},

 {4.0,0.517787833735}}


   % system of ordinary differential equations
depend(y,x);


depend(z,x);


num_odesolve(
    {df(z,x) = y, df(y,x)=  y+x},
    {z=2, y=4},
     x=(0 .. 5),iterations=20);


{{x,z,y},

 {0.0,2.0,4.0},

 {0.25,3.13887708344,5.17012708344},

 {0.5,4.61860635349,6.74360635349},

 {0.75,6.55375008305,8.83500008305},

 {1.0,9.09140914227,11.5914091423},

 {1.25,12.4204647873,15.2017147873},

 {1.5,16.7834453516,19.9084453516},

 {1.75,22.4917633799,26.0230133799},

 {2.0,29.9452804945,33.9452804945},

 {2.25,39.6574291816,44.1886791816},

 {2.5,52.2874698032,57.4124698032},

 {2.75,68.6819094205,74.4631594205},

 {3.0,89.9276846154,96.4276846154},

 {3.25,117.420449585,124.701699585},

 {3.5,152.952259792,161.077259792},

 {3.75,198.824159999,207.855409999},

 {4.0,257.990750164,267.990750164},

 {4.25,334.245811731,345.277061731},

 {4.5,432.460656499,444.585656499},

 {4.75,558.890172631,572.171422631},

 {5.0,721.565795506,736.065795506}}



%----------------- Chebyshev fit -------------------------

on rounded;



func := x**2 * (x**2 - 2) * sin x;


                2   2
func := sin(x)*x *(x  - 2)

ord := 15;


ord := 15


cx:=chebyshev_fit(func,x=(0 .. 2),ord)$


cp:=first cx;


                            13                       12                      11
cp := 0.000000620105096185*x   + 0.0000168737305191*x   - 0.000269014332288*x

                      10                     9                      8
 + 0.000155646029006*x   + 0.00848163760265*x  + 0.000272748604876*x

                   7                     6                  5
 - 0.183540091904*x  + 0.00010680840443*x  + 1.33329694616*x

                        4                  3                          2
 + 0.00000770692780683*x  - 2.00000091554*x  + 0.0000000501515695639*x

 - 0.000000000784989184766*x - 4.86055640181e-13

cc:=second cx;


cc := {2.69320512829,

       2.76751928466,

       2.25642507569,

       0.955452569949,

       0.0509075944268,

        - 0.0868248678183,

        - 0.0170919216091,

       0.00104527137626,

       0.000349190502034,

        - 0.00000253521592323,

        - 0.00000280798840641,

        - 0.0000000157676044858,

       0.0000000121753402195,

       0.000000000118269801846,

        - 0.0000000000331230439026}


for u:=0 step 0.2 until 2 do write
    "x:",u," true value:",sub(x=u,func),
           " Chebyshev eval:", chebyshev_eval(cc,x=(0 .. 2),x=u),
           " Chebyshev polynomial:",sub(x=u,cp);


x:0 true value:0 Chebyshev eval: - 4.85167461761e-13 Chebyshev polynomial:

 - 4.86055640181e-13

x:0.2 true value: - 0.0155756755343 Chebyshev eval: - 0.0155756755339

 Chebyshev polynomial: - 0.015575675548

x:0.4 true value: - 0.114644759976 Chebyshev eval: - 0.114644759976

 Chebyshev polynomial: - 0.114644759974

x:0.6 true value: - 0.333364916292 Chebyshev eval: - 0.333364916292

 Chebyshev polynomial: - 0.333364916295

x:0.8 true value: - 0.624386741519 Chebyshev eval: - 0.624386741519

 Chebyshev polynomial: - 0.624386741504

x:1 true value: - 0.841470984808 Chebyshev eval: - 0.841470984808

 Chebyshev polynomial: - 0.841470984841

x:1.2 true value: - 0.751596318924 Chebyshev eval: - 0.751596318924

 Chebyshev polynomial: - 0.751596318876

x:1.4 true value: - 0.0772592588311 Chebyshev eval: - 0.0772592588311

 Chebyshev polynomial: - 0.0772592588864

x:1.6 true value:1.43298871732 Chebyshev eval:1.43298871732

 Chebyshev polynomial:1.43298871738

x:1.8 true value:3.91253024182 Chebyshev eval:3.91253024182

 Chebyshev polynomial:3.91253024177

x:2.0 true value:7.27437941461 Chebyshev eval:7.27437941461

 Chebyshev polynomial:7.27437941467


% integral

   % integrate coefficients
ci := chebyshev_int(cc,x=(0 .. 2));


ci := {0.0310113015322,

       0.2183900263,

       0.453016678678,

       0.367586246877,

       0.130284679721,

       0.00679995160359,

        - 0.00732251159954,

        - 0.00124579372222,

       0.0000654879120115,

       0.0000195554716911,

        - 0.000000125972415937,

        - 0.000000128189261211,

        - 0.000000000661911428653,

       0.000000000469556279362,

       4.22392149449e-12}


   % compare with true values (normalized absolute term)
ci0:=chebyshev_eval(ci,x=(0 .. 2),x=0)$


ifunc := int(func,x)$

 if0 := sub(x=0,ifunc);


if0 :=  - 28.0


for u:=0 step 0.2 until 2 do write
       {u,sub(x=u,ifunc) - if0,
          chebyshev_eval(ci,x=(0 .. 2),x=u) - ci0};


{0,0,0}

{0.2, - 0.000785836355117, - 0.00078583635293}

{0.4, - 0.0119047051867, - 0.0119047051858}

{0.6, - 0.0548116700418, - 0.0548116700408}

{0.8, - 0.150297976106, - 0.150297976105}

{1, - 0.299838223412, - 0.29983822341}

{1.2, - 0.466528961073, - 0.466528961072}

{1.4, - 0.561460555384, - 0.561460555383}

{1.6, - 0.441445769516, - 0.441445769514}

{1.8,0.0768452822437,0.0768452822437}

{2.0,1.18309971762,1.18309971762}


% derivative

   % differentiate coefficients
cd := chebyshev_df(cc,x=(0 .. 2))$


   % compute coefficients of derivative
cds := second chebyshev_fit(df(func,x),x=(0 .. 2),ord)$


   % compare coefficients
for i:=1:ord do write {part(cd,i),part(cds,i)};


{10.4140931324,10.4140931324}

{9.23338917839,9.2333891784}

{4.87905456308,4.87905456307}

{0.207688875651,0.207688875654}

{ - 0.853660856614, - 0.853660856625}

{ - 0.199571879764, - 0.19957187976}

{0.0145878215687,0.0145878215579}

{0.00553117954514,0.00553117954883}

{ - 0.000045977698902, - 0.0000459777097756}

{ - 0.0000558684874082, - 0.0000558684837245}

{ - 0.00000034381228384, - 0.00000034382315144}

{0.000000291280720039,0.00000029128440905}

{0.00000000307501484799,0.00000000306414359587}

{ - 0.000000000927445229273, - 0.000000000923750392908}

{0, - 0.0000000000109242472928}


clear func,ord,cc,cx,cd,cds,ci,ci0;



% One from ISSAC '97 -- should be  ~ 1.10*10^36300

f := x^(x^x);


       x
      x
f := x
 num_int(f,x= (1 .. 6),iterations=40);


*** ROUNDBF turned on to increase accuracy 

1.10267709584e+36300


off rounded;



end;


Time for test: 1808 ms, plus GC time: 73 ms