Subversion Repositories programming

Rev

Rev 147 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 147 Rev 148
Line -... Line 1...
- 
 
 1 
/* Ira Snyder
- 
 
 2 
 * Project #2
- 
 
 3 
 * Due 2005-10-28
- 
 
 4 
 */
- 
 
 5 
 
- 
 
 6 
/* Copyright: Ira W. Snyder
- 
 
 7 
 * Start Date: 2005-10-23
- 
 
 8 
 * End Date: 2005-10-27
- 
 
 9 
 * License: Public Domain
- 
 
 10 
 */
- 
 
 11 
 
1 
#include <iostream>
 12 
#include <iostream>
2 
#include <iomanip>
 13 
#include <iomanip>
3 
#include <cmath>
 14 
#include <cmath>
- 
 
 15 
#include <cstdio>
4 
using namespace std;
 16 
using namespace std;
5 
 
 17 
 
6 
// Global Variables
 18 
// Global Variables
7 
enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20}; // Constants' Indices
 19 
enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20}; // Constants' Indices
8 
 
 20 
 
Line 20... Line 32...
20 
 
 32 
 
21 
// Global Defines
 33 
// Global Defines
22 
#define ARRAY_SIZE 21
 34 
#define ARRAY_SIZE 21
23 
 
 35 
 
24 
// Function Prototypes
 36 
// Function Prototypes
25 
float Bessel_BottomUp (int n, float x, float *constants);
 37 
float *Bessel_BottomUp (int n, float x, float *constants);
26 
float Bessel_TopDown (int n, float x, float *constants);
 38 
float *Bessel_TopDown (int n, float x, float *constants);
27 
void print_table_bottomup (float x);
 39 
void print_table_bottomup (float x);
28 
void print_table_topdown (float x);
 40 
void print_table_topdown (float x);
29 
void print_equals (int num);
 41 
void print_equals (int num);
30 
 
 42 
 
31 
/* Find the Bessel Function.
 43 
/* Find the Bessel Function.
32 
 * Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)
 44 
 * Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)
33 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
 45 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
34 
 * by simple substitution.
 46 
 * by simple substitution.
35 
 *
 47 
 *
36 
 * This is the recurrence relation for PART #1.
 48 
 * This is the recurrence relation for PART #1.
- 
 
 49 
 *
- 
 
 50 
 * Return values:
- 
 
 51 
 * this returns an array of values, in the form:
- 
 
 52 
 * array[0] = calculated value
- 
 
 53 
 * array[1] = part before the subtraction: I[n-2](x)
- 
 
 54 
 * array[2] = part after  the subtraction: (2(n-1)/x) * I[n-1](x)
37 
 */
 55 
 */
38 
float Bessel_BottomUp (int n, float x, float *constants)
 56 
float *Bessel_BottomUp (int n, float x, float *constants)
39 
{
 57 
{
40 
    int i = 0;
 58 
    int i = 0;
41 
    float  answer;
 59 
    float *answer = new float[3];
42 
    float *storage = new float[ARRAY_SIZE];
 60 
    float *storage = new float[ARRAY_SIZE];
43 
 
 61 
 
44 
    // Prime the array
 62 
    // Prime the array
45 
    storage[0] = constants[IDX0];
 63 
    storage[0] = constants[IDX0];
46 
    storage[1] = constants[IDX1];
 64 
    storage[1] = constants[IDX1];
47 
 
 65 
 
48 
    // Calculate each value in turn
 66 
    // Calculate each value in turn
49 
    for (i=2; i<=n; i++)
 67 
    for (i=2; i<=n; i++)
- 
 
 68 
    {
50 
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];
 69 
        storage[i] = storage[i-2] - (2.0 * (i-1) / x) * storage[i-1];
- 
 
 70 
#ifdef DEBUG
- 
 
 71 
        // This will print out some useful information to find out why this
- 
 
 72 
        // is going right or wrong in it's calculation.
- 
 
 73 
        printf ("storage[%d] = %e - %e * %e = %e - %e = %e\n", i, storage[i-2],
- 
 
 74 
                (2.0 * (i-1)  / x), storage[i-1],
- 
 
 75 
                storage[i-2], (2.0 * (i-1) / x) * storage[i-1], storage[i]);
- 
 
 76 
#endif
- 
 
 77 
    }
51 
 
 78 
 
52 
    // Store the answer, delete the memory
 79 
    // Store the answer, delete the memory
53 
    answer = storage[n];
 80 
    answer[0] = storage[n];
- 
 
 81 
    answer[1] = storage[n-2];
- 
 
 82 
    answer[2] = (2.0 * (n-1) / x) * storage[n-1];
54 
    delete[] storage;
 83 
    delete[] storage;
55 
 
 84 
 
56 
    // Return the result
 85 
    // Return the result
57 
    return answer;
 86 
    return answer;
58 
}
 87 
}
Line 61... Line 90...
61 
 * Algorithm: I[n-1](x) = (2*n/x)*I[n](x) + I[n+1](x)
 90 
 * Algorithm: I[n-1](x) = (2*n/x)*I[n](x) + I[n+1](x)
