iconEuler Home

Exact Computation

In this notebook, we demonstrate how to use exact arithmetic in Euler. The core functions are available on the following page.

Exact Arithmetic

Suppose we want to find the solution of

14 - Exact Computation

Looking at the the plot, we see that there is a solution in [1,2].

>plot2d("x^x",0,2):

14 - Exact Computation

Thus we can use the bisection method to find it.

>bisect("x^x-2",1,2)
1.55961046946

If we need the zero in an Euler program, we want a fast solution. The Newton method is very good here. We need the derivative of the function which we compute with Maxima automatically.

>mxmnewton("x^x-2",2)
1.55961046946

The secant method is implemented in solve. It is as good as the Newton method. The solver functions accept an optional value y to solve for.

>solve("x^x",2,y=2)
1.55961046946

We can implement an inversion operator for any function or expression using solve.

>function map invf(y;f) := solve(f,y,y=y)

Here is an example, where we compute the inverse of y=x^x for values 1,2,3,4.

>invf(1:4,"x^x")
[1,  1.55961,  1.82546,  2]

Maybe we want the fastest possible solver for this, the Newton solver. Then we need a derivative.

Since we need the derivative more than once, we will precompile it calling Maxima at compile time, or passing the derivative to our function as a parameter.

The following demonstrates the parameter method. We implement a function for the inverse of f(x)=x^x.

>function f(x) &= x^x;

We call the Newton method using the function parameters. For simplicity, we always start at x=y.

>function map invf(y;f,fder) := newton(f,fder,y,y=y);

Due to the mapping, this works for vectors too.

>invf(1:4,&f(x),&diff(f(x),x))
[1,  1.55961,  1.82546,  2]

Here is the pre-computation method.

>function map invxx(s) := newton("&:f(x)","&:diff(f(x),x)",1,y=s)

This is a fast solution. The derivative is built into the function.

>type invxx
function map invxx (s)
useglobal; return newton("x^x","x^x*(log(x)+1)",1,y=s) 
endfunction
>invxx(1:4)
[1,  1.55961,  1.82546,  2]

Interval Newton Method

This method works, but we know nothing about its accuracy.

Euler can use the interval Newton method to find a guaranteed inclusion. The interval Newton method uses the Brower fixed point theorem to guarantee a proper solution.

>mxminewton("x^x-2",~1,2~)
~1.559610469462368,1.559610469462371~

Again, it might be wise to pre-compile the derivative and add it as a parameter to the function invfi().

>function map invfi (y,f,fder) := inewton(f,fder,~1,2~,y=y)
>invfi(1:4,&f(x),&diff(f(x),x))
[ ~1,1.0000000000000004~,  ~1.559610469462368,1.559610469462371~,
~1.825455022924829,1.825455022924831~,  ~1.9999999999999993,2~ ]

Exact Integration

Interval inclusions is even more valuable for integration. Let us compute an inclusion of the integral from 1 to 2.

>mxmiint("x^x",1,2)
~2.05044623453469,2.05044623453477~

This inclusion is very narrow. Maxima has to compute derivatives of degree 10 to find it. Then a Taylor series expansion is used in small subintervals to get the inclusion.

The interval simpson method is much faster. However, it can cope with the accuracy, if the number of subintervals is high enough.

>mxmisimpson("x^x",1,2,1000)
~2.0504462345342,2.0504462345353~

The Romberg method is usually very reliable too and much faster.

>romberg("x^x",1,2)
2.05044623453

In the interval [0,1] both method will fail, and we have to resort to the stable and ultra-fast Gauss integration, which does not give any error estimate, however. We split into 10 intervals for more accuracy.

>gauss("x^x",0,1,10)
0.783430300312

It is difficult to get an estimate for the integral of this function in [0,1], since the function is not differentiable in 0. However, we can try splitting up the integral in two regions.

>I1:=mxmiint("x^x",0.1,1,deg=7)
~0.696335037,0.696335102~

