/*
Programmazione II - Esempio di codice
<![CDATA[
Example code: function bracketing + bisection zero search routines
(20000814 prelz@mi.infn.it)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct bracket_s
{
double x1;
double x2;
} bracket;
int get_bracket_expand(double (*f)(double x), bracket *bra);
int get_bracket_partition(double (*f)(double x), bracket in_bra,
unsigned int n_partitions, bracket **out_bra);
double get_zero_with_bisection(double (*f)(double x), bracket bra,
double tolerance);
int
main(int argc, char *argv[])
{
double argument;
double result1, result2;
double zero;
bracket test,*ret_bra;
int result;
int i;
test.x1 = 0;
test.x2 = M_PI/4;
result = get_bracket_expand(cos, &test);
printf ("get_bracket_expand cos with initial interval (0,PI/4) returns %d:\n",
result);
printf ("x1 = %g PI, x2 = %g PI", test.x1/M_PI, test.x2/M_PI);
zero = get_zero_with_bisection(cos, test, 1E-6);
printf (" zero = %.6g PI (+/- 1E-6)\n",zero/M_PI);
test.x1 = 0;
test.x2 = 4*M_PI;
result = get_bracket_partition(cos, test, 20, &ret_bra);
printf ("get_bracket_partition cos with initial interval (0,4*PI) returns %d:\n",
result);
for (i=0;i<result;i++)
{
printf ("Interval # %d) x1 = %g PI, x2 = %g PI", i+1, ret_bra[i].x1/M_PI,
ret_bra[i].x2/M_PI);
zero = get_zero_with_bisection(cos, ret_bra[i], 1E-6);
printf (" zero = %.6g PI (+/- 1E-6)\n",zero/M_PI);
}
exit(0);
}
#define EXPANSION_FACTOR 1.6
#define EXPANSION_MAX_RETRIES 50
int
get_bracket_expand(double (*f)(double x), bracket *bra)
{
int i;
double fx1, fx2;
for (i=0;i<EXPANSION_MAX_RETRIES;i++)
{
fx1 = (*f)(bra->x1);
fx2 = (*f)(bra->x2);
if (fx1*fx2 < 0) return(1);
if (fabs(fx1) < fabs(fx2))
{
/* Try expanding the interval at the x1 side */
bra->x1 = bra->x2 + (bra->x1-bra->x2) * EXPANSION_FACTOR;
}
else
{
/* Try expanding the interval at the x2 side */
bra->x2 = bra->x1 + (bra->x2-bra->x1) * EXPANSION_FACTOR;
}
}
return(0); /* 0 means no intervals were found */
}
int
get_bracket_partition(double (*f)(double x), bracket in_bra,
unsigned int n_partitions, bracket **out_bra)
{
static bracket *results = NULL;
int n_results = 0;
double cur_x1,cur_x2;
double segment_size;
int i;
segment_size = (in_bra.x2 - in_bra.x1) / (double)n_partitions;
cur_x1 = in_bra.x1;
cur_x2 = cur_x1 + segment_size;
for (i=0; i<n_partitions; i++)
{
/* Last segment */
if (i == (n_partitions-1)) cur_x2 = in_bra.x2;
if ((*f)(cur_x1)*(*f)(cur_x2) < 0)
{
results = (bracket *)realloc(results, (n_results+1) * sizeof(bracket));
if (results == NULL)
{
fprintf(stderr,"Out of Memory.\n");
return(-1);
}
results[n_results].x1 = cur_x1;
results[n_results].x2 = cur_x2;
n_results++;
}
cur_x1 += segment_size;
cur_x2 += segment_size;
}
*out_bra = results;
return(n_results);
}
#define MAXIMUM_BISECTIONS 100
double
get_zero_with_bisection(double (*f)(double x), bracket bra,
double tolerance)
{
int i;
double xmiddle;
double fleft, fmiddle, fright;
/* Computing f might be expensive */
fleft = (*f)(bra.x1);
fright = (*f)(bra.x2);
if (fleft*fright >= 0)
{
fprintf(stderr,"Error: input interval is not a valid zero bracket.\n");
return(0);
}
for (i=0; i<MAXIMUM_BISECTIONS; i++)
{
xmiddle = (bra.x1 + bra.x2) / 2;
fmiddle = (*f)(xmiddle);
if (fmiddle == 0) return(xmiddle); /* We cannot expect */
/* to be this lucky too often. */
if (fabs(bra.x2 - bra.x1) < fabs(tolerance)) return(xmiddle);
if (fleft * fmiddle < 0)
{
bra.x2 = xmiddle;
fright = fmiddle;
}
else
{
bra.x1 = xmiddle;
fleft = fmiddle;
}
}
/* We get here only if we received a bad tolerance value */
fprintf(stderr,"Error: Too many bisections.\n");
return(0);
}
]]></plaintext></span>