62 
 * Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)
 91 
 * Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)
63 
 * by simple substitution.
 92 
 * by simple substitution.
64 
 *
 93 
 *
65 
 * This is the recurrence relation for PART #2.
 94 
 * This is the recurrence relation for PART #2.
- 
 
 95 
 *
- 
 
 96 
 * Return values:
- 
 
 97 
 * this returns an array of values in the form:
- 
 
 98 
 * array[0] = calculated value
- 
 
 99 
 * array[1] = part before the addition: I[n+2](x)
- 
 
 100 
 * array[2] = part after  the addition: (2(n+1)/x) * I[n+1](x)
66 
 */
 101 
 */
67 
float Bessel_TopDown (int n, float x, float *constants)
 102 
float *Bessel_TopDown (int n, float x, float *constants)
68 
{
 103 
{
69 
    int i = 18;
 104 
    int i = 18;
70 
    float answer;
 105 
    float *answer = new float[3];
71 
    float *storage = new float[ARRAY_SIZE];
 106 
    float *storage = new float[ARRAY_SIZE];
72 
 
 107 
 
73 
    // Prime the array
 108 
    // Prime the array
74 
    storage[20] = constants[IDX20];
 109 
    storage[20] = constants[IDX20];
75 
    storage[19] = constants[IDX19];
 110 
    storage[19] = constants[IDX19];
76 
 
 111 
 
77 
    // Calculate each value in turn
 112 
    // Calculate each value in turn
78 
    for (i=18; i>=n; i--)
 113 
    for (i=18; i>=n; i--)
- 
 
 114 
    {
79 
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];
 115 
        storage[i] = (2.0 * (i+1) / x) * storage[i+1] + storage[i+2];
- 
 
 116 
#ifdef DEBUG
- 
 
 117 
        // This will print out some useful information to find out why this
- 
 
 118 
        // is going right or wrong in it's calculation.
- 
 
 119 
        printf ("storage[%d] = %e * %e + %e = %e + %e = %e\n", i,
- 
 
 120 
                (2.0 * (i+1) / x), storage[i+1], storage[i+2],
- 
 
 121 
                (2.0 * (i+1) / x) * storage[i+1], storage[i+2], storage[i]);
- 
 
 122 
#endif
- 
 
 123 
    }
80 
 
 124 
 
81 
    // Store the answer, free the memory
 125 
    // Store the answer, free the memory
82 
    answer = storage[n];
 126 
    answer[0] = storage[n];
- 
 
 127 
    answer[1] = storage[n+2];
- 
 
 128 
    answer[2] = (2.0 * (n+1) / x) * storage[n+1];
83 
    delete[] storage;
 129 
    delete[] storage;
84 
 
 130 
 
85 
    // Return the result
 131 
    // Return the result
86 
    return answer;
 132 
    return answer;
87 
}
 133 
}
88 
 
 134 
 
89 
/* Print out a table for a given x value, and function.
 135 
/* Print a table of the Bessel function (for PART #2)
90 
 * Call like this: print_table (1.0, &BesselFunc);
 - 
 
91 
 *
 - 
 
92 
 * It chooses the correct constants appropriately.
 136 
 * for the given value of x */
93 
 */
 - 
 
94 
void print_table_bottomup (float x)
 137 
