@@ -9,18 +9,17 @@ algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal
9
9
struct Alefeld <: AbstractBracketingAlgorithm end
10
10
11
11
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... )
16
15
f = Base. Fix2 (prob. f, prob. p)
17
16
a, b = prob. tspan
18
17
c = a - (b - a) / (f (b) - f (a)) * f (a)
19
-
18
+
20
19
fc = f (c)
21
20
if iszero (fc)
22
21
return SciMLBase. build_solution (prob, alg, c, fc;
23
- retcode = ReturnCode. Success,
22
+ retcode = ReturnCode. Success,
24
23
left = a,
25
24
right = b)
26
25
end
@@ -33,47 +32,47 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
33
32
f₁, f₂, f₃, f₄ = f (a), f (b), f (d), f (e)
34
33
if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
35
34
c = _newton_quadratic (f, a, b, d, 2 )
36
- else
35
+ else
37
36
c = _ipzero (f, a, b, d, e)
38
37
if (c - a) * (c - b) ≥ 0
39
38
c = _newton_quadratic (f, a, b, d, 2 )
40
39
end
41
- end
40
+ end
42
41
ē, fc = d, f (c)
43
- (a == c || b == c) &&
42
+ (a == c || b == c) &&
44
43
return SciMLBase. build_solution (prob, alg, c, fc;
45
44
retcode = ReturnCode. FloatingPointLimit,
46
- left = a,
47
- right = b)
45
+ left = a,
46
+ right = b)
48
47
iszero (fc) &&
49
48
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)
49
+ retcode = ReturnCode. Success,
50
+ left = a,
51
+ right = b)
52
+ ā, b̄, d̄ = _bracket (f, a, b, c)
54
53
55
54
# The second bracketing block
56
55
f₁, f₂, f₃, f₄ = f (ā), f (b̄), f (d̄), f (ē)
57
56
if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
58
57
c = _newton_quadratic (f, ā, b̄, d̄, 3 )
59
- else
58
+ else
60
59
c = _ipzero (f, ā, b̄, d̄, ē)
61
60
if (c - ā) * (c - b̄) ≥ 0
62
61
c = _newton_quadratic (f, ā, b̄, d̄, 3 )
63
62
end
64
63
end
65
64
fc = f (c)
66
- (ā == c || b̄ == c) &&
65
+ (ā == c || b̄ == c) &&
67
66
return SciMLBase. build_solution (prob, alg, c, fc;
68
67
retcode = ReturnCode. FloatingPointLimit,
69
- left = ā,
68
+ left = ā,
70
69
right = b̄)
71
70
iszero (fc) &&
72
71
return SciMLBase. build_solution (prob, alg, c, fc;
73
- retcode = ReturnCode. Success,
74
- left = ā,
75
- right = b̄)
76
- ā, b̄, d̄ = _bracket (f, ā, b̄, c)
72
+ retcode = ReturnCode. Success,
73
+ left = ā,
74
+ right = b̄)
75
+ ā, b̄, d̄ = _bracket (f, ā, b̄, c)
77
76
78
77
# The third bracketing block
79
78
if abs (f (ā)) < abs (f (b̄))
@@ -86,17 +85,17 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
86
85
c = 0.5 * (ā + b̄)
87
86
end
88
87
fc = f (c)
89
- (ā == c || b̄ == c) &&
88
+ (ā == c || b̄ == c) &&
90
89
return SciMLBase. build_solution (prob, alg, c, fc;
91
90
retcode = ReturnCode. FloatingPointLimit,
92
- left = ā,
91
+ left = ā,
93
92
right = b̄)
94
93
iszero (fc) &&
95
94
return SciMLBase. build_solution (prob, alg, c, fc;
96
- retcode = ReturnCode. Success,
97
- left = ā,
98
- right = b̄)
99
- ā, b̄, d = _bracket (f, ā, b̄, c)
95
+ retcode = ReturnCode. Success,
96
+ left = ā,
97
+ right = b̄)
98
+ ā, b̄, d = _bracket (f, ā, b̄, c)
100
99
101
100
# The last bracketing block
102
101
if b̄ - ā < 0.5 * (b - a)
@@ -105,14 +104,14 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
105
104
e = d
106
105
c = 0.5 * (ā + b̄)
107
106
fc = f (c)
108
- (ā == c || b̄ == c) &&
107
+ (ā == c || b̄ == c) &&
109
108
return SciMLBase. build_solution (prob, alg, c, fc;
110
109
retcode = ReturnCode. FloatingPointLimit,
111
- left = ā,
110
+ left = ā,
112
111
right = b̄)
113
112
iszero (fc) &&
114
113
return SciMLBase. build_solution (prob, alg, c, fc;
115
- retcode = ReturnCode. Success,
114
+ retcode = ReturnCode. Success,
116
115
left = ā,
117
116
right = b̄)
118
117
a, b, d = _bracket (f, ā, b̄, c)
@@ -133,43 +132,45 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
133
132
end
134
133
135
134
# Define subrotine function bracket, check fc before bracket to return solution
136
- function _bracket (f:: F , a, b, c) where F
135
+ function _bracket (f:: F , a, b, c) where {F}
137
136
if iszero (f (c))
138
137
ā, b̄, d = a, b, c
139
138
else
140
- if f (a) * f (c) < 0
139
+ if f (a) * f (c) < 0
141
140
ā, b̄, d = a, c, b
142
141
elseif f (b) * f (c) < 0
143
142
ā, b̄, d = c, b, a
144
143
end
145
144
end
146
145
147
- return ā, b̄, d
146
+ return ā, b̄, d
148
147
end
149
148
150
149
# 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)
150
+ function _newton_quadratic (f:: F , a, b, d, k) where {F}
151
+ A = ((f (d) - f (b)) / (d - b) - (f (b) - f (a)) / (b - a)) / (d - a)
153
152
B = (f (b) - f (a)) / (b - a)
154
153
155
154
if iszero (A)
156
155
return a - (1 / B) * f (a)
157
156
elseif A * f (a) > 0
158
- rᵢ₋₁ = a
159
- else
157
+ rᵢ₋₁ = a
158
+ else
160
159
rᵢ₋₁ = b
161
- end
160
+ end
162
161
163
162
for i in 1 : k
164
- rᵢ = rᵢ₋₁ - (f (a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) / (B + A * (2 * rᵢ₋₁ - a - b))
163
+ rᵢ = rᵢ₋₁ -
164
+ (f (a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) /
165
+ (B + A * (2 * rᵢ₋₁ - a - b))
165
166
rᵢ₋₁ = rᵢ
166
167
end
167
168
168
169
return rᵢ₋₁
169
170
end
170
171
171
172
# Define subrotine function ipzero, also return the approximation of zero
172
- function _ipzero (f:: F , a, b, c, d) where F
173
+ function _ipzero (f:: F , a, b, c, d) where {F}
173
174
Q₁₁ = (c - d) * f (c) / (f (d) - f (c))
174
175
Q₂₁ = (b - c) * f (b) / (f (c) - f (b))
175
176
Q₃₁ = (a - b) * f (a) / (f (b) - f (a))
@@ -181,4 +182,4 @@ function _ipzero(f::F, a, b, c, d) where F
181
182
Q₃₃ = (D₃₂ - Q₂₂) * f (a) / (f (d) - f (a))
182
183
183
184
return a + Q₃₁ + Q₃₂ + Q₃₃
184
- end
185
+ end
0 commit comments