iconEuler Examples

Convergence of the Newton Method

by R. Grothmann

Let us test the convergence of the well known iteration to sqrt(2).
The iterate function applies the function or expression until the
limit is reached numerically.
>longestformat; iterate("(x+2/x)/2",2)
      1.414213562373095 
An additional parameter sets the number of iterations.
>iterate("(x+2/x)/2",2,6)
Column 1 to 2:
                    1.5       1.416666666666667 
Column 3 to 4:
       1.41421568627451        1.41421356237469 
Column 5 to 6:
      1.414213562373095       1.414213562373095 
We can try to prove the existence of a limit with interval arithmetic.
First, we define a small interval h around 0.
>h=~-1e-13,1e-13~
                         ~-1e-013,1e-013~ 
For easy of use, we define the interation operator as a Maxima
expression.
>function f(x) &= (x+2/x)/2
                                    2
                                x + -
                                    x
                                -----
                                  2

Now we use the mean value theorem to compute the image of a
small interval around sqrt(2). We use

f(x+h) = f(x) + f'(xi)*h

and evaluate the right hand side in interval fashion.
>x0=sqrt(2); res=mxmeval("f(x)",x=x0)+mxmeval("diff(f(x),x)",x=x0+h)*h
  ~1.4142135623730947,1.4142135623730951~ 
New the fixed point theorem tells us, that there must be a fixed
point in the interval f(x+h), if f(x+h) is contained in x+h.
>res << x0+h
                      1 
There is a shortcut for this kind of interval evaluation (using
the mean value theorem, resp. the Taylor series of degree 1).
>mxmieval("(x+2/x)/2",x0+h)
    ~1.414213562373094,1.414213562373096~ 
Now continue with the Newton example.

The good convergence is no surprise, since the iteration is
in fact the Newton iteration.
>&expand(x - (x^2-2)/diff(x^2-2,x))
                                x   1
                                - + -
                                2   x

There is a function to start the Newton iteration automatically.
>mxminewton("x^2-2",~1,2~)
  ~1.4142135623730947,1.4142135623730954~ 
To test the speed of the convergence, the easiest method is
to look at the Taylor expansion of the iterator around the fix
point.

In our case, the the iterator clearly converges qudratically
to the limit.
>&taylor(f(x),x,sqrt(2),2)
                                    2
                       (x - sqrt(2))
                       -------------- + sqrt(2)
                             3/2
                            2

Since our case is so easy, we can compute the Taylor expansion
explicitely.
>&ratsimp(f(x)-sqrt(2))
                            2    3/2
                           x  - 2    x + 2
                           ---------------
                                 2 x

If we try to apply the newton method to a function with a double
zero, we do not get a good result.
>mxmnewton("x^3-2*x^2+x",2)
      1.000000005479078 
Let us define a function like this in Maxima.
>function g(x) &= x^3-2*x^2+x
                             3      2
                            x  - 2 x  + x

Then plot it in Euler.
>plot2d(mxm("g(x)"),0,2); insimg;

Convergence of the Newton Method

Now, the newton interation is no longer good.
>iterate(&x-g(x)/diff(g(x),x),2,10)
Column 1 to 2:
                    1.6       1.347368421052632 
Column 3 to 4:
      1.193516663631397       1.104014284855782 
Column 5 to 6:
      1.054346842021409       1.027856158851746 
Column 7 to 8:
      1.014114290149139       1.007105915824406 
Column 9 to 10:
      1.003565448288775       1.001785885343114 
The iteration speed can be fixed by multiplying with the order
of the zero.
>iterate(&x-2*g(x)/diff(g(x),x),2,10)
Column 1 to 2:
                    1.2       1.015384615384616 
Column 3 to 4:
      1.000115673799884        1.00000000668709 
Column 5 to 6:
     0.9999999734821227       1.000000002788999 
Column 7 to 8:
      1.000000002788999       1.000000002788999 
Column 9 to 10:
      1.000000002788999       1.000000002788999 
Alternatively, we can apply the Newton method to the function
g/g'. We use to define instead of := to get an evaluated expression
for dg.
>function dg(x) &= g(x)/diff(g(x),x)
                             3      2
                            x  - 2 x  + x
                            --------------
                               2
                            3 x  - 4 x + 1

dg has a single zero at 1 now.
>plot2d(&dg(x),0.8,2); insimg;

Convergence of the Newton Method

Indeed the iteration is good, but not the result.
>iterate(&x-dg(x)/diff(dg(x),x),2,10)
Column 1 to 2:
     0.8888888888888888      0.9922480620155019 
Column 3 to 4:
     0.9999694833531766      0.9999999995340271 
Column 5 to 6:
     0.9999999995340271      0.9999999995340271 
Column 7 to 8:
     0.9999999995340271      0.9999999995340271 
Column 9 to 10:
     0.9999999995340271      0.9999999995340271 
Why is the result not good? The reason is that the evaluation
of the polynomial is already 0, even if we are away from 1.
>t=1.000000001; polyval([0,1,-2,1],t)
                      0 
It does not help to scale the polynomial.
>polyval(10^10*[0,1,-2,1],t)/10^10
                      0 
However, the foolowing trick helps. We use the exact evaluation of
polynomial in Euler. This feature makes use of the exact scalar
product.
>xpolyval(10^10*[0,1,-2,1],t)/10^10
 1.000000166480752e-018 
This is not the exact result, of course.
>&g(x) with x=1+1/100000000
                              100000001
                      -------------------------
                      1000000000000000000000000

Or with Maxima big float numbers.
>fpprec::=50; &bfloat(g(x)) with x=1.000000001b0
       1.0000000010000000000000000000000006141818534508651b-18

Still, we try to base our Newton method on this trick.
>function g(x) := xpolyval(10^10*[0,1,-2,1],x)/10^10
The following function evaluates the derivative at compile time.
>function dg(x) &= diff(g(x),x)
                               2
                            3 x  - 4 x + 1

This gives a very good Newton method.
>iterate("x-2*g(x)/dg(x)",2,5)
Column 1 to 2:
                    1.2       1.015384615384616 
Column 3 to 4:
      1.000115673799888       1.000000006689053 
Column 5:
                      1 
Note that we could also transfer the zero to 0.
>function h(x) &= expand(subst(x+1,x,g(x)))
                                3    2
                               x  + x

This works very well too, even without any exact evaluation.
>iterate(&x-2*h(x)/diff(h(x),x),1,5)
Column 1 to 2:
                    0.2     0.01538461538461536 
Column 3 to 4:
  0.0001156737998843256  6.689053367505012e-009 
Column 5:
  2.23717167432336e-017 

Examples Homepage