ot_root_find.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       ot_root_find.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:53:59 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include <stdio.h>
#include <math.h>
 
#pragma hdrstop
 
#include "ot_root_find.h"
 
#pragma package(smart_init)
 
#define SMALL_NUM_1 1E-14
#define SMALL_NUM_2 1E-28
 
unsigned OT_ROOT_FIND_1D::ZeroBracOut (Function fn, double &x1, double &x2, void* pvData) {
    const double FACTOR = 1.6;
    const unsigned NTRY = 50;
 
    double f1 = fn (x1, pvData), f2 = fn (x2, pvData);
    unsigned j;
 
    if (x1 == x2)
    {
        //    printf ("you have to guess an initial range");
        return 1e8;
    }
    for (j = 1; j <= NTRY; j++) {
        if (f1*f2 < 0.) return 1;           // found a root
        if (fabs (f1) < fabs (f2)) {
            x1 += FACTOR*(x1-x2);             // expand towards x1
            f1 = fn (x1, pvData);
        }
        else {
            x2 += FACTOR*(x2-x1);             // ... or towards x2
            f2 = fn (x2, pvData);
        }
    }
    return 1E8;
}                               // no root bracketed
 
unsigned OT_ROOT_FIND_1D::ZeroBracOut (Function fn, double xStart, double xEnd, double &x1, double &x2, void* pvData) {
    const double FACTOR = 1.6;
    const unsigned NTRY = 50;
 
    if (x1>x2) {
        double tmp=x1;
        x1=x2;
        x2=tmp;
    }
    if (xStart>xEnd) {
        double tmp=xStart;
        xEnd=xStart;
        xEnd=tmp;
    }
    double deltaT=xEnd-xStart;
 
    if (x1 < xStart)
        x1 = xStart;
    else if (x1 > xEnd)
        x1 = xEnd;
    if (x2 < xStart)
        x2 = xStart;
    else if (x2 > xEnd)
        x2 = xEnd;
 
    if (x1 == x2)
    {
        x1 -= .1*deltaT;
        x1 = (x1>xStart)?x1:xStart;
        x2 += .1*deltaT;
        x2 = (x2<xEnd)?x2:xEnd;
    }
 
    double f1 = fn (x1, pvData), f2 = fn (x2, pvData);
    unsigned j;
 
    if (x1 == x2)
    {
        //    printf ("you have to guess an initial range");
        return 1E8;
    }
    for (j = 1; j <= NTRY; j++) {
        if (f1*f2 < 0.) return 1;           // found a root
        if (fabs (f1) < fabs (f2) && x1 > xStart )
        {
            if (x1+FACTOR*(x1-x2) > xStart)
            {
                x1 += FACTOR*(x1-x2);             // expand towards x1
            }
            else
            {
                x1 = xStart;
            }
            f1 = fn (x1, pvData);
        }
        else if (x2 < xEnd ) {
            if (x2+FACTOR*(x2-x1) < xEnd)
            {
                x2 += FACTOR*(x2-x1);             // ... or towards x2
            }
            else
            {
                x2 = xEnd;
            }
            f2 = fn (x2, pvData);
        }
        else
            break;
    }
    return 0;
}                               // no root bracketed
 
