/*12.010 Homework #3 Question 1: Error function calculation. Here we compute the error function and we compare it with the C erf result Compile and link with cc -lm hw3_1_06.c -o hw3_1_06 Accuracy may be passed in runstring i.e., hw3_1_06 1.e-9 will evaluate erf with 9 significant digits. */ /* Prototype definitions */ double factorial( int ) ; double myerf(double, double ); double dmyerf(double ); #define PI 3.1415926535897932 #include #include int main(argc, argv) int argc ; // Number of arguments char *argv[] ; // Pointer to arguments { double eps ; // Accuracy of error function calculation double a ; // argument for erf double cerf, merf, derf ; // C-function ERF, My coding of Erf, derivative of Erf /* See how many arguments were passed */ if ( argc > 1 ) { sscanf(argv[1],"%lf",&eps); } else { eps = 1.e-6; // Use default value } /* Now output headers */ printf("---------------------------------------------------------------\n"); printf("TABLE OF ERROR FUNCTION (ERF) within %8.2e AND FIRST DERIVATIVE\n",eps); printf("---------------------------------------------------------------\n"); printf("| Argument | ERF | d(ERF)/dx | ERFC | Error |\n"); printf("|___________|__________|____________|_____________|___________|\n"); /* Loop over argument range */ for ( a = -3.0; a <= 3.0 ; a += 0.25 ) { cerf = erf(a); merf = myerf(a,eps); derf = dmyerf(a); printf("| %6.3f | %8.5f | %9.5f | %8.5f | %9.2e |\n",a, merf, derf, cerf, merf-cerf); } printf("|___________|__________|____________|_____________|___________|\n"); } double myerf(double x, double eps) { // Routine to compute Erf with accuaracy erf double err =1.0 ; // Error in term double erf_sum = 0.0 ; // Sum of terms in ERF function int max_n = 1 ; // Maximum number of terms in Maclaurin series int n ; // Loop variable /* Compute howmany terms we need */ while ( err > eps && max_n < 100 ) { max_n += 1 ; err = fabs(pow(x,(2*max_n+1))/(factorial(max_n)*(2*max_n+1))); } /* Now compute the function value */ for ( n = 0 ; n <= max_n ; ++n ) { erf_sum = erf_sum + pow(-1.0,n)*pow(x,(2*n+1)) / (factorial(n)*(2*n+1)); } return(2*erf_sum/sqrt(PI)); } double dmyerf(double x) { // Routine to compute first derivative of ERF return(2/sqrt(PI)*exp(-x*x)); } double factorial( int n) { // Routine to compute factorial double fs = 1.0 ; // Factorial value int i ; // Loop counter // Loop doing multiplication for ( i = 2; i <= n ; ++i ) { fs = fs*i; } return(fs); }