Rev 147 | Blame | Compare with Previous | Last modification | View Log | RSS feed
/* Ira Snyder* Project #2* Due 2005-10-28*//* Copyright: Ira W. Snyder* Start Date: 2005-10-23* End Date: 2005-10-27* License: Public Domain*/#include <iostream>#include <iomanip>#include <cmath>#include <cstdio>using namespace std;// Global Variablesenum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20}; // Constants' Indicesfloat x1_const[7] = { 1.266065978, 5.651591040E-01, 2.714631560E-04,2.752948040E-10, 2.370463051E-17, 1.587678369E-23,3.966835986E-25 };float x2_const[7] = { 2.279585302, 1.590636855, 9.825679323E-03,3.016963879E-07, 8.139432531E-13, 8.641603385E-18,4.310560576E-19 };float x5_const[7] = { 2.723987182E+01, 2.433564214E+01, 2.157974547,4.580044419E-03, 1.047977675E-6 , 4.078415017E-10,5.024239358E-11 };// Global Defines#define ARRAY_SIZE 21// Function Prototypesfloat *Bessel_BottomUp (int n, float x, float *constants);float *Bessel_TopDown (int n, float x, float *constants);void print_table_bottomup (float x);void print_table_topdown (float x);void print_equals (int num);/* Find the Bessel Function.* Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)* Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)* by simple substitution.** This is the recurrence relation for PART #1.** Return values:* this returns an array of values, in the form:* array[0] = calculated value* array[1] = part before the subtraction: I[n-2](x)* array[2] = part after the subtraction: (2(n-1)/x) * I[n-1](x)*/float *Bessel_BottomUp (int n, float x, float *constants){int i = 0;float *answer = new float[3];float *storage = new float[ARRAY_SIZE];// Prime the arraystorage[0] = constants[IDX0];storage[1] = constants[IDX1];// Calculate each value in turnfor (i=2; i<=n; i++){storage[i] = storage[i-2] - (2.0 * (i-1) / x) * storage[i-1];#ifdef DEBUG// This will print out some useful information to find out why this// is going right or wrong in it's calculation.printf ("storage[%d] = %e - %e * %e = %e - %e = %e\n", i, storage[i-2],(2.0 * (i-1) / x), storage[i-1],storage[i-2], (2.0 * (i-1) / x) * storage[i-1], storage[i]);#endif}// Store the answer, delete the memoryanswer[0] = storage[n];answer[1] = storage[n-2];answer[2] = (2.0 * (n-1) / x) * storage[n-1];delete[] storage;// Return the resultreturn answer;}/* Find the Bessel Function, from the Top -> Down.* Algorithm: I[n-1](x) = (2*n/x)*I[n](x) + I[n+1](x)* Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)* by simple substitution.** This is the recurrence relation for PART #2.** Return values:* this returns an array of values in the form:* array[0] = calculated value* array[1] = part before the addition: I[n+2](x)* array[2] = part after the addition: (2(n+1)/x) * I[n+1](x)*/float *Bessel_TopDown (int n, float x, float *constants){int i = 18;float *answer = new float[3];float *storage = new float[ARRAY_SIZE];// Prime the arraystorage[20] = constants[IDX20];storage[19] = constants[IDX19];// Calculate each value in turnfor (i=18; i>=n; i--){storage[i] = (2.0 * (i+1) / x) * storage[i+1] + storage[i+2];#ifdef DEBUG// This will print out some useful information to find out why this// is going right or wrong in it's calculation.printf ("storage[%d] = %e * %e + %e = %e + %e = %e\n", i,(2.0 * (i+1) / x), storage[i+1], storage[i+2],(2.0 * (i+1) / x) * storage[i+1], storage[i+2], storage[i]);#endif}// Store the answer, free the memoryanswer[0] = storage[n];answer[1] = storage[n+2];answer[2] = (2.0 * (n+1) / x) * storage[n+1];delete[] storage;// Return the resultreturn answer;}/* Print a table of the Bessel function (for PART #2)* for the given value of x */void print_table_bottomup (float x){int n, tblwidth = 6+16+16+16+16+16+20;float true_val, calc_val, abs_err, rel_err, part1, part2;float *constants, *answer;// Find the right constants to use for the x value we're using.if (x == 1.0) constants = x1_const;if (x == 2.0) constants = x2_const;if (x == 5.0) constants = x5_const;// Reset the output flagscout << resetiosflags (ios_base::scientific | ios_base::showpoint);// A line of equal signsprint_equals (tblwidth);// Print the x value we're usingcout << "Using x = " << x<< " and working Bottom-up (PART #1)" << endl;// A line of equal signsprint_equals (tblwidth);// Set the output flags we're usingcout << setiosflags (ios_base::scientific | ios_base::showpoint);// Print the table headercout << setw(6) << "n Val."<< setw(16) << "TRUE VAL"<< setw(16) << "COMPUTED"<< setw(16) << "ABS ERROR"<< setw(16) << "REL ERROR"<< setw(16) << "I[n-2](x)"<< setw(20) << "(2(n-1)/x)I[n-1](x)"<< endl;// A line of equal signsprint_equals (tblwidth);// Run once for n = 5,10,15,20for (n=5; n<=20; n+=5){// Choose the true valuesswitch (n){case 5: true_val = constants[IDX5]; break;case 10: true_val = constants[IDX10]; break;case 15: true_val = constants[IDX15]; break;case 20: true_val = constants[IDX20]; break;}// Get all of the pieces to printanswer = Bessel_BottomUp(n, x, constants);calc_val = answer[0];abs_err = abs (true_val - calc_val);rel_err = abs_err / abs (true_val);part1 = answer[1];part2 = answer[2];delete[] answer;// Print a table linecout << setw(6) << n<< setw(16) << true_val<< setw(16) << calc_val<< setw(16) << abs_err<< setw(16) << rel_err<< setw(16) << part1<< setw(20) << part2<< endl;}cout << endl;}/* Print a table of the Bessel function (for PART #2)* for the given value of x */void print_table_topdown (float x){int n, tblwidth = 6+16+16+16+16+16+20;float true_val, calc_val, abs_err, rel_err, part1, part2;float *constants, *answer;bool done = false;// Choose the correct constants for the value of x we're usingif (x == 1.0) constants = x1_const;if (x == 2.0) constants = x2_const;if (x == 5.0) constants = x5_const;// Reset the output flagscout << resetiosflags (ios_base::scientific | ios_base::showpoint);// A line of equal signsprint_equals (tblwidth);// Print the x value we're usingcout << "Using x = " << x<< " and working Top-down (PART #2)" << endl;// A line of equal signsprint_equals (tblwidth);// Set the output flags we're usingcout << setiosflags (ios_base::scientific | ios_base::showpoint);// Print the table headercout << setw(6) << "n Val."<< setw(16) << "TRUE VAL"<< setw(16) << "COMPUTED"<< setw(16) << "ABS ERROR"<< setw(16) << "REL ERROR"<< setw(16) << "I[n+2](x)"<< setw(20) << "(2(n+1)/x)I[n+1](x)"<< endl;// A line of equal signsprint_equals (tblwidth);n = 0;// Run once for n = 0,1,5,10,15while (!done){// Choose the correct true values of the functionswitch (n){case 0: true_val = constants[IDX0]; break;case 1: true_val = constants[IDX1]; break;case 5: true_val = constants[IDX5]; break;case 10: true_val = constants[IDX10]; break;case 15: true_val = constants[IDX15]; break;}// Get all of the values to be printedanswer = Bessel_TopDown(n, x, constants);calc_val = answer[0];abs_err = abs (true_val - calc_val);rel_err = abs_err / abs (true_val);part1 = answer[1];part2 = answer[2];delete[] answer;// Print a line in the tablecout << setw(6) << n<< setw(16) << true_val<< setw(16) << calc_val<< setw(16) << abs_err<< setw(16) << rel_err<< setw(16) << part1<< setw(20) << part2<< endl;switch (n) // Set n for the next loop{case 0: n = 1; break;case 1: n = 5; break;case 5: n = 10; break;case 10: n = 15; break;case 15: done = true; break;}}cout << endl;}/* Print a line of equal signs on the screen */void print_equals (int num){int i;for (i=0; i<num; i++)cout << "=";cout << endl;}int main (void){#ifdef DEBUG// Print out lots of good information to find out why these functions are// going right or wrong.float *answer;cout << "Bottom-up" << endl << endl;answer = Bessel_BottomUp(20, 1.0, x1_const);cout << endl << endl << "Top-down" << endl << endl;answer = Bessel_TopDown (0, 1.0, x1_const);#else// Print the first line of output as my namecout << "Name: Ira Snyder" << endl << endl;// Print all of the tables we needprint_table_bottomup (1.0);print_table_bottomup (2.0);print_table_bottomup (5.0);print_table_topdown (1.0);print_table_topdown (2.0);print_table_topdown (5.0);#endifreturn 0;}