• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 <boost/numeric/interval/io.hpp>
44 #include <list>
45 #include <deque>
46 #include <vector>
47 #include <fstream>
48 #include <iostream>
49 
50 
51 template<class T>
52 struct test_func2d
53 {
operator ()test_func2d54   T operator()(T x, T y) const
55   {
56     return sin(x)*cos(y) - exp(x*y)/45.0 * (pow(x+y, 2)+100.0) -
57       cos(sin(y))*y/4.0;
58   }
59 };
60 
61 template <class T>
62 struct test_func1d
63 {
operator ()test_func1d64   T operator()(T x) const
65   {
66     return sin(x)/(x*x+1.0);
67   }
68 };
69 
70 template<class T>
71 struct test_func1d_2
72 {
operator ()test_func1d_273   T operator()(T x) const
74   {
75     using std::sqrt;
76     return sqrt(x*x-1.0);
77   }
78 };
79 
80 template<class Function, class I>
find_zeros(std::ostream & os,Function f,I searchrange)81 void find_zeros(std::ostream & os, Function f, I searchrange)
82 {
83   std::list<I> l, done;
84   l.push_back(searchrange);
85   while(!l.empty()) {
86     I range = l.front();
87     l.pop_front();
88     I val = f(range);
89     if (zero_in(val)) {
90       if(width(range) < 1e-6) {
91         os << range << '\n';
92         continue;
93       }
94       // there's still a solution hidden somewhere
95       std::pair<I,I> p = bisect(range);
96       l.push_back(p.first);
97       l.push_back(p.second);
98     }
99   }
100 }
101 
102 template<class T>
operator <<(std::ostream & os,const std::pair<T,T> & x)103 std::ostream &operator<<(std::ostream &os, const std::pair<T, T> &x) {
104   os << "(" << x.first << ", " << x.second << ")";
105   return os;
106 }
107 
108 template<class T, class Policies>
operator <<(std::ostream & os,const boost::numeric::interval<T,Policies> & x)109 std::ostream &operator<<(std::ostream &os,
110                          const boost::numeric::interval<T, Policies> &x) {
111   os << "[" << x.lower() << ", " << x.upper() << "]";
112   return os;
113 }
114 
115 static const double epsilon = 5e-3;
116 
117 template<class Function, class I>
find_zeros(std::ostream & os,Function f,I rx,I ry)118 void find_zeros(std::ostream & os, Function f, I rx, I ry)
119 {
120   typedef std::pair<I, I> rectangle;
121   typedef std::deque<rectangle> container;
122   container l, done;
123   // l.reserve(50);
124   l.push_back(std::make_pair(rx, ry));
125   for(int i = 1; !l.empty(); ++i) {
126     rectangle rect = l.front();
127     l.pop_front();
128     I val = f(rect.first, rect.second);
129     if (zero_in(val)) {
130       if(width(rect.first) < epsilon && width(rect.second) < epsilon) {
131         os << median(rect.first) << " " << median(rect.second) << " "
132            << lower(rect.first) << " " << upper(rect.first) << " "
133            << lower(rect.second) << " " << upper(rect.second)
134            << '\n';
135       } else {
136         if(width(rect.first) > width(rect.second)) {
137           std::pair<I,I> p = bisect(rect.first);
138           l.push_back(std::make_pair(p.first, rect.second));
139           l.push_back(std::make_pair(p.second, rect.second));
140         } else {
141           std::pair<I,I> p = bisect(rect.second);
142           l.push_back(std::make_pair(rect.first, p.first));
143           l.push_back(std::make_pair(rect.first, p.second));
144         }
145       }
146     }
147     if(i % 10000 == 0)
148       std::cerr << "\rIteration " << i << ", l.size() = " << l.size();
149   }
150   std::cerr << '\n';
151 }
152 
main()153 int main()
154 {
155   using namespace boost;
156   using namespace numeric;
157   using namespace interval_lib;
158 
159   typedef interval<double,
160                    policies<save_state<rounded_transc_opp<double> >,
161                             checking_base<double> > > I;
162 
163   std::cout << "Zero points of sin(x)/(x*x+1)\n";
164   find_zeros(std::cout, test_func1d<I>(), I(-11, 10));
165   std::cout << "Zero points of sqrt(x*x-1)\n";
166   find_zeros(std::cout, test_func1d_2<I>(), I(-5, 6));
167   std::cout << "Zero points of Van Iwaarden's 2D function\n";
168   std::ofstream f("func2d.data");
169   find_zeros(f, test_func2d<I>(), I(-20, 20), I(-20, 20));
170   std::cout << "Use gnuplot, command 'plot \"func2d.data\" with dots'   to plot\n";
171 }
172