]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /* Boost example/findroot_demo.cpp |
2 | * find zero points of some function by dichotomy | |
3 | * | |
4 | * Copyright 2000 Jens Maurer | |
5 | * Copyright 2002-2003 Guillaume Melquiond | |
6 | * | |
7 | * Distributed under the Boost Software License, Version 1.0. | |
8 | * (See accompanying file LICENSE_1_0.txt or | |
9 | * copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | * | |
11 | * The idea and the 2D function are based on RVInterval, | |
12 | * which contains the following copyright notice: | |
13 | ||
14 | This file is copyrighted 1996 by Ronald Van Iwaarden. | |
15 | ||
16 | Permission is hereby granted, without written agreement and | |
17 | without license or royalty fees, to use, copy, modify, and | |
18 | distribute this software and its documentation for any | |
19 | purpose, subject to the following conditions: | |
20 | ||
21 | The above license notice and this permission notice shall | |
22 | appear in all copies or substantial portions of this software. | |
23 | ||
24 | The name "RVInterval" cannot be used for any modified form of | |
25 | this software that does not originate from the authors. | |
26 | Nevertheless, the name "RVInterval" may and should be used to | |
27 | designate the optimization software implemented and described | |
28 | in this package, even if embedded in any other system, as long | |
29 | as any portion of this code remains. | |
30 | ||
31 | The authors specifically disclaim any warranties, including, | |
32 | but not limited to, the implied warranties of merchantability | |
33 | and fitness for a particular purpose. The software provided | |
34 | hereunder is on an "as is" basis, and the authors have no | |
35 | obligation to provide maintenance, support, updates, | |
36 | enhancements, or modifications. In no event shall the authors | |
37 | be liable to any party for direct, indirect, special, | |
38 | incidental, or consequential damages arising out of the use of | |
39 | this software and its documentation. | |
40 | */ | |
41 | ||
42 | #include <boost/numeric/interval.hpp> // must be first for <limits> workaround | |
43 | #include <list> | |
44 | #include <deque> | |
45 | #include <vector> | |
46 | #include <fstream> | |
47 | #include <iostream> | |
48 | ||
49 | ||
50 | template<class T> | |
51 | struct test_func2d | |
52 | { | |
53 | T operator()(T x, T y) const | |
54 | { | |
55 | return sin(x)*cos(y) - exp(x*y)/45.0 * (pow(x+y, 2)+100.0) - | |
56 | cos(sin(y))*y/4.0; | |
57 | } | |
58 | }; | |
59 | ||
60 | template <class T> | |
61 | struct test_func1d | |
62 | { | |
63 | T operator()(T x) const | |
64 | { | |
65 | return sin(x)/(x*x+1.0); | |
66 | } | |
67 | }; | |
68 | ||
69 | template<class T> | |
70 | struct test_func1d_2 | |
71 | { | |
72 | T operator()(T x) const | |
73 | { | |
74 | using std::sqrt; | |
75 | return sqrt(x*x-1.0); | |
76 | } | |
77 | }; | |
78 | ||
79 | template<class Function, class I> | |
80 | void find_zeros(std::ostream & os, Function f, I searchrange) | |
81 | { | |
82 | std::list<I> l, done; | |
83 | l.push_back(searchrange); | |
84 | while(!l.empty()) { | |
85 | I range = l.front(); | |
86 | l.pop_front(); | |
87 | I val = f(range); | |
88 | if (zero_in(val)) { | |
89 | if(width(range) < 1e-6) { | |
90 | os << range << '\n'; | |
91 | continue; | |
92 | } | |
93 | // there's still a solution hidden somewhere | |
94 | std::pair<I,I> p = bisect(range); | |
95 | l.push_back(p.first); | |
96 | l.push_back(p.second); | |
97 | } | |
98 | } | |
99 | } | |
100 | ||
101 | template<class T> | |
102 | std::ostream &operator<<(std::ostream &os, const std::pair<T, T> &x) { | |
103 | os << "(" << x.first << ", " << x.second << ")"; | |
104 | return os; | |
105 | } | |
106 | ||
107 | template<class T, class Policies> | |
108 | std::ostream &operator<<(std::ostream &os, | |
109 | const boost::numeric::interval<T, Policies> &x) { | |
110 | os << "[" << x.lower() << ", " << x.upper() << "]"; | |
111 | return os; | |
112 | } | |
113 | ||
114 | static const double epsilon = 5e-3; | |
115 | ||
116 | template<class Function, class I> | |
117 | void find_zeros(std::ostream & os, Function f, I rx, I ry) | |
118 | { | |
119 | typedef std::pair<I, I> rectangle; | |
120 | typedef std::deque<rectangle> container; | |
121 | container l, done; | |
122 | // l.reserve(50); | |
123 | l.push_back(std::make_pair(rx, ry)); | |
124 | for(int i = 1; !l.empty(); ++i) { | |
125 | rectangle rect = l.front(); | |
126 | l.pop_front(); | |
127 | I val = f(rect.first, rect.second); | |
128 | if (zero_in(val)) { | |
129 | if(width(rect.first) < epsilon && width(rect.second) < epsilon) { | |
130 | os << median(rect.first) << " " << median(rect.second) << " " | |
131 | << lower(rect.first) << " " << upper(rect.first) << " " | |
132 | << lower(rect.second) << " " << upper(rect.second) | |
133 | << '\n'; | |
134 | } else { | |
135 | if(width(rect.first) > width(rect.second)) { | |
136 | std::pair<I,I> p = bisect(rect.first); | |
137 | l.push_back(std::make_pair(p.first, rect.second)); | |
138 | l.push_back(std::make_pair(p.second, rect.second)); | |
139 | } else { | |
140 | std::pair<I,I> p = bisect(rect.second); | |
141 | l.push_back(std::make_pair(rect.first, p.first)); | |
142 | l.push_back(std::make_pair(rect.first, p.second)); | |
143 | } | |
144 | } | |
145 | } | |
146 | if(i % 10000 == 0) | |
147 | std::cerr << "\rIteration " << i << ", l.size() = " << l.size(); | |
148 | } | |
149 | std::cerr << '\n'; | |
150 | } | |
151 | ||
152 | int main() | |
153 | { | |
154 | using namespace boost; | |
155 | using namespace numeric; | |
156 | using namespace interval_lib; | |
157 | ||
158 | typedef interval<double, | |
159 | policies<save_state<rounded_transc_opp<double> >, | |
160 | checking_base<double> > > I; | |
161 | ||
162 | std::cout << "Zero points of sin(x)/(x*x+1)\n"; | |
163 | find_zeros(std::cout, test_func1d<I>(), I(-11, 10)); | |
164 | std::cout << "Zero points of sqrt(x*x-1)\n"; | |
165 | find_zeros(std::cout, test_func1d_2<I>(), I(-5, 6)); | |
166 | std::cout << "Zero points of Van Iwaarden's 2D function\n"; | |
167 | std::ofstream f("func2d.data"); | |
168 | find_zeros(f, test_func2d<I>(), I(-20, 20), I(-20, 20)); | |
169 | std::cout << "Use gnuplot, command 'plot \"func2d.data\" with dots' to plot\n"; | |
170 | } |