Tue Feb 10 12:28:01 2004 run on Linux
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.3095781041*sin(5*x) + 7.16375576968*sin(4*x) - 18.5490185227*sin(3*x)
+ 26.5601713287*sin(2*x) - 19.4492188858*sin(x) - 0.00197199701919,
{ - 0.00197199701919, - 19.4492188858,26.5601713287, - 18.5490185227,
7.16375576968, - 1.3095781041}}
% 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.000000620104755125*x + 0.000016873735186*x - 0.000269014360576*x
10 9 8
+ 0.000155646128984*x + 0.00848163737447*x + 0.00027274895671*x
7 6 5
- 0.183540092277*x + 0.000106808674377*x + 1.33329694603*x
4 3 2
+ 0.00000770696844232*x - 2.00000091554*x + 0.0000000501522801066*x
- 0.000000000785017162386*x - 4.85611550971e-13
cc:=second cx;
cc := {2.69320512829,
2.76751928466,
2.25642507569,
0.955452569949,
0.0509075944268,
- 0.0868248678183,
- 0.0170919216091,
0.00104527137626,
0.000349190502034,
- 0.00000253521592324,
- 0.00000280798840646,
- 0.0000000157676044831,
0.0000000121753403333,
0.000000000118269806472,
- 0.0000000000331229560099}
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.85611550971e-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.000000125972415938,
- 0.000000128189261218,
- 0.00000000066191142873,
0.000000000469556280358,
4.2239216597e-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.000055868487404, - 0.0000558684837246}
{ - 0.000000343812283659, - 0.000000343823151536}
{0.000000291280725231,0.000000291284409096}
{0.00000000307501496826,0.00000000306414403996}
{ - 0.000000000927442768278, - 0.000000000923750378105}
{0, - 0.0000000000109240104452}
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: 1570 ms, plus GC time: 30 ms