The function is monoton decreasing in [0,0.1]. Thus we get an estimate using a simple upper and lower sum. The error is between 0 and f(a)-f(b) for such functions.

Adding the integral over the other part yields an exact inclusion.

>h:=0.1/100000; x:=0:h:0.1-h; y:=sum(x^x)*h; I1+~y-h*(1-0.1^0.1),y~
~0.78343039,0.78343067~

Even for a function like sinc(x) defined as sin(x)/x, which uses a special case for x~=0, we can use the Romberg method.

>romberg("sinc",0,pi)
1.85193705198

The Gauß method used in integrate with 10 sub-intervals gets very good results too. Note, that sinc is an entire function.

>gauss("sinc",0,pi)
1.85193705198

An inclusion is not easy, since mxmiint fails for values close to 0. So we use the Taylor series.

First we integrate from 1 to pi exactly.

>I1:=mxmiint("sin(x)/x",1,pi)
~0.9058539816132,0.9058539816174~

Then we add an estimate for the integral from 0 to 1. This estimate comes from integrating the power series and from estimating the remainder term of the power series.

>I1+mxmeval("integrate(taylor(sin(x),x,0,20)/x,x,0,1)")±2*pi^21/21!
~1.85193705,1.851937054~

Inclusions for Differential Equations

There is also a high accuracy solver for differential equations. The derivatives are again computed by Maxima. We try the example y'=2xy, y(0)=1, with the solution exp(-x^2).

>x:=0:0.05:1; y:=mxmidgl("2*x*y",x,1); y[-1]
~2.71828182845901,2.71828182845908~

Inclusions for Linear Systems.

For linear systems with badly conditioned matrices, we often get very inaccurate results.

>A:=hilb(10); b:=sum(A); A\b
            1 
            1 
            1 
     0.999996 
      1.00002 
      0.99995 
      1.00008 
     0.999921 
      1.00004 
     0.999991 

Euler has an interval solver, which uses iterative process and yields a nice inclusion.

>ilgs(A,b)
 ~0.99999999999999978,1.0000000000000002~ 
 ~0.99999999999999956,1.0000000000000004~ 
 ~0.99999999999999822,1.0000000000000018~ 
 ~0.99999999999999334,1.0000000000000067~ 
  ~0.99999999999995237,1.000000000000048~ 
   ~0.9999999999996404,1.000000000000359~ 
     ~0.999999999999152,1.00000000000085~ 
   ~0.9999999999995126,1.000000000000486~ 
    ~0.9999999999997816,1.00000000000022~ 
 ~0.99999999999996947,1.0000000000000302~ 

If only a good solution is needed, we can use the residue iteration with xlgs.

>xlgs(A,b)
            1 
            1 
            1 
            1 
            1 
            1 
            1 
            1 
            1 
            1 

To get such improvements, Euler uses an exact scalar product. The first of the following commands delivers a wrong result, since the 1 is swallowed by the large value 1e20. The second command computes an exact result. "residuum(A,b,r)" simply computes the matrix A.b-r, but in an exact way.

