outil/src/ot_root_find.cpp

//---------------------------------------------------------------------------


#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;
}

Revision:	1019
Committed:	Tue Jun 4 21:16:50 2019 UTC (6 years ago) by francois
File size:	8119 byte(s)
Log Message:	restructuration de magic outil est sorti de lib pour pouvoir etre utiliser en dehors de lib template est merge avec outil poly_occ et un sous projet de magic qui utilise le nouveau outil
#	User	Rev	Content
1	francois	283	//---------------------------------------------------------------------------
2
3
4			#include <stdio.h>
5			#include <math.h>
6
7			#pragma hdrstop
8
9			#include "ot_root_find.h"
10
11			//---------------------------------------------------------------------------
12			#pragma package(smart_init)
13
14			#define SMALL_NUM_1 1E-14
15			#define SMALL_NUM_2 1E-28
16
17			unsigned OT_ROOT_FIND_1D::ZeroBracOut (Function fn, double &x1, double &x2, void* pvData) {
18			const double FACTOR = 1.6;
19			const unsigned NTRY = 50;
20
21			double f1 = fn (x1, pvData), f2 = fn (x2, pvData);
22			unsigned j;
23
24			if (x1 == x2)
25			{
26			// printf ("you have to guess an initial range");
27			return 1e8;
28			}
29			for (j = 1; j <= NTRY; j++) {
30			if (f1*f2 < 0.) return 1; // found a root
31			if (fabs (f1) < fabs (f2)) {
32			x1 += FACTOR*(x1-x2); // expand towards x1
33			f1 = fn (x1, pvData);
34			}
35			else {
36			x2 += FACTOR*(x2-x1); // ... or towards x2
37			f2 = fn (x2, pvData);
38			}
39			}
40			return 1E8;
41			} // no root bracketed
42
43			unsigned OT_ROOT_FIND_1D::ZeroBracOut (Function fn, double xStart, double xEnd, double &x1, double &x2, void* pvData) {
44			const double FACTOR = 1.6;
45			const unsigned NTRY = 50;
46
47			if (x1>x2) {
48			double tmp=x1;
49			x1=x2;
50			x2=tmp;
51			}
52			if (xStart>xEnd) {
53			double tmp=xStart;
54			xEnd=xStart;
55			xEnd=tmp;
56			}
57			double deltaT=xEnd-xStart;
58
59			if (x1 < xStart)
60			x1 = xStart;
61			else if (x1 > xEnd)
62			x1 = xEnd;
63			if (x2 < xStart)
64			x2 = xStart;
65			else if (x2 > xEnd)
66			x2 = xEnd;
67
68			if (x1 == x2)
69			{
70			x1 -= .1*deltaT;
71			x1 = (x1>xStart)?x1:xStart;
72			x2 += .1*deltaT;
73			x2 = (x2<xEnd)?x2:xEnd;
74			}
75
76			double f1 = fn (x1, pvData), f2 = fn (x2, pvData);
77			unsigned j;
78
79			if (x1 == x2)
80			{
81			// printf ("you have to guess an initial range");
82			return 1E8;
83			}
84			for (j = 1; j <= NTRY; j++) {
85			if (f1*f2 < 0.) return 1; // found a root
86			if (fabs (f1) < fabs (f2) && x1 > xStart )
87			{
88			if (x1+FACTOR*(x1-x2) > xStart)
89			{
90			x1 += FACTOR*(x1-x2); // expand towards x1
91			}
92			else
93			{
94			x1 = xStart;
95			}
96			f1 = fn (x1, pvData);
97			}
98			else if (x2 < xEnd ) {
99			if (x2+FACTOR*(x2-x1) < xEnd)
100			{
101			x2 += FACTOR*(x2-x1); // ... or towards x2
102			}
103			else
104			{
105			x2 = xEnd;
106			}
107			f2 = fn (x2, pvData);
108			}
109			else
110			break;
111			}
112			return 0;
113			} // no root bracketed
114
115			unsigned OT_ROOT_FIND_1D::ZeroBracIn (Function fn, const double &x1, const double &x2,
116			unsigned n, double xb1[], double xb2[], unsigned nb, void* pvData) {
117			unsigned i, nbb = nb;
118			double x = x1, dx = (x2-x1)/n, // determine appropriate spacing
119			fp = fn (x, pvData);
120
121			nb = 0;
122			for (i = 1; i <= n; i++) { // loop over all intervals
123			double fc;
124
125			x += dx;
126			fc = fn (x, pvData);
127			if (fc*fp < 0.) { // if a sign change occurs, record
128			xb1[nb] = x-dx; // values for the bounds.
129			xb2[nb] = x;
130			nb++;
131			}
132			fp = fc;
133			if (nbb == nb) return nb;
134			}
135			return nb;
136			}
137
138
139			double OT_ROOT_FIND_1D::RootBisect (Function fn, const double &x1, const double &x2, const double &xacc, void* pvData)
140			{
141			const unsigned JMAX = 50;
142			double fmid = fn (x2, pvData), f = fn (x1, pvData), rtbis, dx;
143			unsigned j;
144
145			if (fmid*f >= 0)
146			{
147			// printf ("root must be bracketed for bisection");
148			return 1E308;
149			}
150
151			if (f < 0.) { // orient the search so that f > 0
152			rtbis = x1; // lies at x+dx.
153			dx = x2-x1;
154			}
155			else {
156			rtbis = x2;
157			dx = x1-x2;
158			}
159
160			for (j = 1; j <= JMAX; j++) { // bisection loop
161			double xmid;
162
163			dx *= .5; // here we halve the interval
164			xmid = rtbis+dx; // and update boundaries and
165			fmid = fn (xmid, pvData); // function value
166			if (fmid <= 0.) rtbis = xmid;
167			if ((fabs (dx) <= xacc)\|\|(fmid == 0.)) return rtbis;
168			}
169			// printf ("too many bisections");
170			return 1E308;
171			} // this should happen rarely (2^-50 ~ 10^-15)
172
173
174			double OT_ROOT_FIND_1D::RootFlsPos (Function fn, const double &x1, const double &x2,
175			const double &xacc, void* pvData) {
176			const unsigned MAXIT = 40;
177			double flo = fn (x1, pvData), fhi = fn (x2, pvData), xl, xh, dx, dummy;
178			unsigned j;
179
180			if (flo*fhi >= 0)
181			{
182			// printf ("root must be bracketed for false position");
183			return 1E308;
184			}
185
186			if (flo < 0.) { // let xl correspond to the low side
187			xl = x1;
188			xh = x2;
189			}
190			else {
191			{
192			dummy=flo;
193			flo=fhi;
194			fhi=dummy;
195			}
196			xl = x2;
197			xh = x1;
198			}
199			dx = xh-xl;
200
201			for (j = 1; j <= MAXIT; j++) { // false position loop
202			double rtflsp = xl+dx*flo/(flo-fhi),// increment w/ respect to latest
203			f = fn (rtflsp, pvData), // value
204			del = (f < 0.)? xl-rtflsp: xh-rtflsp;
205
206			if (f < 0.) { // replace appropriate limit
207			xl = rtflsp;
208			flo = f;
209			}
210			else {
211			xh = rtflsp;
212			fhi = f;
213			}
214			dx = xh-xl;
215
216			if ((fabs (del) <= xacc)\|\|(f == 0.)) return rtflsp;
217			} // convergence
218
219			// printf ("RootFlsPos exceeded maximum iterations");
220			return 1E308;
221			}
222
223
224			double OT_ROOT_FIND_1D::RootNewton (Function fn, Function df,
225			const double &x1, const double &x2, const double &xacc, void* pvData) {
226			const unsigned JMAX = 20;
227			double rtnewt = .5*(x1+x2); // initial guess
228			unsigned j;
229
230			for (j = 1; j <= JMAX; j++) {
231			double f = fn (rtnewt, pvData), fs = df (rtnewt, pvData),
232			dx = f/fs;
233
234			rtnewt -= dx;
235			if ((x1-rtnewt)*(rtnewt-x2) < 0.)
236			{
237			// printf ("jumped out of brackets");
238			return 1E308;
239			}
240			if (fabs (dx) < xacc) return rtnewt;
241			} // convergence
242
243			// printf("RootNewton exceeded maximum iterations");
244			return 1E308;
245			}
246
247			double OT_ROOT_FIND_1D::RootSafe (Function fn, Function df,
248			const double &x1, const double &x2, const double &xacc, void* pvData) {
249			const unsigned JMAX = 100;
250
251			double fl = fn (x1, pvData), fh = fn (x2, pvData), xl, xh;
252
253			if (fl*fh > SMALL_NUM_2)
254			{
255			// printf ("root must be bracketed");
256			return 1E308;
257			}
258			if (fl < 0) { // orient the search so that f (xl) < 0
259			xl = x1;
260			xh = x2;
261			}
262			else {
263			xl = x2;
264			xh = x1;
265			}
266
267			double rtsafe = .5*(x1+x2); // initial guess
268			double dxold = fabs (x2-x1),
269			dx = dxold;
270			double f = fn (rtsafe, pvData), fp = df (rtsafe, pvData);
271
272			for (unsigned j = 1; j <= JMAX; j++) {
273			if ( ( ((rtsafe-xh)fp-f)((rtsafe-xl)*fp-f) >= 0 ) \|\| // bisect if newton out of brackets
274			( fabs (2f) > fabs (dxoldfp) ) ) { // or not decreasing fast enough
275			dxold = dx;
276			dx = .5*(xh -xl);
277			rtsafe = xl+dx;
278			if (xl == rtsafe) return rtsafe;
279			} // change in root is negligible
280			else {
281			dxold = dx;
282			dx = f/fp;
283			double temp = rtsafe;
284			rtsafe -= dx;
285			if (temp == rtsafe) return rtsafe;
286			}
287			if (fabs (dx) < xacc) return rtsafe; // convergence
288			f = fn (rtsafe, pvData);
289			fp = df (rtsafe, pvData);
290			if (f < 0) xl = rtsafe;
291			else xh = rtsafe;
292			}
293			// printf("too many iterations");
294			return 1E308;
295			}
296