Skip to content
This repository was archived by the owner on May 15, 2025. It is now read-only.

Commit 8eb8afc

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #56 from huiyuxie/main
Fix Issue #34: Add new method Alefeld, Potra, Shi (1995)
2 parents df85c3f + e44bab2 commit 8eb8afc

File tree

4 files changed

+200
-3
lines changed

4 files changed

+200
-3
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,4 +43,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
4343
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
4444

4545
[targets]
46-
test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
46+
test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]

‎src/SimpleNonlinearSolve.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ include("brent.jl")
3838
include("dfsane.jl")
3939
include("ad.jl")
4040
include("halley.jl")
41+
include("alefeld.jl")
4142

4243
import SnoopPrecompile
4344

@@ -59,13 +60,13 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
5960
=#
6061

6162
prob_brack = IntervalNonlinearProblem{false}((u, p) -> u * u - p, T.((0.0, 2.0)), T(2))
62-
for alg in (Bisection, Falsi, Ridder, Brent)
63+
for alg in (Bisection, Falsi, Ridder, Brent, Alefeld)
6364
solve(prob_brack, alg(), abstol = T(1e-2))
6465
end
6566
end end
6667

6768
# DiffEq styled algorithms
6869
export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Halley, Klement,
69-
Ridder, SimpleNewtonRaphson, SimpleTrustRegion
70+
Ridder, SimpleNewtonRaphson, SimpleTrustRegion, Alefeld
7071

7172
end # module

‎src/alefeld.jl

