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

Commit 2a43225

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #59 from ranocha/hr/format
format
2 parents f4e1188 + 78f1461 commit 2a43225

File tree

1 file changed

+43
-42
lines changed

1 file changed

+43
-42
lines changed

‎src/alefeld.jl‎

Lines changed: 43 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,13 @@ algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal
99
struct Alefeld <: AbstractBracketingAlgorithm end
1010

1111
function SciMLBase.solve(prob::IntervalNonlinearProblem,
12-
alg::Alefeld, args...; abstol = nothing,
13-
reltol = nothing,
14-
maxiters = 1000, kwargs...)
15-
12+
alg::Alefeld, args...; abstol = nothing,
13+
reltol = nothing,
14+
maxiters = 1000, kwargs...)
1615
f = Base.Fix2(prob.f, prob.p)
1716
a, b = prob.tspan
1817
c = a - (b - a) / (f(b) - f(a)) * f(a)
19-
18+
2019
fc = f(c)
2120
(a == c || b == c) &&
2221
return SciMLBase.build_solution(prob, alg, c, fc;
@@ -25,7 +24,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
2524
right = b)
2625
iszero(fc) &&
2726
return SciMLBase.build_solution(prob, alg, c, fc;
28-
retcode = ReturnCode.Success,
27+
retcode = ReturnCode.Success,
2928
left = a,
3029
right = b)
3130
a, b, d = _bracket(f, a, b, c)
@@ -37,47 +36,47 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
3736
f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
3837
if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
3938
c = _newton_quadratic(f, a, b, d, 2)
40-
else
39+
else
4140
c = _ipzero(f, a, b, d, e)
4241
if (c - a) * (c - b) ≥ 0
4342
c = _newton_quadratic(f, a, b, d, 2)
4443
end
45-
end
44+
end
4645
ē, fc = d, f(c)
47-
(a == c || b == c) &&
46+
(a == c || b == c) &&
4847
return SciMLBase.build_solution(prob, alg, c, fc;
4948
retcode = ReturnCode.FloatingPointLimit,
50-
left = a,
51-
right = b)
49+
left = a,
50+
right = b)
5251
iszero(fc) &&
5352
return SciMLBase.build_solution(prob, alg, c, fc;
54-
retcode = ReturnCode.Success,
55-
left = a,
56-
right = b)
57-
ā, b̄, d̄ = _bracket(f, a, b, c)
53+
retcode = ReturnCode.Success,
54+
left = a,
55+
right = b)
56+
ā, b̄, d̄ = _bracket(f, a, b, c)
5857

5958
# The second bracketing block
6059
f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
6160
if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
6261
c = _newton_quadratic(f, ā, b̄, d̄, 3)
63-
else
62+
else
6463
c = _ipzero(f, ā, b̄, d̄, ē)
6564
if (c - ā) * (c - b̄) ≥ 0
6665
c = _newton_quadratic(f, ā, b̄, d̄, 3)
6766
end
6867
end
6968
fc = f(c)
70-
(ā == c || b̄ == c) &&
69+
(ā == c || b̄ == c) &&
7170
return SciMLBase.build_solution(prob, alg, c, fc;
7271
retcode = ReturnCode.FloatingPointLimit,
73-
left = ā,
72+
left = ā,
7473
right = b̄)
7574
iszero(fc) &&
7675
return SciMLBase.build_solution(prob, alg, c, fc;
77-
retcode = ReturnCode.Success,
78-
left = ā,
79-
right = b̄)
80-
ā, b̄, d̄ = _bracket(f, ā, b̄, c)
76+
retcode = ReturnCode.Success,
77+
left = ā,
78+
right = b̄)
79+
ā, b̄, d̄ = _bracket(f, ā, b̄, c)
8180

