Repeated Roots

This short worksheet demonstrates some of the numerical problems associated with finding roots where there is a multiple root.

%auto
def nestlist(f,start,times):
    result=[start]
    for i in range(times):
        result.append(f(result[-1]))
    return result
def consecutive_numbers(start,num):
    return nestlist(lambda x: x.nextabove(), start, num-1)

consecutive_numbers(RR('1'),20)

[1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000, 1.00000000000000]

[i.str(truncate=False) for i in consecutive_numbers(RR('1'),20)]

['1.0000000000000000', '1.0000000000000002', '1.0000000000000004', '1.0000000000000007', '1.0000000000000009', '1.0000000000000011', '1.0000000000000013', '1.0000000000000016', '1.0000000000000018', '1.0000000000000020', '1.0000000000000022', '1.0000000000000024', '1.0000000000000027', '1.0000000000000029', '1.0000000000000031', '1.0000000000000033', '1.0000000000000036', '1.0000000000000038', '1.0000000000000040', '1.0000000000000042']

f(x)=x^3-2*x^2+4/3*x-8/27
f(x); f(x).factor()

x^3 - 2*x^2 + 4/3*x - 8/27
1/27*(3*x - 2)^3

Note that 2/3 is an exact zero (we give an exact input, so Sage uses exact arithmetic).

f(2/3)

0

However, note that in floating point arithmetic, 0.66666410000000009 is also a root.

a=RR('0.66666410000000009')
b=f(a)
print RR(b).str(truncate=False), b.is_zero()

0.00000000000000000 True

In fact, let's see how many roots there are in the 20 floating point numbers following $a$ .

data=[(x,f(x)) for x in consecutive_numbers(a,20)]

html.table([['$x$', '$f(x)$', 'Root?']]+[[RR(x).str(truncate=False),RR(y).str(truncate=False), RR(y).is_zero()] for x,y in data], header=True)

x	f(x)	Root?
0.66666410000000009	0.00000000000000000	{\rm True}
0.66666410000000020	-1.1102230246251565e-16	{\rm False}
0.66666410000000031	-1.1102230246251565e-16	{\rm False}
0.66666410000000043	0.00000000000000000	{\rm True}
0.66666410000000054	-1.1102230246251565e-16	{\rm False}
0.66666410000000065	-1.1102230246251565e-16	{\rm False}
0.66666410000000076	0.00000000000000000	{\rm True}
0.66666410000000087	-1.1102230246251565e-16	{\rm False}
0.66666410000000098	-1.1102230246251565e-16	{\rm False}
0.66666410000000109	0.00000000000000000	{\rm True}
0.66666410000000120	-1.1102230246251565e-16	{\rm False}
0.66666410000000131	-1.1102230246251565e-16	{\rm False}
0.66666410000000142	0.00000000000000000	{\rm True}
0.66666410000000154	-1.1102230246251565e-16	{\rm False}
0.66666410000000165	-1.1102230246251565e-16	{\rm False}
0.66666410000000176	0.00000000000000000	{\rm True}
0.66666410000000187	-1.1102230246251565e-16	{\rm False}
0.66666410000000198	-1.1102230246251565e-16	{\rm False}
0.66666410000000209	0.00000000000000000	{\rm True}
0.66666410000000220	-1.1102230246251565e-16	{\rm False}

points([p[0:2] for p in data],ticks=[[p[0] for p in data[::6]],None])

Wow, that is a lot of roots for a cubic polynomial!

Here is the same example, just written using Horner's form.

f(x)=(((x-2)*x+4/3)*x-8/27)
f(x); f(x).factor()

1/3*(3*(x - 2)*x + 4)*x - 8/27
1/27*(3*x - 2)^3

f(2/3)

0

a=RR('0.66666410000000009')
b=f(a)
print RR(b).str(truncate=False), b.is_zero()

1.1102230246251565e-16 False

data=[(x,f(x)) for x in consecutive_numbers(a,20)]

html.table([['$x$', '$f(x)$', 'Root?']]+[[RR(x).str(truncate=False),RR(y).str(truncate=False), RR(y).is_zero()] for x,y in data], header=True)

x	f(x)	Root?
0.66666410000000009	1.1102230246251565e-16	{\rm False}
0.66666410000000020	5.5511151231257827e-17	{\rm False}
0.66666410000000031	-1.1102230246251565e-16	{\rm False}
0.66666410000000043	-5.5511151231257827e-17	{\rm False}
0.66666410000000054	0.00000000000000000	{\rm True}
0.66666410000000065	-5.5511151231257827e-17	{\rm False}
0.66666410000000076	0.00000000000000000	{\rm True}
0.66666410000000087	5.5511151231257827e-17	{\rm False}
0.66666410000000098	1.1102230246251565e-16	{\rm False}
0.66666410000000109	-5.5511151231257827e-17	{\rm False}
0.66666410000000120	0.00000000000000000	{\rm True}
0.66666410000000131	5.5511151231257827e-17	{\rm False}
0.66666410000000142	1.1102230246251565e-16	{\rm False}
0.66666410000000154	5.5511151231257827e-17	{\rm False}
0.66666410000000165	-1.1102230246251565e-16	{\rm False}
0.66666410000000176	-5.5511151231257827e-17	{\rm False}
0.66666410000000187	0.00000000000000000	{\rm True}
0.66666410000000198	-5.5511151231257827e-17	{\rm False}
0.66666410000000209	0.00000000000000000	{\rm True}
0.66666410000000220	5.5511151231257827e-17	{\rm False}

points([p[0:2] for p in data],ticks=[[p[0] for p in data[::6]],None])

Notice that we still have numerical problems!

line([p[0:2] for p in data],ticks=[[p[0] for p in data[::6]],None])

This work by Jason Grout is licensed under a Creative Commons Attribution-Share Alike 3.0 United States License.