>1e20+1-1e20, residuum([1e20,1,-1e20],[1,1,1]',0),
0
1

An Example

The following example (Rump et al.) is a polynomial, which evaluates very badly.

>p:=[-945804881,1753426039,-1083557822,223200658]; ...
 t:=linspace(1.61801916,1.61801917,100); ...
 plot2d(t-1.61801916,evalpoly(t,p)):

14 - Exact Computation

As we see, the plot is completely wrong. Euler can compute a very accurate evaluation using a residuum iteration. This works, since the evaluation of a polynomial can be rewritten as a linear system.

>plot2d(t-1.61801916,xevalpoly(t,p,eps=1e-17)):

14 - Exact Computation

Inclusion for the Solution of a System of Equations

For two dimensional problems, Euler has some functions and methods to find guaranteed inclusions.

Let us start with the following example. We wish to solve

14 - Exact Computation

>expr1 &= x*sin(y)-y*cos(x)-5; expr2 &= x^2+cos(y)-8;

First of all, we can plot an image of the zero sets of both functions.

>plot2d(expr1,level=0,a=-10,b=10,c=-10,d=10,n=80); ...
 plot2d(expr2,level=0,color=2,add=1):

14 - Exact Computation

In this example, we can solve one equation first.

>:: solve(expr2,x)
            [x = - sqrt(8 - cos(y)), x = sqrt(8 - cos(y))]

Then insert one of the solutions into the other.

>&at(expr1,solve(expr2,x)[2])
       - y cos(sqrt(8 - cos(y))) + sqrt(8 - cos(y)) sin(y) - 5

Now solve for y.

>sy:=bisect(&subst(x,y,at(expr1,solve(expr2,x)[2])),6,8)
6.13171103227

Now we can also find an inclusion for this using the interval Newton method.

>syi:=mxminewton(&subst(x,y,at(expr1,solve(expr2,x)[2])),~6.1,6.2~)
~6.131711032268751,6.131711032268759~

Substitute this y back into the formula for x.

>sxi:=&"rhs(solve(expr2,x)[2])"(y=syi)
~2.647914331962739,2.647914331962741~

We add this point to our plot.

>plot2d(middle(sxi),middle(syi),points=1,style="[]",add=1):

14 - Exact Computation

Euler has also a higher dimensional Newton method. There is a service function for two dimensional problems in x and y. It uses Maxima to compute the partial derivatives.

>mxmnewtonfxy(expr1,expr2,[2.5,6])
[2.64791,  6.13171]

The higher dimensional interval Newton method uses functions. So we define the functions for it. First the function f.

>function f([x,y]) &= [expr1,expr2]
                                                  2
              [x sin(y) - cos(x) y - 5, cos(y) + x  - 8]

Likewise a function for the Jocobian matrix of f.

>function df([x,y]) &= jacobian(f(x,y),[x,y])
               [ sin(y) + sin(x) y  x cos(y) - cos(x) ]
               [                                      ]
               [        2 x             - sin(y)      ]

The multi-dimensional interval Newton method returns an interval inclusion plus a flag. The flag is 1, if the inclusion is a guaranteed inclusion.

>{res,valid}:=inewton2("f","df",[~2.6,2.7~,~6.1,6.2~]); res,
[ ~2.647914331962738,2.647914331962742~,
~6.131711032268749,6.13171103226876~ ]
>valid
1

For two dimensions, there is a general interval method, splitting the area into small rectangles. It returns a list of possible intervals. Everything outside this list is guaranteed not to contain a zero.

We print only the midpoints of the intervals.

>res:=mxmibisectfxy(expr1,expr2,~-10,10~,~-10,10~,0.01); middle(res)
      -2.9834        3.6377 
     -2.66113       6.61621 
     -2.66113       6.62598 
     -2.65137       6.58691 
     -2.65137       6.59668 
     -2.65137       6.60645 
     -2.86621       8.09082 
      2.65137       6.12793 
      2.65137        6.1377 

We can use any of this as a starting point for the interval Newton method and keep only those, which are guaranteed inclusions (maybe loosing some).

>function test(res) ...
 h=[];
 loop 1 to rows(res)
 {y,valid}=inewton2("f","df",res[#]);
 if valid then h=h_y; endif;
 end
 return h
 endfunction
>r:=test(res)
Interval 4 x 2 matrix

~-2.9800345111787907,-2.9800345111787876~     ...
 ~-2.655428339950558,-2.6554283399505545~     ...
~-2.8692456081316013,-2.8692456081315973~     ...
  ~2.6479143319627383,2.6479143319627418~     ...

Then add these points to the plot.

>h:=middle(r)'; plot2d(h[1],h[2],points=1,add=1,style="[]"):

14 - Exact Computation

Euler Home