more/num/nonlinear.h Source File

00001 
00002 //  Copyright (C) 1999--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)
00003 
00004 //  This file is free software; you can redistribute it and/or modify
00005 //  it under the terms of the GNU General Public License as published by
00006 //  the Free Software Foundation; either version 2 of the License, or
00007 //  (at your option) any later version.
00008 
00009 //  This file is distributed in the hope that it will be useful,
00010 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 //  GNU General Public License for more details.
00013 
00014 //  You should have received a copy of the GNU General Public License
00015 //  along with this program; if not, write to the Free Software
00016 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018 //  As a special exception, you may use this file as part of a free
00019 //  software library without restriction.  Specifically, if other files
00020 //  instantiate templates or use macros or inline functions from this
00021 //  file, or you compile this file and link it with other files to
00022 //  produce an executable, this file does not by itself cause the
00023 //  resulting executable to be covered by the GNU General Public
00024 //  License.  This exception does not however invalidate any other
00025 //  reasons why the executable file might be covered by the GNU General
00026 //  Public License.
00027 
00028 //  $Id: nonlinear.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $
00029 
00030 
00031 #ifndef MORE_NONLINEAR_H
00032 #define MORE_NONLINEAR_H
00033 
00034 #include <utility>
00035 #include <stdexcept>
00036 #include <cmath>
00037 
00038 namespace more {
00039 namespace num {
00040 
00041   template<typename T, typename Function>
00042     std::pair<T, bool>
00043     bisection_root(Function f, T xl, T xh, T precx, T precf, int maxit=128) {
00044         maxit /= 8;
00045         const bool signl = f(xl) > 0;
00046         T fm, xm;
00047         if(signl == (f(xh) > 0))
00048             throw std::logic_error("bisection_roots: f(xl)*f(xh) > 0");
00049         do {
00050             xm = .5*(xl+xh);
00051             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00052             xm = .5*(xl+xh);
00053             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00054             xm = .5*(xl+xh);
00055             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00056             xm = .5*(xl+xh);
00057             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00058             xm = .5*(xl+xh);
00059             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00060             xm = .5*(xl+xh);
00061             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00062             xm = .5*(xl+xh);
00063             if((f(xm) > 0) == signl) xl = xm; else xh = xm;
00064             xm = .5*(xl+xh);
00065             fm = f(xm);
00066             if(xm == xh || xm == xl) break;
00067             if(!--maxit) return std::pair<T, bool>(xm, false);
00068             if((fm > 0) == signl) xl = xm; else xh = xm;
00069         } while(std::fabs(fm) > precf || std::fabs(xh-xl) > precx);
00070         return std::pair<T, bool>(xm, true);
00071     }
00072 
00073   template<typename T, typename Function, typename OutputIterator>
00074     std::pair<OutputIterator, bool>
00075     bisection_roots(Function f, T xl, T xh, T dx, int n,
00076                     T precx, T precf, OutputIterator it_out) {
00077         bool isconv = true;
00078         bool signdx = dx > 0;
00079         bool signl = f(xl) > 0;
00080         while((xl < xh) == signdx) {
00081             xl += dx;
00082             bool signm = f(xl) > 0;
00083             if(signl != signm) {
00084                 std::pair<T, bool>
00085                     root = bisection_root(f, xl-dx, xl, precx, precf);
00086                 isconv = isconv && root.second;
00087                 (*it_out) = root.first;
00088                 ++it_out;
00089                 signl = signm;
00090             }
00091         }
00092         return std::pair<OutputIterator, bool>(it_out, isconv);
00093     }
00094 
00095   template<typename T, typename Function, typename OutputIterator>
00096     std::pair<OutputIterator, bool>
00097     bisection_roots(Function f, T xl, T xh, T dx, int n,
00098                     OutputIterator it_out) {
00099         return bisection_roots(f, xl, xh, dx, n, 0.0, 1.0, it_out);
00100     }
00101 }}
00102 
00103 #endif