Skip to content

Commit 1412b1f

Browse files
eeshan9815eeshan9815
committed
Major robustness improvements and hard test cases added
1 parent aa50d91 commit 1412b1f

File tree

2 files changed

+111
-42
lines changed

2 files changed

+111
-42
lines changed

‎src/newton1d.jl‎

Lines changed: 87 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,16 @@ Optional keyword arguments give the tolerances `reltol` and `abstol`.
55
and a `debug` boolean argument that prints out diagnostic information."""
66

77
function newton1d{T}(f::Function, f′::Function, x::Interval{T};
8-
reltol=10eps(T), abstol=eps(T), debug=false)
8+
reltol=eps(T), abstol=eps(T), debug=false, debugroot=false)
99

1010
L = Interval{T}[]
1111

1212
R = Root{Interval{T}}[]
13+
reps = reps1 = 0
1314

1415
push!(L, x)
15-
16+
initial_width = ∞
17+
X = emptyinterval(T)
1618
while !isempty(L)
1719
X = pop!(L)
1820

@@ -28,56 +30,111 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
2830
debug && println("0 ∉ f′(X)")
2931

3032
while true
31-
m = mid(X)
32-
N = m - (f(Interval(m)) / f′(X))
33-
34-
debug && (print("Newton step: "); @show (X, X ∩ N))
3533

34+
m = mid(X)
35+
N = m - (f(interval(m)) / f′(X))
36+
37+
debug && (print("Newton step1: "); @show (X, X ∩ N))
38+
if X == X ∩ N
39+
reps1 += 1
40+
if reps1 > 20
41+
reps1 = 0
42+
break
43+
end
44+
end
3645
X = X ∩ N
3746

38-
if isempty(X)
47+
if (isempty(X) || diam(X) == 0)
3948
break
4049

41-
elseif 0 ∈ f(Interval(prevfloat(m), nextfloat(m)))
42-
push!(R, Root(X, :unique))
43-
44-
debug && @show "Root found", X
50+
elseif 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X))))
51+
n = fa = fb = 0
52+
root_exist = true
53+
while (n < 4 && (fa == 0 || fb == 0))
54+
if fa == 0
55+
if 0 ∈ f(interval(prevfloat(X.lo), nextfloat(X.lo)))
56+
fa = 1
57+
else
58+
N = X.lo - (f(interval(X.lo)) / f′(X))
59+
X = X ∩ N
60+
if (isempty(X) || diam(X) == 0)
61+
root_exist = false
62+
break
63+
end
64+
end
65+
end
66+
if fb == 0
67+
if 0 ∈ f(interval(prevfloat(X.hi), nextfloat(X.hi)))
68+
fb = 1
69+
else
70+
if 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X))))
71+
N = X.hi - (f(interval(X.hi)) / f′(X))
72+
else
73+
N = mid(X) - (f(interval(mid(X))) / f′(X))
74+
end
75+
X = X ∩ N
76+
if (isempty(X) || diam(X) == 0)
77+
root_exist = false
78+
break
79+
end
80+
end
81+
end
82+
N = mid(X) - (f(interval(mid(X))) / f′(X))
83+
X = X ∩ N
84+
if (isempty(X) || diam(X) == 0)
85+
root_exist = false
86+
break
87+
end
88+
n += 1
89+
end
90+
if root_exist
91+
push!(R, Root(X, :unique))
92+
debugroot && @show "Root found", X
93+
end
4594

4695
break
4796
end
4897
end
4998

5099
else
51-
# 0 ∈ f'(X)
52-
100+
if diam(X) == initial_width
101+
reps += 1
102+
if reps > 10
103+
push!(R, Root(X, :unknown))
104+
debugroot && @show "Repititive root found", X
105+
reps = 0
106+
continue
107+
end
108+
end
109+
initial_width = diam(X)
53110
debug && println("0 ∈ f'(X)")
54111

55112
expansion_pt = Inf
56113
# expansion point for the newton step might be m, X.lo or X.hi according to some conditions
57114

58-
if 0 ∈ f(Interval(mid(X)))
115+
if 0 ∈ f(interval(prevfloat(mid(X)), nextfloat(mid(X))))
59116
# 0 ∈ fⁱ(x)
60117

61118
debug && println("0 ∈ fⁱ(x)")
62119

63-
if 0 ∉ f(Interval(X.lo))
120+
if 0 ∉ f(interval(prevfloat(X.lo), nextfloat(X.lo)))
64121
expansion_pt = X.lo
65122

66-
elseif 0 ∉ f(Interval(X.hi))
123+
elseif 0 ∉ f(interval(prevfloat(X.hi), nextfloat(X.hi)))
67124
expansion_pt = X.hi
68125

69126
else
70-
x1 = mid(Interval(X.lo, mid(X)))
71-
x2 = mid(Interval(mid(X), X.hi))
72-
if 0 ∉ f(Interval(x1)) || 0 ∉ f(Interval(x2))
73-
push!(L, Interval(X.lo, m))
74-
push!(L, Interval(m, X.hi))
127+
x1 = mid(interval(X.lo, mid(X)))
128+
x2 = mid(interval(mid(X), X.hi))
129+
if 0 ∉ f(interval(prevfloat(x1), nextfloat(x1))) || 0 ∉ f(interval(prevfloat(x2), nextfloat(x2)))
130+
push!(L, interval(X.lo, m))
131+
push!(L, interval(m, X.hi))
75132
continue
76133

77134
else
78135
push!(R, Root(X, :unknown))
79136

80-
debug && @show "Multiple root found", X
137+
debugroot && @show "Multiple root found", X
81138

82139
continue
83140
end
@@ -91,7 +148,7 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
91148
if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol
92149
push!(R, Root(X, :unknown))
93150

94-
debug && @show "Tolerance root found", X
151+
debugroot && @show "Tolerance root found", X
95152

96153
continue
97154
end
@@ -102,28 +159,27 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
102159
expansion_pt = mid(X)
103160
end
104161

105-
initial_width = diam(X)
106162

107-
a = f(Interval(expansion_pt))
163+
a = f(interval(expansion_pt))
108164
b = f′(X)
109165

110166
if 0 < b.hi && 0 > b.lo && 0 ∉ a
111167
if a.hi < 0
112-
push!(L, X ∩ (expansion_pt - Interval(-Inf, a.hi / b.hi)))
113-
push!(L, X ∩ (expansion_pt - Interval(a.hi / b.lo, Inf)))
168+
push!(L, X ∩ (expansion_pt - interval(-Inf, a.hi / b.hi)))
169+
push!(L, X ∩ (expansion_pt - interval(a.hi / b.lo, Inf)))
114170

115171
elseif a.lo > 0
116-
push!(L, X ∩ (expansion_pt - Interval(-Inf, a.lo / b.lo)))
117-
push!(L, X ∩ (expansion_pt - Interval(a.lo / b.hi, Inf)))
172+
push!(L, X ∩ (expansion_pt - interval(-Inf, a.lo / b.lo)))
173+
push!(L, X ∩ (expansion_pt - interval(a.lo / b.hi, Inf)))
118174

119175
end
120176

121177
continue
122178

123179
else
124-
N = expansion_pt - (f(Interval(expansion_pt))/f′(X))
180+
N = expansion_pt - (f(interval(expansion_pt))/f′(X))
125181

126-
debug && (print("Newton step: "); @show (X, X ∩ N))
182+
debug && (print("Newton step2: "); @show (X, X ∩ N))
127183

128184
X = X ∩ N
129185
m = mid(X)
@@ -133,11 +189,6 @@ function newton1d{T}(f::Function, f′::Function, x::Interval{T};
133189
end
134190
end
135191

136-
if diam(X) > initial_width/2
137-
push!(L, Interval(m, X.hi))
138-
X = Interval(X.lo, m)
139-
end
140-
141192
push!(L, X)
142193
end
143194
end

‎test/newton1d.jl‎

Lines changed: 24 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -18,34 +18,42 @@ three_halves_pi = 3*big_pi/2
1818

1919
@testset "Testing newton1d" begin
2020

21-
f(x) = exp(x^2) - cos(x)
21+
f(x) = exp(x^2) - cos(x) #double root
2222
f′(x) = 2*x*exp(x^2) + sin(x)
23-
f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24
23+
f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 #four unique roots
2424
f1′(x) = 4x^3 - 30x^2 + 70x - 50
25-
f2(x) = 4567x^2 - 9134x + 4567
25+
f2(x) = 4567x^2 - 9134x + 4567 #double root
2626
f2′(x) = 9134x - 9134
27-
f3(x) = (x^2 - 2)^2
27+
f3(x) = (x^2 - 2)^2 #two double roots
2828
f3′(x) = 4x * (x^2 - 2)
29+
f4(x) = sin(x) - x #triple root at 0
30+
f4′(x) = cos(x) - 1
31+
f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots
32+
f5′(x) = 8x * (-3 + 2x^2) * (2 - 3x^2 + x^4)^3
2933
for autodiff in (false, true)
3034
if autodiff
3135
rts1 = newton1d(sin, -5..5)
3236
rts2 = newton1d(f, -∞..∞)
3337
rts3 = newton1d(f1, -10..10)
3438
rts4 = newton1d(f2, -10..11)
3539
rts5 = newton1d(f3, -10..10)
40+
rts6 = newton1d(f4, -10..10)
41+
rts7 = newton1d(f5, -10..10)
3642

3743
else
3844
rts1 = newton1d(sin, cos, -5..5)
3945
rts2 = newton1d(f, f′, -∞..∞)
4046
rts3 = newton1d(f1, f1′, -10..10)
4147
rts4 = newton1d(f2, f2′, -10..11)
4248
rts5 = newton1d(f3, f3′, -10..10)
49+
rts6 = newton1d(f4, f4′, -10..10)
50+
rts7 = newton1d(f5, f5′, -10..10, reltol=0)
4351
end
4452

4553
@test length(rts1) == 3
4654
L = [ -pi_interval(Float64), 0..0, pi_interval(Float64)]
4755
for i = 1:length(rts1)
48-
@test L[i] == rts1[i].interval && :unique == rts1[i].status
56+
@test L[i] in rts1[i].interval && :unique == rts1[i].status
4957
end
5058

5159
@test length(rts2) == 1
@@ -60,9 +68,19 @@ three_halves_pi = 3*big_pi/2
6068
@test length(rts4) == 1
6169
@test 1 in rts4[1].interval && :unknown == rts4[1].status
6270

63-
L1 = [-sqrt(2), -sqrt(2), sqrt(2), sqrt(2)]
71+
L1 = [-sqrt(2), sqrt(2)]
6472
for i = 1:length(rts5)
6573
@test L1[i] in rts5[i].interval && :unknown == rts5[i].status
6674
end
75+
76+
@test length(rts6) == 1
77+
@test 0 in rts6[1].interval && :unknown == rts6[1].status
78+
79+
@test length(rts7) == 4
80+
L = [-sqrt(2), -1, 1, sqrt(2)]
81+
for i = 1:length(rts7)
82+
@test L[i] in rts7[i].interval && :unknown == rts7[i].status
83+
end
84+
6785
end
6886
end

0 commit comments

Comments
 (0)