/*
Programmazione II - Esempio di codice
<![CDATA[
Example code: "safe" Newton-Raphson method demonstration
(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_partition(double (*f)(double x), bracket in_bra,
unsigned int n_partitions, bracket **out_bra);
double get_zero_with_newtraph(double (*f)(double x),
double (*df)(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 = -M_PI/2;
test.x2 = 4*M_PI;
result = get_bracket_partition(sin, test, 20, &ret_bra);
printf ("get_bracket_partition sin 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_newtraph(sin, cos, ret_bra[i], 1E-6);
printf (" zero = %.6g PI (+/- 1E-6)\n",zero/M_PI);
}
exit(0);
}
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_ITERATIONS 100
double
get_zero_with_newtraph(double (*f)(double x),
double (*df)(double x),
bracket bra, double tolerance)
{
int i;
double xmiddle, newx;
double fleft, fmiddle, fright;
double increment,old_increment;
double dfmiddle;
/* 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);
}
xmiddle = (bra.x1 + bra.x2) / 2;
fmiddle = (*f)(xmiddle);
old_increment = fabs(bra.x1 - bra.x2);
increment = old_increment;
for (i=0; i<MAXIMUM_ITERATIONS; i++)
{
if (fmiddle == 0) return(xmiddle); /* We cannot expect */
/* to be this lucky too often. */
/* Compute first derivative at xmiddle */
dfmiddle = (*df)(xmiddle);
if (dfmiddle != 0) newx = xmiddle - fmiddle/dfmiddle;
else newx = xmiddle;
/* If we are trying to move out of the bracket boundaries */
/* (that is, if the tangent line has the same sign at the boundaries)*/
/* or if the tangent line does not decrease at least by twice */
/* as the function value over the current interval, revert to */
/* normal bisection. */
if ( ((xmiddle - bra.x1)*dfmiddle - fmiddle) *
((xmiddle - bra.x2)*dfmiddle - fmiddle) > 0 ||
(fabs(2*fmiddle) > fabs(old_increment*dfmiddle)) )
{
old_increment = increment;
if (fleft * fmiddle < 0)
{
bra.x2 = xmiddle;
fright = fmiddle;
}
else
{
bra.x1 = xmiddle;
fleft = fmiddle;
}
xmiddle = (bra.x1 + bra.x2) / 2;
increment = fabs(bra.x1 - bra.x2) / 2;
fmiddle = (*f)(xmiddle);
}
else
{
old_increment = increment;
increment = fabs(newx - xmiddle);
xmiddle = newx;
fmiddle = (*f)(xmiddle);
if (fleft * fmiddle < 0)
{
bra.x2 = xmiddle;
fright = fmiddle;
}
else
{
bra.x1 = xmiddle;
fleft = fmiddle;
}
}
if (fabs(bra.x2 - bra.x1) < fabs(tolerance)) return(xmiddle);
}
/* We get here only if we received a bad tolerance value */
fprintf(stderr,"Error: Too many bisections.\n");
return(0);
}
]]></plaintext></span>