* * Programmed by Jack W. Crenshaw, Crenshaw Associates * CIS #72325,1327 * * This file defines certain fundamental math functions. * These functions replace functions max, min, abs, and * sqrt as defined in math.h. * * They are more general, and are defined via templates * so can be used without name changes for any numeric * type. * */ #include #include #ifndef JMATH #define JMATH #define FALSE 0 #define TRUE (!FALSE) // Return the maximum of two numbers template T max(T x, T y) { return (x > y)? x: y; }; // Return the minimum of two numbers template T min(T x, T y) { return (x < y)? x: y; }; // Return the absolute value of a number template T abs(T x) { return (x > (T)0)? x: -x; }; // Return a number with the magnitude of the first // argument, and the sign of the second template T signum(T x, T y) { return (y > (T)0)? abs(x): -abs(x); }; // Return a number modulo the second argument template T mod(T x, T y){ if(y == (T)0) return x; int i = (int)(x/y); if(x*y < (T)0)--i; x = x-((T) i)*y; if(x==y)x -= y; return x; } // return the square root of a number double square_root(double x); long double square_root(long double x); #endif /* jmath.cpp * * Programmed by Jack W. Crenshaw, Crenshaw Associates * CIS #72325,1327 * * This file, with jmath.h, defines certain fundamental * math functions. These functions replace functions max, * fmax, min, fmin, abs, fabs, and sqrt as defined in math.h. * * They are more general, and are defined via templates * so can be used without name changes for any numeric * type. * */ #include // Return the square root of a number float square_root(float x){ if(x < 0.0) return 0.0; else return sqrt(x); } double square_root(double x){ if(x < 0.0) return 0.0; else return sqrt(x); } long double square_root(long double x){ if(x < 0.0L) return 0.0L; else return sqrtl(x); } /* jutils.h * * Copyright (c) 2004 Crenshaw Technologies * * This file contains some simple functions that I find useful * in writing and testing software. */ void pause(void); // use with Borland compiler to hold screen in view void nl(void); // newline (function name follows Small C convention) void error(char *s); // show error and continue void errorStop(char *s); // show error and terminate /* constant.cpp * * Copyright (c) 2004 Crenshaw Technologies * * This file initiates a number of useful constants. It calculates the * values of the constants by use of executable assignments at run time. */ #include using namespace std; // Pi and relatives double Pi = 4*atan(1.0); // 180 degrees double TwoPi = 2*Pi; // 360 degrees double PiOver2 = Pi/2; // 90 degrees double PiOver4 = Pi/4; // 45 degrees double PiOver6 = Pi/6; // 30 degrees // conversion factors double RPD = Pi/180; double DPR = 180/Pi; // Square roots double Root2 = sqrt(2.0); double Root3 = sqrt(3.0); double Root5 = sqrt(5.0); double Golden = (Root5 + 1)/2; // Trig values double Sin30 = 0.5; double Cos30 = Root3/2; double Sin45 = Root2/2; // Synonyms #define Cos45 (Sin45) #define Sin60 (Cos30) #define Cos60 (Sin30) #define pi (Pi) #define HalfPi (PiOver2) #define halfpi (PiOver2) #define Radians_Per_Degree (RPD) #define Degrees_Per_Radian (DPR) /* constant.cpp * * Copyright (c) 2004 Crenshaw Technologies * * This file initiates a number of useful constants. It calculates the * values of the constants by use of executable assignments at run time. */ // Pi and relatives extern double Pi; // 180 degrees extern double TwoPi; // 360 degrees extern double PiOver2; // 90 degrees extern double PiOver4; // 45 degrees extern double PiOver6; // 30 degrees // conversion factors extern double RPD; extern double DPR; // Square roots extern double Root2; extern double Root3; extern double Root5; extern double Golden; // Trig values extern double Sin30; extern double Cos30; extern double Sin45; // Synonyms #define Cos45 (Sin45) #define Sin60 (Cos30) #define Cos60 (Sin30) #define pi (Pi) #define HalfPi (PiOver2) #define halfpi (PiOver2) #define Radians_Per_Degree (RPD) #define Degrees_Per_Radian (DPR) /* jutils.cpp * * Copyright (c) 2004 Crenshaw Technologies * * This file contains some simple functions that I find useful * in writing and testing software. */ #include #include using namespace std; typedef unsigned int uint; // let the customer see the screen (use with Borland compilers only) void pause(void){ cout << "Press any key to continue ..."; while(!kbhit()); } // make a new line (function name follows Small C) void nl(void){ cout << endl; } // show error and continue void error(char *s){ cout << s << endl; pause(); } // show error and terminate void errorStop(char *s){ error(s); exit(1); }//--------------------------------------------------------------------------- #ifndef MullerTestH #define MullerTestH //--------------------------------------------------------------------------- #endif //--------------------------------------------------------------------------- #pragma hdrstop #include #include #include "MullerTest.h" #include "constant.h" #include "jutils.h" #include "jmath.h" using namespace std; // The TWBRF test function double f(double x){ return(4*sin(x)+exp(-x/6)); // return x-1; } // This function evaluates a new point, sets the y range, // and tests for convergence double get_y(double x, double (*f)(double), double eps, double &ymax, double &ymin, bool &converged){ double y = f(x); ymax = max(ymax, y); ymin = min(ymin, y); converged = (abs(y) < eps*(ymax-ymin)); return y; } // The iterator // Based upon the old IBM subroutine rtmi.for // It uses alternating bisection and inverse // parabolic interpolation // restructured by Jack Crenshaw 03/2004 double iterate(double (*f)(double), double x1, double x3, double eps, int imax){ double x2, y2; // the middle point in bisections double xm, ym; // the estimated root in interpolations double xeps = eps*abs(x3-x1); // convergence criterion in x double ymin = 1.0e-300; // set up for criterion in y double ymax = -ymin; bool converged = false; // evaluate function at limits double y1 = get_y(x1, f, eps, ymax, ymin, converged); if(y1 == 0) return x1; // only exact zero accepted here double y3 = get_y(x3, f, eps, ymax, ymin, converged); if(y3 == 0) return x3; // and here // initial points must straddle root if(y1 * y3 > 0) errorStop("Iterate: Initial bounds do not straddle root"); // begin main bisection loop for(int i=0; i0){ swap(x1, x3); swap(y1, y3); } // here's where we try parabolic interpolation. // there are two criteria for accepting the point. // if any one of them fails, just bisect again. // note that y21 and y31 cannot be zero here // but y32 might be. double y21=y2-y1; double y32=y3-y2; if(y32 != 0){ double y31=y3-y1; double b=(x2-x1)/y21; double c=(y21-y32)/(y32*y31); // Test for the RTMI condition if(y3*y31 >= 2.0*y2*y21){ xm=x1-b*y1*(1.0-c*y2); ym = get_y(xm, f, eps, ymax, ymin, converged); // perhaps we nailed the root? if(converged){ return xm; } // Relabel the points 9as needed // to keep root between x1 and x2 if(ym*y1<0.0){ x2=xm; y2=ym; } else{ x1=xm; y1=ym; } } } // we didn't do parabolic interpolation. // just relabel points and continue. x3 = x2; y3 = y2; } // end of bisection loop. If we get here, we did not converge. error("Iterate: Failed to converge in imax tries."); return x2; } //--------------------------------------------------------------------------- #pragma package(smart_init) int main(int argc, char *argv[]){ cout << iterate(f, 0, 6, 1.0e-8, 40) << endl; pause(); return 0; }