iconEuler Home

Exact Computation

In this notebook, we demonstrate how to use exact arithmetic in Euler.

Suppose we want to find the solution of

14 - Introduction to Exact Computation

Looking at the the plot, we see that there is a solution in [1,2].
>plot2d("x^x",0,2):

14 - Introduction to 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 optiona y value 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.55961046946  1.82545502292  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.55961046946  1.82545502292  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.55961046946  1.82545502292  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 precompile the derivative.
>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.
>integrate("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.00000000022 
     0.999999980607 
      1.00000041457 
      0.99999621873 
      1.00001808878 
     0.999950144682 
      1.00008198655 
     0.999920606001 
      1.00004175823 
     0.999990801457 
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 - Introduction to 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 - Introduction to 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 - Introduction to 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,niveau=0,r=10); ...
 plot2d(expr2,niveau=0,color=2,add=1):

14 - Introduction to 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 - Introduction to 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.64791433196  6.13171103227 ]
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.9833984375        3.6376953125 
      -2.6611328125        6.6162109375 
      -2.6611328125        6.6259765625 
      -2.6513671875        6.5869140625 
      -2.6513671875        6.5966796875 
      -2.6513671875        6.6064453125 
      -2.8662109375        8.0908203125 
       2.6513671875        6.1279296875 
       2.6513671875        6.1376953125 
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)
Column 1:
~-2.9800345111787907,-2.9800345111787876~ 
 ~-2.655428339950558,-2.6554283399505545~ 
~-2.8692456081316013,-2.8692456081315973~ 
  ~2.6479143319627383,2.6479143319627418~ 
Column 2:
  ~3.6352500667129695,3.6352500667129752~ 
  ~6.6048819930826532,6.6048819930826674~ 
  ~8.0887013125709597,8.0887013125709739~ 
  ~6.1317110322687496,6.1317110322687594~ 
Then add these points to the plot.
>h:=middle(r)'; plot2d(h[1],h[2],points=1,add=1,style="[]"):

14 - Introduction to Exact Computation

Euler Home