[Boost-users] RE Is there a false-position root finding method in boost

Subject: [Boost-users] RE Is there a false-position root finding method in boost
From: frederic.bron_at_[hidden]
Date: 2008-09-11 02:59:09

Look here: http://www.boost.org/doc/libs/1_36_0/libs/math/doc/sf_and_dist/html/math_toolkit/toolkit/internals1/roots2.html

I have also written the following inspire by Numerical Recipes:

template <class F, class T>
T false_position_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t &max_iter) {
        T fmin=f(min), fmax=f(max) ;

        if (not is_of_opposite_sign(fmin, fmax)){ // not bracketed
                max_iter=0 ;
                return (std::abs(fmin)>std::abs(fmax)) ? max : min ;
        }
        if (fmax<0) { // have fmin<fmax even if min becomes > max
                std::swap(min, max) ;
                std::swap(fmin, fmax) ;
        }

        T fmed, diff, med=guess ;
        T factor=static_cast<T>(ldexp(1., 1-digits)) ;
        boost::uintmax_t count=max_iter ;
        do {
                fmed=f(med) ;
                if (fmed==0) { --count ; break ; }
                if (fmed<0) {
                        diff=std::abs(med-min) ;
                        min=med ;
                        fmin=fmed ;
                } else {
                        diff=std::abs(max-med) ;
                        max=med ;
                        fmax=fmed ;
                }
                med=min-(fmin/(fmax-fmin))*(max-min) ;
        } while(--count and diff>std::abs(med*factor)) ;
        max_iter-=count ;
        return med ;
}

template <class F, class T>
T false_position_iterate(F f, T guess, T min, T max, int digits) {
        boost::uintmax_t m=std::numeric_limits<boost::uintmax_t>::max() ;
        return false_position_iterate(f, guess, min, max, digits, m) ;
}

template <class F, class T>
T ridders_iterate(F f, T guess, T x1, T x2, int digits, boost::uintmax_t &max_iter) {
        T f1=f(x1), f2=f(x2), fm, fn ;
        T xm, xn, s, diff, root=std::numeric_limits<T>::max() ;

        if (not is_of_opposite_sign(f1, f2)) { // not bracketed
                max_iter=0 ;
                return (std::abs(f1)>std::abs(f2)) ? x2 : x1 ;
        }

        if (f1==0) {
                max_iter=0 ;
                return x1 ;
        } else if (f2==0) {
                max_iter=0 ;
                return x2 ;
        }

        T factor=static_cast<T>(ldexp(1., 1-digits)) ;
        boost::uintmax_t count=max_iter ;
        do {
                fm=f(xm=0.5f*(x1+x2)) ;
                s=std::sqrt(fm*fm-f1*f2) ;
                if (not s) { root=xm ; --count ; break ; }
                xn=xm ;
                if (f1>=f2) xn+=(xm-x1)*fm/s ; else xn-=(xm-x1)*fm/s ;
                diff=std::abs(xn-root) ;
                if (diff<=std::abs(root*factor)) { --count ; break ; }
                fn=f(root=xn) ;
                if (fn==0) { --count ; break ; }
                if (is_of_opposite_sign(fn, fm)) {
                        x1=xm ;
                        f1=fm ;
                        x2=xn ;
                        f2=fn ;
                } else if (is_of_opposite_sign(fn, f1)) {
                        x2=xn ;
                        f2=fn ;
                } else {
                        x1=xn ;
                        f1=fn ;
                }
                diff=std::min(diff, std::abs(x2-x1)) ;
        } while (--count and diff>std::abs(root*factor)) ;
        max_iter-=count ;
        return root ;
}

template <class F, class T>
T ridders_iterate(F f, T guess, T min, T max, int digits){
        boost::uintmax_t m=std::numeric_limits<boost::uintmax_t>::max();
        return ridders_iterate(f, guess, min, max, digits, m) ;
}

F. Bron

boost-users-bounces_at_[hidden] a écrit sur 11/09/2008 05:57:02 :

> #include <boost/math/tools/roots.hpp>
> There are a few root finding algorithm in the above file. I'm not sure
> if there is the false-position method in there?
>
> What I want is give an interval which only has one root. I want to
> always find the root in that interval. I'm wondering if any function
> in that header can do so or not.

Avis :
Ce message et toute pièce jointe sont la propriété d'Alcan et sont destinés seulement aux personnes ou à l'entité à qui le message est adressé. Si vous avez reçu ce message par erreur, veuillez le
détruire et en aviser l'expéditeur par courriel. Si vous n'êtes pas le destinataire du message, vous n'êtes pas autorisé à utiliser, à copier ou à divulguer le contenu du message ou ses pièces jointes
en tout ou en partie.

Notice:
This message and any attachments are the property of Alcan and are intended solely for the named recipients or entity to whom this message is addressed. If you have received this message in error
please inform the sender via e-mail and destroy the message. If you are not the intended recipient you are not allowed to use, copy or disclose the contents or attachments in whole or in part.