Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

Example 1

Suppose we wish to find a root of the polynomial f(x) = x⁷ + 3x⁶ + 7x⁵ + x⁴ + 5x³ + 2x² + 5x + 5. There exist closed form solutions to the roots of polynomials for quartics and below, and this is a degree seven polynomial, so thus we must use a numerical technique. We will use x₀ = 0, x₁ = -0.1, and x₂ = -0.2 as our three initial approximations. We will let the two values ε_step = 0.0001 and ε_abs = 0.0001 and we will halt after N = 100 iterations.

We will use six decimal digit arithmetic to find a solution.

The easiest way to present such information is in tabular format:

n	x_{n - 2}	x_{n - 1}	x_n	x_{n + 1}	\|f(x_{n + 1})\|	\|x_{n + 1} - x_n\|
2	0.0	-0.1	-0.2	-1.14874	13.6910	0.94874
3	-0.1	-0.2	-1.14874	-0.568100	1.65995	0.580640
4	-0.2	-1.14874	-0.568100	-0.669622	0.51613	0.101522
5	-1.14874	-0.568100	-0.669622	-0.702849	0.05806	0.033227
6	-0.568100	-0.669622	-0.702849	-0.706858	0.00047	0.004009
7	-0.669622	-0.702849	-0.706858	-0.706826	0.00001	0.000032

Thus, with the last step, both halting conditions are met, and therefore, after six iterations, our approximation to the root is -0.706826 .

Example 2

The following demonstrates the first six iterations of Müller's method in Matlab. Suppose we wish to find a root of the same polynomial

p(x) = x⁷ + 3x⁶ + 7x⁵ + x⁴ + 5x³ + 2x² + 5x + 5

starting with the same three initial approximations x₀ = 0, x₁ = -0.1, and x₂ = -0.2.

The first formula in red is the root of the quadratic polynomial which is added onto the middle approximation x(2).

>> p = [1 3 7 1 5 2 5 5]
p =
1      3      7      1      5      2      5      5

>> x = [0.0 -0.1 -0.2]'      % first three approximations
x =
   0.00000
  -0.10000
  -0.20000

>> % 1st iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
   0.01000   0.10000   1.00000
   0.00000   0.00000   1.00000
   0.01000  -0.10000   1.00000

>> y = polyval( p, x )
y =
  5.00000
  4.51503
  4.03954

>> c = M \ y
c =
  0.47367
  4.80230
  4.51503

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'
x =
  -0.10000
  -0.20000
  -1.14864

>> % 2nd iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
   0.01000   0.10000   1.00000
   0.00000   0.00000   1.00000
   0.89992  -0.94864   1.00000

>> y = polyval( p, x )                                                          
y =
    4.5150
    4.0395
  -13.6858

>> c = M \ y                                                                    
c =
  -13.2838
    6.0833
    4.0395

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'         
x =
  -0.20000
  -1.14864
  -0.56812

>> % 3rd iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
  0.89992  0.94864  1.00000
  0.00000  0.00000  1.00000
  0.33701  0.58052  1.00000

>> y = polyval( p, x )                                                          
y =
    4.0395
  -13.6858
    1.6597

>> c = M \ y                                                                    
c =
  -21.0503
   38.6541
  -13.6858

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'         
x =
  -1.14864
  -0.56812
  -0.66963

>> % 4th iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
   0.33701  -0.58052   1.00000
   0.00000   0.00000   1.00000
   0.01030  -0.10151   1.00000

>> y = polyval( p, x )                                                          
y =
  -13.6858
    1.6597
    0.5160

>> c = M \ y                                                                    
c =
  -31.6627
    8.0531
    1.6597

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'         
x =
  -0.56812
  -0.66963
  -0.70285

>> % 5th iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
   0.01030   0.10151   1.00000
   0.00000   0.00000   1.00000
   0.00110  -0.03322   1.00000

>> y = polyval( p, x )                                                          
y =
  1.65973
  0.51602
  0.05802

>> c = M \ y                                                                    
c =
  -18.6991
   13.1653
    0.5160

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'         
x =
  -0.66963
  -0.70285
  -0.70686

>> % 6th iteration ---------------------------------------
>> M = [(x(1) - x(2))^2,  x(1) - x(2),  1
      0,               0,       1
(x(3) - x(2))^2,  x(3) - x(2),  1]
M =
   0.00110   0.03322   1.00000
   0.00000   0.00000   1.00000
   0.00002  -0.00401   1.00000

>> y = polyval( p, x )                                                          
y =
   0.51602
   0.05802
  -0.00046

>> c = M \ y                                                                    
c =
  -21.8018
   14.5107
    0.0580

>> x = [x(2), x(3), x(2) - 2*c(3)/(c(2) + sqrt(c(2)^2 - 4*c(1)*c(3)))]'         
x =
  -0.70285
  -0.70686
  -0.70683

The list of all iterations using Matlab are:

   0.000000000000000
  -0.100000000000000
  -0.200000000000000
  -1.148643697414111
  -0.568122032631211
  -0.669630566165950
  -0.702851144883234
  -0.706857484921269
  -0.706825973130949
  -0.706825980788168
  -0.706825980788170

Note how convergence speeds up after none of the first three initial approximations are used to calculate the next iterate?

Topic 10.6: Müller's Method (Examples)

Example 1

Example 2