void print_table_bottomup (float x)
95 
{
 138 
{
96 
    int n, tblwidth = 6+16+16+16+16+16+20;
 139 
    int n, tblwidth = 6+16+16+16+16+16+20;
97 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
 140 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
98 
    float *constants;
 141 
    float *constants, *answer;
99 
 
 142 
 
- 
 
 143 
    // Find the right constants to use for the x value we're using.
100 
    if (x == 1.0) constants = x1_const;
 144 
    if (x == 1.0) constants = x1_const;
101 
    if (x == 2.0) constants = x2_const;
 145 
    if (x == 2.0) constants = x2_const;
102 
    if (x == 5.0) constants = x5_const;
 146 
    if (x == 5.0) constants = x5_const;
103 
 
 147 
 
104 
    // Reset the output flags
 148 
    // Reset the output flags
Line 107... Line 151...
107 
    // A line of equal signs
 151 
    // A line of equal signs
108 
    print_equals (tblwidth);
 152 
    print_equals (tblwidth);
109 
 
 153 
 
110 
    // Print the x value we're using
 154 
    // Print the x value we're using
111 
    cout << "Using x = " << x
 155 
    cout << "Using x = " << x
- 
 
 156 
         << " and working Bottom-up (PART #1)" << endl;
112 
         << "          "
 157 
    // A line of equal signs
113 
         << "Working Bottom-up" << endl;
 158 
    print_equals (tblwidth);
114 
 
 159 
 
115 
    // Set the output flags we're using
 160 
    // Set the output flags we're using
116 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
 161 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
117 
 
 162 
 
118 
    // Print the table header
 163 
    // Print the table header
Line 121... Line 166...
121 
         << setw(16) << "COMPUTED"
 166 
         << setw(16) << "COMPUTED"
122 
         << setw(16) << "ABS ERROR"
 167 
         << setw(16) << "ABS ERROR"
123 
         << setw(16) << "REL ERROR"
 168 
         << setw(16) << "REL ERROR"
124 
         << setw(16) << "I[n-2](x)"
 169 
         << setw(16) << "I[n-2](x)"
125 
         << setw(20) << "(2(n-1)/x)I[n-1](x)"
 170 
         << setw(20) << "(2(n-1)/x)I[n-1](x)"
126 
         << setw(16) << "IraVal"
 - 
 
127 
         << endl;
 171 
         << endl;
128 
 
 172 
 
129 
    // A line of equal signs
 173 
    // A line of equal signs
130 
    print_equals (tblwidth);
 174 
    print_equals (tblwidth);
131 
 
 175 
 
132 
    // Run once for n = 5,10,15,20
 176 
    // Run once for n = 5,10,15,20
133 
    for (n=5; n<=20; n+=5)
 177 
    for (n=5; n<=20; n+=5)
134 
    {
 178 
    {
- 
 
 179 
        // Choose the true values
- 
 
 180 
        switch (n)
- 
 
 181 
        {
135 
        if (n==5 ) true_val = constants[IDX5];
 182 
            case 5:  true_val = constants[IDX5];  break;
136 
        if (n==10) true_val = constants[IDX10];
 183 
            case 10: true_val = constants[IDX10]; break;
137 
        if (n==15) true_val = constants[IDX15];
 184 
            case 15: true_val = constants[IDX15]; break;
138 
        if (n==20) true_val = constants[IDX20];
 185 
            case 20: true_val = constants[IDX20]; break;
- 
 
 186 
        }
139 
 
 187 
 
- 
 
 188 
        // Get all of the pieces to print
140 
        calc_val = Bessel_BottomUp(n, x, constants);
 189 
        answer = Bessel_BottomUp(n, x, constants);
- 
 
 190 
        calc_val = answer[0];
141 
        abs_err  = abs (true_val - calc_val);
 191 
        abs_err  = abs (true_val - calc_val);
142 
        rel_err  = abs_err / abs (true_val);
 192 
        rel_err  = abs_err / abs (true_val);
143 
        part1    = Bessel_BottomUp(n-2, x, constants);
 193 
        part1    = answer[1];
144 
        part2    = (2.0 * (n-1) / x) * Bessel_BottomUp(n-1, x, constants);
 194 
        part2    = answer[2];
- 
 
 195 
        delete[] answer;
145 
 
 196 
 
- 
 
 197 
        // Print a table line
146 
        cout << setw(6) << n
 198 
        cout << setw(6) << n
147 
             << setw(16) << true_val
 199 
             << setw(16) << true_val
148 
             << setw(16) << calc_val
 200 
             << setw(16) << calc_val
149 
             << setw(16) << abs_err
 201 
             << setw(16) << abs_err
150 
             << setw(16) << rel_err
 202 
             << setw(16) << rel_err
151 
             << setw(16) << part1
 203 
             << setw(16) << part1
152 
             << setw(20) << part2
 204 
             << setw(20) << part2
153 
             << setw(16) << part1 - part2
 - 
 
154 
             << endl;
 205 
             << endl;
155 
    }
 206 
    }
156 
 
 207 
 
157 
    cout << endl;
 208 
    cout << endl;
158 
}
 209 
}
159 
 
 210 
 