8281
# The third bracketing block
8382
if abs(f(ā)) < abs(f(b̄))
@@ -90,17 +89,17 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
9089
c = 0.5 * (ā + b̄)
9190
end
9291
fc = f(c)
93-
(ā == c || b̄ == c) &&
92+
(ā == c || b̄ == c) &&
9493
return SciMLBase.build_solution(prob, alg, c, fc;
9594
retcode = ReturnCode.FloatingPointLimit,
96-
left = ā,
95+
left = ā,
9796
right = b̄)
9897
iszero(fc) &&
9998
return SciMLBase.build_solution(prob, alg, c, fc;
100-
retcode = ReturnCode.Success,
101-
left = ā,
102-
right = b̄)
103-
ā, b̄, d = _bracket(f, ā, b̄, c)
99+
retcode = ReturnCode.Success,
100+
left = ā,
101+
right = b̄)
102+
ā, b̄, d = _bracket(f, ā, b̄, c)
104103

105104
# The last bracketing block
106105
if b̄ - ā < 0.5 * (b - a)
@@ -109,14 +108,14 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
109108
e = d
110109
c = 0.5 * (ā + b̄)
111110
fc = f(c)
112-
(ā == c || b̄ == c) &&
111+
(ā == c || b̄ == c) &&
113112
return SciMLBase.build_solution(prob, alg, c, fc;
114113
retcode = ReturnCode.FloatingPointLimit,
115-
left = ā,
114+
left = ā,
116115
right = b̄)
117116
iszero(fc) &&
118117
return SciMLBase.build_solution(prob, alg, c, fc;
119-
retcode = ReturnCode.Success,
118+
retcode = ReturnCode.Success,
120119
left = ā,
121120
right = b̄)
122121
a, b, d = _bracket(f, ā, b̄, c)
@@ -137,43 +136,45 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
137136
end
138137

139138
# Define subrotine function bracket, check fc before bracket to return solution
140-
function _bracket(f::F, a, b, c) where F
139+
function _bracket(f::F, a, b, c) where {F}
141140
if iszero(f(c))
142141
ā, b̄, d = a, b, c
143142
else
144-
if f(a) * f(c) < 0
143+
if f(a) * f(c) < 0
145144
ā, b̄, d = a, c, b
146145
elseif f(b) * f(c) < 0
147146
ā, b̄, d = c, b, a
148147
end
149148
end
150149

151-
return ā, b̄, d
150+
return ā, b̄, d
152151
end
153152

154153
# Define subrotine function newton quadratic, return the approximation of zero
155-
function _newton_quadratic(f::F, a, b, d, k) where F
156-
A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a)
154+
function _newton_quadratic(f::F, a, b, d, k) where {F}
155+
A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a)
157156
B = (f(b) - f(a)) / (b - a)
158157

159158
if iszero(A)
160159
return a - (1 / B) * f(a)
161160
elseif A * f(a) > 0
162-
rᵢ₋₁ = a
163-
else
161+
rᵢ₋₁ = a
162+
else
164163
rᵢ₋₁ = b
165-
end
164+
end
166165

167166
for i in 1:k
168-
rᵢ = rᵢ₋₁ - (f(a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) / (B + A * (2 * rᵢ₋₁ - a - b))
167+
rᵢ = rᵢ₋₁ -
168+
(f(a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) /
169+
(B + A * (2 * rᵢ₋₁ - a - b))
169170
rᵢ₋₁ = rᵢ
170171
end
171172

172173
return rᵢ₋₁
173174
end
174175

175176
# Define subrotine function ipzero, also return the approximation of zero
176-
function _ipzero(f::F, a, b, c, d) where F
177+
function _ipzero(f::F, a, b, c, d) where {F}
177178
Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
178179
Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
179180
Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
@@ -185,4 +186,4 @@ function _ipzero(f::F, a, b, c, d) where F
185186
Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
186187

187188
return a + Q₃₁ + Q₃₂ + Q₃₃
188-
end
189+
end

0 commit comments

Comments
 (0)