Newton method with double root

Consider the function $$ f(x) = (x-1)^2 \sin(x) $$ for which $x=1$ is a double root.

In [35]:
from numpy import sin,cos,linspace,zeros,abs
from matplotlib.pyplot import plot,xlabel,ylabel
In [36]:
def f(x):
    return (x-1.0)**2 * sin(x)

def df(x):
    return 2.0*(x-1.0)*sin(x) + (x-1.0)**2 * cos(x)
In [37]:
x = linspace(0.0,2.0,100)
plot(x,f(x))
xlabel('x')
ylabel('f(x)')
Out[37]:
<matplotlib.text.Text at 0x10dcf6850>
In [38]:
def newton(x0,m):
    n = 50
    x = zeros(50)
    x[0] = x0
    print "%6d %24.14e" % (0,x[0])
    for i in range(1,50):
        x[i] = x[i-1] - m*f(x[i-1])/df(x[i-1])
        if i > 1:
            r = (x[i] - x[i-1])/(x[i-1]-x[i-2])
        else:
            r = 0.0
        print "%6d %24.14e %14.6e" % (i,x[i],r)
        if abs(f(x[i])) < 1.0e-16:
            break

We first apply the standard newton method.

In [39]:
newton(2.0,1)
     0     2.00000000000000e+00
     1     1.35163555744248e+00   0.000000e+00
     2     1.18244356861394e+00   2.609520e-01
     3     1.09450383817604e+00   5.197630e-01
     4     1.04837639635805e+00   5.245347e-01
     5     1.02452044172239e+00   5.171749e-01
     6     1.01235093391487e+00   5.101245e-01
     7     1.00619920246051e+00   5.055037e-01
     8     1.00310567444820e+00   5.028711e-01
     9     1.00155437342780e+00   5.014666e-01
    10     1.00077757303341e+00   5.007412e-01
    11     1.00038888338214e+00   5.003726e-01
    12     1.00019446594325e+00   5.001868e-01
    13     1.00009723903915e+00   5.000935e-01
    14     1.00004862103702e+00   5.000468e-01
    15     1.00002431089794e+00   5.000234e-01
    16     1.00001215554384e+00   5.000117e-01
    17     1.00000607779564e+00   5.000059e-01
    18     1.00000303890375e+00   5.000029e-01
    19     1.00000151945336e+00   5.000015e-01
    20     1.00000075972705e+00   5.000007e-01
    21     1.00000037986362e+00   5.000004e-01
    22     1.00000018993183e+00   5.000002e-01
    23     1.00000009496592e+00   5.000001e-01
    24     1.00000004748296e+00   5.000000e-01
    25     1.00000002374148e+00   5.000000e-01
    26     1.00000001187074e+00   5.000000e-01
    27     1.00000000593537e+00   5.000000e-01

Newton method is converging linearly. Now try the modified Newton method.

In [40]:
newton(2.0,2)
     0     2.00000000000000e+00
     1     7.03271114884954e-01   0.000000e+00
     2     1.06293359720705e+00  -2.773614e-01
     3     1.00108318855754e+00  -1.719679e-01
     4     1.00000037565568e+00   1.750696e-02
     5     1.00000000000005e+00   3.469257e-04