123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171 |
- /* Boost example/findroot_demo.cpp
- * find zero points of some function by dichotomy
- *
- * Copyright 2000 Jens Maurer
- * Copyright 2002-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)
- *
- * The idea and the 2D function are based on RVInterval,
- * which contains the following copyright notice:
- This file is copyrighted 1996 by Ronald Van Iwaarden.
- Permission is hereby granted, without written agreement and
- without license or royalty fees, to use, copy, modify, and
- distribute this software and its documentation for any
- purpose, subject to the following conditions:
-
- The above license notice and this permission notice shall
- appear in all copies or substantial portions of this software.
-
- The name "RVInterval" cannot be used for any modified form of
- this software that does not originate from the authors.
- Nevertheless, the name "RVInterval" may and should be used to
- designate the optimization software implemented and described
- in this package, even if embedded in any other system, as long
- as any portion of this code remains.
-
- The authors specifically disclaim any warranties, including,
- but not limited to, the implied warranties of merchantability
- and fitness for a particular purpose. The software provided
- hereunder is on an "as is" basis, and the authors have no
- obligation to provide maintenance, support, updates,
- enhancements, or modifications. In no event shall the authors
- be liable to any party for direct, indirect, special,
- incidental, or consequential damages arising out of the use of
- this software and its documentation.
- */
- #include <boost/numeric/interval.hpp> // must be first for <limits> workaround
- #include <boost/numeric/interval/io.hpp>
- #include <list>
- #include <deque>
- #include <vector>
- #include <fstream>
- #include <iostream>
- template<class T>
- struct test_func2d
- {
- T operator()(T x, T y) const
- {
- return sin(x)*cos(y) - exp(x*y)/45.0 * (pow(x+y, 2)+100.0) -
- cos(sin(y))*y/4.0;
- }
- };
- template <class T>
- struct test_func1d
- {
- T operator()(T x) const
- {
- return sin(x)/(x*x+1.0);
- }
- };
- template<class T>
- struct test_func1d_2
- {
- T operator()(T x) const
- {
- using std::sqrt;
- return sqrt(x*x-1.0);
- }
- };
- template<class Function, class I>
- void find_zeros(std::ostream & os, Function f, I searchrange)
- {
- std::list<I> l, done;
- l.push_back(searchrange);
- while(!l.empty()) {
- I range = l.front();
- l.pop_front();
- I val = f(range);
- if (zero_in(val)) {
- if(width(range) < 1e-6) {
- os << range << '\n';
- continue;
- }
- // there's still a solution hidden somewhere
- std::pair<I,I> p = bisect(range);
- l.push_back(p.first);
- l.push_back(p.second);
- }
- }
- }
- template<class T>
- std::ostream &operator<<(std::ostream &os, const std::pair<T, T> &x) {
- os << "(" << x.first << ", " << x.second << ")";
- return os;
- }
- 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;
- }
- static const double epsilon = 5e-3;
- template<class Function, class I>
- void find_zeros(std::ostream & os, Function f, I rx, I ry)
- {
- typedef std::pair<I, I> rectangle;
- typedef std::deque<rectangle> container;
- container l, done;
- // l.reserve(50);
- l.push_back(std::make_pair(rx, ry));
- for(int i = 1; !l.empty(); ++i) {
- rectangle rect = l.front();
- l.pop_front();
- I val = f(rect.first, rect.second);
- if (zero_in(val)) {
- if(width(rect.first) < epsilon && width(rect.second) < epsilon) {
- os << median(rect.first) << " " << median(rect.second) << " "
- << lower(rect.first) << " " << upper(rect.first) << " "
- << lower(rect.second) << " " << upper(rect.second)
- << '\n';
- } else {
- if(width(rect.first) > width(rect.second)) {
- std::pair<I,I> p = bisect(rect.first);
- l.push_back(std::make_pair(p.first, rect.second));
- l.push_back(std::make_pair(p.second, rect.second));
- } else {
- std::pair<I,I> p = bisect(rect.second);
- l.push_back(std::make_pair(rect.first, p.first));
- l.push_back(std::make_pair(rect.first, p.second));
- }
- }
- }
- if(i % 10000 == 0)
- std::cerr << "\rIteration " << i << ", l.size() = " << l.size();
- }
- std::cerr << '\n';
- }
- int main()
- {
- using namespace boost;
- using namespace numeric;
- using namespace interval_lib;
- typedef interval<double,
- policies<save_state<rounded_transc_opp<double> >,
- checking_base<double> > > I;
- std::cout << "Zero points of sin(x)/(x*x+1)\n";
- find_zeros(std::cout, test_func1d<I>(), I(-11, 10));
- std::cout << "Zero points of sqrt(x*x-1)\n";
- find_zeros(std::cout, test_func1d_2<I>(), I(-5, 6));
- std::cout << "Zero points of Van Iwaarden's 2D function\n";
- std::ofstream f("func2d.data");
- find_zeros(f, test_func2d<I>(), I(-20, 20), I(-20, 20));
- std::cout << "Use gnuplot, command 'plot \"func2d.data\" with dots' to plot\n";
- }
|