1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
/* Boost example/newton-raphson.cpp
* Newton iteration for intervals
*
* Copyright 2003 Guillaume Melquiond
*
* Distributed under the Boost Software License, Version 1.0.
* (See accompanying file LICENSE_1_0.txt or
* copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <boost/numeric/interval.hpp>
#include <vector>
#include <algorithm>
#include <utility>
#include <iostream>
#include <iomanip>
template <class I> I f(const I& x)
{ return x * (x - 1.) * (x - 2.) * (x - 3.) * (x - 4.); }
template <class I> I f_diff(const I& x)
{ return (((5. * x - 40.) * x + 105.) * x - 100.) * x + 24.; }
static const double max_width = 1e-10;
static const double alpha = 0.75;
using namespace boost;
using namespace numeric;
using namespace interval_lib;
// First method: no empty intervals
typedef interval<double> I1_aux;
typedef unprotect<I1_aux>::type I1;
std::vector<I1> newton_raphson(const I1& xs) {
std::vector<I1> l, res;
I1 vf, vd, x, x1, x2;
l.push_back(xs);
while (!l.empty()) {
x = l.back();
l.pop_back();
bool x2_used;
double xx = median(x);
vf = f<I1>(xx);
vd = f_diff<I1>(x);
if (zero_in(vf) && zero_in(vd)) {
x1 = I1::whole();
x2_used = false;
} else {
x1 = xx - division_part1(vf, vd, x2_used);
if (x2_used) x2 = xx - division_part2(vf, vd);
}
if (overlap(x1, x)) x1 = intersect(x, x1);
else if (x2_used) { x1 = x2; x2_used = false; }
else continue;
if (x2_used)
if (overlap(x2, x)) x2 = intersect(x, x2);
else x2_used = false;
if (x2_used && width(x2) > width(x1)) std::swap(x1, x2);
if (!zero_in(f(x1)))
if (x2_used) { x1 = x2; x2_used = false; }
else continue;
if (width(x1) < max_width) res.push_back(x1);
else if (width(x1) > alpha * width(x)) {
std::pair<I1, I1> p = bisect(x);
if (zero_in(f(p.first))) l.push_back(p.first);
x2 = p.second;
x2_used = true;
} else l.push_back(x1);
if (x2_used && zero_in(f(x2)))
if (width(x2) < max_width) res.push_back(x2);
else l.push_back(x2);
}
return res;
}
// Second method: with empty intervals
typedef change_checking<I1_aux, checking_no_nan<double> >::type I2_aux;
typedef unprotect<I2_aux>::type I2;
std::vector<I2> newton_raphson(const I2& xs) {
std::vector<I2> l, res;
I2 vf, vd, x, x1, x2;
l.push_back(xs);
while (!l.empty()) {
x = l.back();
l.pop_back();
double xx = median(x);
vf = f<I2>(xx);
vd = f_diff<I2>(x);
if (zero_in(vf) && zero_in(vd)) {
x1 = x;
x2 = I2::empty();
} else {
bool x2_used;
x1 = intersect(x, xx - division_part1(vf, vd, x2_used));
x2 = intersect(x, xx - division_part2(vf, vd, x2_used));
}
if (width(x2) > width(x1)) std::swap(x1, x2);
if (empty(x1) || !zero_in(f(x1)))
if (!empty(x2)) { x1 = x2; x2 = I2::empty(); }
else continue;
if (width(x1) < max_width) res.push_back(x1);
else if (width(x1) > alpha * width(x)) {
std::pair<I2, I2> p = bisect(x);
if (zero_in(f(p.first))) l.push_back(p.first);
x2 = p.second;
} else l.push_back(x1);
if (!empty(x2) && zero_in(f(x2)))
if (width(x2) < max_width) res.push_back(x2);
else l.push_back(x2);
}
return res;
}
template<class T, class Policies>
std::ostream &operator<<(std::ostream &os,
const boost::numeric::interval<T, Policies> &x) {
os << "[" << x.lower() << ", " << x.upper() << "]";
return os;
}
int main() {
{
I1_aux::traits_type::rounding rnd;
std::vector<I1> res = newton_raphson(I1(-1, 5.1));
std::cout << "Results: " << std::endl << std::setprecision(12);
for(std::vector<I1>::const_iterator i = res.begin(); i != res.end(); ++i)
std::cout << " " << *i << std::endl;
std::cout << std::endl;
}
{
I2_aux::traits_type::rounding rnd;
std::vector<I2> res = newton_raphson(I2(-1, 5.1));
std::cout << "Results: " << std::endl << std::setprecision(12);
for(std::vector<I2>::const_iterator i = res.begin(); i != res.end(); ++i)
std::cout << " " << *i << std::endl;
std::cout << std::endl;
}
}
|