Lines changed: 184 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,184 @@
1+
"""
2+
`Alefeld()`
3+
4+
An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
5+
6+
The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than
7+
algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
8+
"""
9+
struct Alefeld <: AbstractBracketingAlgorithm end
10+
11+
function SciMLBase.solve(prob::IntervalNonlinearProblem,
12+
alg::Alefeld, args...; abstol = nothing,
13+
reltol = nothing,
14+
maxiters = 1000, kwargs...)
15+
16+
f = Base.Fix2(prob.f, prob.p)
17+
a, b = prob.tspan
18+
c = a - (b - a) / (f(b) - f(a)) * f(a)
19+
20+
fc = f(c)
21+
if iszero(fc)
22+
return SciMLBase.build_solution(prob, alg, c, fc;
23+
retcode = ReturnCode.Success,
24+
left = a,
25+
right = b)
26+
end
27+
a, b, d = _bracket(f, a, b, c)
28+
e = zero(a) # Set e as 0 before iteration to avoid a non-value f(e)
29+
30+
# Begin of algorithm iteration
31+
for i in 2:maxiters
32+
# The first bracketing block
33+
f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
34+
if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
35+
c = _newton_quadratic(f, a, b, d, 2)
36+
else
37+
c = _ipzero(f, a, b, d, e)
38+
if (c - a) * (c - b) ≥ 0
39+
c = _newton_quadratic(f, a, b, d, 2)
40+
end
41+
end
42+
ē, fc = d, f(c)
43+
(a == c || b == c) &&
44+
return SciMLBase.build_solution(prob, alg, c, fc;
45+
retcode = ReturnCode.FloatingPointLimit,
46+
left = a,
47+
right = b)
48+
iszero(fc) &&
49+
return SciMLBase.build_solution(prob, alg, c, fc;
50+
retcode = ReturnCode.Success,
51+
left = a,
52+
right = b)
53+
ā, b̄, d̄ = _bracket(f, a, b, c)
54+
55+
# The second bracketing block
56+
f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
57+
if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
58+
c = _newton_quadratic(f, ā, b̄, d̄, 3)
59+
else
60+
c = _ipzero(f, ā, b̄, d̄, ē)
61+
if (c - ā) * (c - b̄) ≥ 0
62+
c = _newton_quadratic(f, ā, b̄, d̄, 3)
63+
end
64+
end
65+
fc = f(c)
66+
(ā == c || b̄ == c) &&
67+
return SciMLBase.build_solution(prob, alg, c, fc;
68+
retcode = ReturnCode.FloatingPointLimit,
69+
left = ā,
70+
right = b̄)
71+
iszero(fc) &&
72+
return SciMLBase.build_solution(prob, alg, c, fc;
73+
retcode = ReturnCode.Success,
74+
left = ā,
75+
right = b̄)
76+
ā, b̄, d̄ = _bracket(f, ā, b̄, c)
77+
78+
# The third bracketing block
79+
if abs(f(ā)) < abs(f(b̄))
80+
u = ā
81+
else
82+
u = b̄
83+
end
84+
c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
85+
if (abs(c - u)) > 0.5 * (b̄ - ā)
86+
c = 0.5 * (ā + b̄)
87+
end
88+
fc = f(c)
89+
(ā == c || b̄ == c) &&
90+
return SciMLBase.build_solution(prob, alg, c, fc;
91+
retcode = ReturnCode.FloatingPointLimit,
92+
left = ā,
93+
right = b̄)
94+
iszero(fc) &&
95+
return SciMLBase.build_solution(prob, alg, c, fc;
96+
retcode = ReturnCode.Success,
97+
left = ā,
98+
right = b̄)
99+
ā, b̄, d = _bracket(f, ā, b̄, c)
100+
101+
# The last bracketing block
102+
if b̄ - ā < 0.5 * (b - a)
103+
a, b, e = ā, b̄, d̄
104+
else
105+
e = d
106+
c = 0.5 * (ā + b̄)
107+
fc = f(c)
108+
(ā == c || b̄ == c) &&
109+
return SciMLBase.build_solution(prob, alg, c, fc;
110+
retcode = ReturnCode.FloatingPointLimit,
111+
left = ā,
112+
right = b̄)
113+
iszero(fc) &&
114+
return SciMLBase.build_solution(prob, alg, c, fc;
115+
retcode = ReturnCode.Success,
116+
left = ā,
117+
right = b̄)
118+
a, b, d = _bracket(f, ā, b̄, c)
119+
end
120+
end
121+
122+
# Reassign the value a, b, and c
123+
if b == c
124+
b = d
125+
elseif a == c
126+
a = d
127+
end
128+
fc = f(c)
129+
130+
# Reuturn solution when run out of max interation
131+
return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
132+
left = a, right = b)
133+
end
134+
135+
# Define subrotine function bracket, check fc before bracket to return solution
136+
function _bracket(f::F, a, b, c) where F
137+
if iszero(f(c))
138+
ā, b̄, d = a, b, c
139+
else
140+
if f(a) * f(c) < 0
141+
ā, b̄, d = a, c, b
142+
elseif f(b) * f(c) < 0
143+
ā, b̄, d = c, b, a
144+
end
145+
end
146+
147+
return ā, b̄, d
148+
end
149+
150+
# Define subrotine function newton quadratic, return the approximation of zero
151+
function _newton_quadratic(f::F, a, b, d, k) where F
152+
A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a)
153+
B = (f(b) - f(a)) / (b - a)
154+
155+
if iszero(A)
156+
return a - (1 / B) * f(a)
157+
elseif A * f(a) > 0
158+
rᵢ₋₁ = a
159+
else
160+
rᵢ₋₁ = b
161+
end
162+
163+
for i in 1:k
164+
rᵢ = rᵢ₋₁ - (f(a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) / (B + A * (2 * rᵢ₋₁ - a - b))
165+
rᵢ₋₁ = rᵢ
166+
end
167+
168+
return rᵢ₋₁
169+
end
170+
171+
# Define subrotine function ipzero, also return the approximation of zero
172+
function _ipzero(f::F, a, b, c, d) where F
173+
Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
174+
Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
175+
Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
176+
D₂₁ = (b - c) * f(c) / (f(c) - f(b))
177+
D₃₁ = (a - b) * f(b) / (f(b) - f(a))
178+
Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
179+
Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
180+
D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
181+
Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
182+
183+
return a + Q₃₁ + Q₃₂ + Q₃₃
184+
end

‎test/basictests.jl

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -210,6 +210,18 @@ for p in 1.1:0.1:100.0
210210
@test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
211211
end
212212

213+
# Alefeld
214+
g = function (p)
215+
probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
216+
sol = solve(probN, Alefeld())
217+
return sol.u
218+
end
219+
220+
for p in 1.1:0.1:100.0
221+
@test g(p) ≈ sqrt(p)
222+
@test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
223+
end
224+
213225
f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
214226
t = (p) -> [sqrt(p[2] / p[1])]
215227
p = [0.9, 50.0]

0 commit comments

Comments
 (0)