unsigned OT_ROOT_FIND_1D::ZeroBracIn (Function fn, const double &x1, const double &x2,
                                      unsigned n, double xb1[], double xb2[], unsigned nb, void* pvData) {
    unsigned i, nbb = nb;
    double x = x1, dx = (x2-x1)/n,      // determine appropriate spacing
                        fp = fn (x, pvData);
 
    nb = 0;
    for (i = 1; i <= n; i++) {          // loop over all intervals
        double fc;
 
        x += dx;
        fc = fn (x, pvData);
        if (fc*fp < 0.) {                       // if a sign change occurs, record
            xb1[nb] = x-dx;                     // values for the bounds.
            xb2[nb] = x;
            nb++;
        }
        fp = fc;
        if (nbb == nb) return nb;
    }
    return nb;
}
 
 
double OT_ROOT_FIND_1D::RootBisect (Function fn, const double &x1, const double &x2, const double &xacc, void* pvData)
{
    const unsigned JMAX = 50;
    double fmid = fn (x2, pvData), f = fn (x1, pvData), rtbis, dx;
    unsigned j;
 
    if (fmid*f >= 0)
    {
        //    printf ("root must be bracketed for bisection");
        return 1E308;
    }
 
    if (f < 0.) {                               // orient the search so that f > 0
        rtbis = x1;                             // lies at x+dx.
        dx = x2-x1;
    }
    else {
        rtbis = x2;
        dx = x1-x2;
    }
 
    for (j = 1; j <= JMAX; j++) {               // bisection loop
        double xmid;
 
        dx *= .5;                           // here we halve the interval
        xmid = rtbis+dx;                    // and update boundaries and
        fmid = fn (xmid, pvData);                   // function value
        if (fmid <= 0.) rtbis = xmid;
        if ((fabs (dx) <= xacc)||(fmid == 0.)) return rtbis;
    }
    //  printf ("too many bisections");
    return 1E308;
}       // this should happen rarely (2^-50 ~ 10^-15)
 
 
double OT_ROOT_FIND_1D::RootFlsPos (Function fn, const double &x1, const double &x2,
                                    const double &xacc, void* pvData) {
    const unsigned MAXIT = 40;
    double flo = fn (x1, pvData), fhi = fn (x2, pvData), xl, xh, dx, dummy;
    unsigned j;
 
    if (flo*fhi >= 0)
    {
        //    printf ("root must be bracketed for false position");
        return 1E308;
    }
 
    if (flo < 0.) {                     // let xl correspond to the low side
        xl = x1;
        xh = x2;
    }
    else {
        {
            dummy=flo;
            flo=fhi;
            fhi=dummy;
        }
        xl = x2;
        xh = x1;
    }
    dx = xh-xl;
 
    for (j = 1; j <= MAXIT; j++) {      // false position loop
        double rtflsp = xl+dx*flo/(flo-fhi),// increment w/ respect to latest
                        f = fn (rtflsp, pvData),                // value
                            del = (f < 0.)? xl-rtflsp: xh-rtflsp;
 
        if (f < 0.) {                   // replace appropriate limit
            xl = rtflsp;
            flo = f;
        }
        else {
            xh = rtflsp;
            fhi = f;
        }
        dx = xh-xl;
 
        if ((fabs (del) <= xacc)||(f == 0.)) return rtflsp;
    }   // convergence
 
    //  printf  ("RootFlsPos exceeded maximum iterations");
    return 1E308;
}
 
 
double OT_ROOT_FIND_1D::RootNewton (Function fn, Function df,
                                    const double &x1, const double &x2, const double &xacc, void* pvData) {
    const unsigned JMAX = 20;
    double rtnewt = .5*(x1+x2);         // initial guess
    unsigned j;
 
    for (j = 1; j <= JMAX; j++) {
        double f = fn (rtnewt, pvData), fs = df (rtnewt, pvData),
                                             dx = f/fs;
 
        rtnewt -= dx;
        if ((x1-rtnewt)*(rtnewt-x2) < 0.)
        {
            //      printf ("jumped out of brackets");
            return 1E308;
        }
        if (fabs (dx) < xacc) return rtnewt;
    }       // convergence
 
    //    printf("RootNewton exceeded maximum iterations");
    return 1E308;
}
 
double OT_ROOT_FIND_1D::RootSafe (Function fn, Function df,
                                  const double &x1, const double &x2, const double &xacc, void* pvData) {
    const unsigned JMAX = 100;
 
    double fl = fn (x1, pvData), fh = fn (x2, pvData), xl, xh;
 
    if (fl*fh > SMALL_NUM_2)
    {
        //    printf ("root must be bracketed");
        return 1E308;
    }
    if (fl < 0) {                         // orient the search so that f (xl) < 0
        xl = x1;
        xh = x2;
    }
    else {
        xl = x2;
        xh = x1;
    }
 
    double rtsafe = .5*(x1+x2);         // initial guess
    double dxold = fabs (x2-x1),
                   dx = dxold;
    double f = fn (rtsafe, pvData), fp = df (rtsafe, pvData);
 
    for (unsigned j = 1; j <= JMAX; j++) {
        if ( ( ((rtsafe-xh)*fp-f)*((rtsafe-xl)*fp-f) >= 0 ) ||  // bisect if newton out of brackets
                ( fabs (2*f) > fabs (dxold*fp) ) ) {                 // or not decreasing fast enough
            dxold = dx;
            dx = .5*(xh -xl);
            rtsafe = xl+dx;
            if (xl == rtsafe) return rtsafe;
        } // change in root is negligible
        else {
            dxold = dx;
            dx = f/fp;
            double temp = rtsafe;
            rtsafe -= dx;
            if (temp == rtsafe) return rtsafe;
        }
        if (fabs (dx) < xacc) return rtsafe;  // convergence
        f = fn (rtsafe, pvData);
        fp = df (rtsafe, pvData);
        if (f < 0) xl = rtsafe;
        else       xh = rtsafe;
    }
    //   printf("too many iterations");
    return 1E308;
}