- 
 
 211 
/* Print a table of the Bessel function (for PART #2)
- 
 
 212 
 * for the given value of x */
160 
void print_table_topdown (float x)
 213 
void print_table_topdown (float x)
161 
{
 214 
{
162 
    int n, tblwidth = 6+16+16+16+16+16+20;
 215 
    int n, tblwidth = 6+16+16+16+16+16+20;
163 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
 216 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
164 
    float *constants;
 217 
    float *constants, *answer;
165 
    bool done = false;
 218 
    bool done = false;
166 
 
 219 
 
- 
 
 220 
    // Choose the correct constants for the value of x we're using
167 
    if (x == 1.0) constants = x1_const;
 221 
    if (x == 1.0) constants = x1_const;
168 
    if (x == 2.0) constants = x2_const;
 222 
    if (x == 2.0) constants = x2_const;
169 
    if (x == 5.0) constants = x5_const;
 223 
    if (x == 5.0) constants = x5_const;
170 
 
 224 
 
171 
    // Reset the output flags
 225 
    // Reset the output flags
Line 174... Line 228...
174 
    // A line of equal signs
 228 
    // A line of equal signs
175 
    print_equals (tblwidth);
 229 
    print_equals (tblwidth);
176 
 
 230 
 
177 
    // Print the x value we're using
 231 
    // Print the x value we're using
178 
    cout << "Using x = " << x
 232 
    cout << "Using x = " << x
- 
 
 233 
         << " and working Top-down (PART #2)" << endl;
179 
         << "          "
 234 
    // A line of equal signs
180 
         << "Working Top-down" << endl;
 235 
    print_equals (tblwidth);
181 
 
 236 
 
182 
    // Set the output flags we're using
 237 
    // Set the output flags we're using
183 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
 238 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
184 
 
 239 
 
185 
    // Print the table header
 240 
    // Print the table header
Line 188... Line 243...
188 
         << setw(16) << "COMPUTED"
 243 
         << setw(16) << "COMPUTED"
189 
         << setw(16) << "ABS ERROR"
 244 
         << setw(16) << "ABS ERROR"
190 
         << setw(16) << "REL ERROR"
 245 
         << setw(16) << "REL ERROR"
191 
         << setw(16) << "I[n+2](x)"
 246 
         << setw(16) << "I[n+2](x)"
192 
         << setw(20) << "(2(n+1)/x)I[n+1](x)"
 247 
         << setw(20) << "(2(n+1)/x)I[n+1](x)"
193 
         << setw(16) << "IraVal"
 - 
 
194 
         << endl;
 248 
         << endl;
195 
 
 249 
 
196 
    // A line of equal signs
 250 
    // A line of equal signs
197 
    print_equals (tblwidth);
 251 
    print_equals (tblwidth);
198 
 
 252 
 
199 
    n = 0;
 253 
    n = 0;
200 
    // Run once for n = 0,1,5,10,15
 254 
    // Run once for n = 0,1,5,10,15
201 
    while (!done)
 255 
    while (!done)
202 
    {
 256 
    {
- 
 
 257 
        // Choose the correct true values of the function
203 
        switch (n)
 258 
        switch (n)
204 
        {
 259 
        {
205 
            case 0:  true_val = constants[IDX0];  break;
 260 
            case 0:  true_val = constants[IDX0];  break;
206 
            case 1:  true_val = constants[IDX1];  break;
 261 
            case 1:  true_val = constants[IDX1];  break;
207 
            case 5:  true_val = constants[IDX5];  break;
 262 
            case 5:  true_val = constants[IDX5];  break;
208 
            case 10: true_val = constants[IDX10]; break;
 263 
            case 10: true_val = constants[IDX10]; break;
209 
            case 15: true_val = constants[IDX15]; break;
 264 
            case 15: true_val = constants[IDX15]; break;
210 
        }
 265 
        }
211 
 
 266 
 
- 
 
 267 
        // Get all of the values to be printed
212 
        calc_val = Bessel_TopDown(n, x, constants);
 268 
        answer = Bessel_TopDown(n, x, constants);
- 
 
 269 
        calc_val = answer[0];
213 
        abs_err  = abs (true_val - calc_val);
 270 
        abs_err  = abs (true_val - calc_val);
214 
        rel_err  = abs_err / abs (true_val);
 271 
        rel_err  = abs_err / abs (true_val);
215 
        part1    = Bessel_TopDown(n+2, x, constants);
 272 
        part1    = answer[1];
216 
        part2    = (2.0 * (n+1) / x) * Bessel_TopDown(n+1, x, constants);
 273 
        part2    = answer[2];
- 
 
 274 
        delete[] answer;
217 
 
 275 
 
- 
 
 276 
        // Print a line in the table
218 
        cout << setw(6) << n
 277 
        cout << setw(6) << n
219 
             << setw(16) << true_val
 278 
             << setw(16) << true_val
220 
             << setw(16) << calc_val
 279 
             << setw(16) << calc_val
221 
             << setw(16) << abs_err
 280 
             << setw(16) << abs_err
222 
             << setw(16) << rel_err
 281 
             << setw(16) << rel_err
223 
             << setw(16) << part1
 282 
             << setw(16) << part1
224 
             << setw(20) << part2
 283 
             << setw(20) << part2
225 
             << setw(16) << (part1 + part2)
 - 
 
226 
             << endl;
 284 
             << endl;
227 
 
 285 
 
228 
        switch (n) // Set n for the next loop
 286 
        switch (n) // Set n for the next loop
229 
        {
 287 
        {
230 
            case 0:  n = 1;  break;
 288 
            case 0:  n = 1;  break;
Line 236... Line 294...
236 
    }
 294 
    }
237 
 
 295 
 
238 
    cout << endl;
 296 
    cout << endl;
239 
}
 297 
}
240 
 
 298 
 
- 
 
 299 
/* Print a line of equal signs on the screen */
241 
void print_equals (int num)
 300 
void print_equals (int num)
242 
{
 301 
{
243 
    int i;
 302 
    int i;
244 
 
 303 
 
245 
    for (i=0; i<num; i++)
 304 
    for (i=0; i<num; i++)
Line 248... Line 307...
248 
    cout << endl;
 307 
    cout << endl;
249 
}
 308 
}
250 
 
 309 
 
251 
int main (void)
 310 
int main (void)
252 
{
 311 
{
- 
 
 312 
#ifdef DEBUG
- 
 
 313 
    // Print out lots of good information to find out why these functions are
- 
 
 314 
    // going right or wrong.
- 
 
 315 
    float *answer;
- 
 
 316 
    cout << "Bottom-up" << endl << endl;
- 
 
 317 
    answer = Bessel_BottomUp(20, 1.0, x1_const);
- 
 
 318 
    cout << endl << endl << "Top-down" << endl << endl;
- 
 
 319 
    answer = Bessel_TopDown (0,  1.0, x1_const);
- 
 
 320 
#else
- 
 
 321 
    // Print the first line of output as my name
- 
 
 322 
    cout << "Name: Ira Snyder" << endl << endl;
- 
 
 323 
 
- 
 
 324 
    // Print all of the tables we need
253 
    print_table_bottomup (1.0);
 325 
    print_table_bottomup (1.0);
254 
    print_table_bottomup (2.0);
 326 
    print_table_bottomup (2.0);
255 
    print_table_bottomup (5.0);
 327 
    print_table_bottomup (5.0);
256 
    print_table_topdown  (1.0);
 328 
    print_table_topdown  (1.0);
257 
    print_table_topdown  (2.0);
 329 
    print_table_topdown  (2.0);
258 
    print_table_topdown  (5.0);
 330 
    print_table_topdown  (5.0);
259 
 
 331 
#endif
260 
    return 0;
 332 
    return 0;
261 
}
 333 
}
262 
 
 334