Subversion Repositories programming

Rev

Rev 146 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 146 Rev 147
Line 1... Line 1...
1 
#include <iostream>
 1 
#include <iostream>
2 
#include <iomanip>
 2 
#include <iomanip>
3 
#include <cmath>
 3 
#include <cmath>
4 
using namespace std;
 4 
using namespace std;
5 
 
 5 
 
- 
 
 6 
// Global Variables
6 
enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20};
 7 
enum {IDX0, IDX1, IDX5, IDX10, IDX15, IDX19, IDX20}; // Constants' Indices
7 
 
 8 
 
8 
#define ARRAY_SIZE 21
 - 
 
9 
float x1_const[7] = { 1.266065978,     5.651591040E-01, 2.714631560E-04,
 9 
float x1_const[7] = { 1.266065978,     5.651591040E-01, 2.714631560E-04,
10 
                      2.752948040E-10, 2.370463051E-17, 1.587678369E-23,
 10 
                      2.752948040E-10, 2.370463051E-17, 1.587678369E-23,
11 
                      3.966835986E-25 };
 11 
                      3.966835986E-25 };
12 
 
 12 
 
13 
float x2_const[7] = { 2.279585302,     1.590636855,     9.825679323E-03,
 13 
float x2_const[7] = { 2.279585302,     1.590636855,     9.825679323E-03,
Line 16... Line 16...
16 
 
 16 
 
17 
float x5_const[7] = { 2.723987182E+01, 2.433564214E+01, 2.157974547,
 17 
float x5_const[7] = { 2.723987182E+01, 2.433564214E+01, 2.157974547,
18 
                      4.580044419E-03, 1.047977675E-6 , 4.078415017E-10,
 18 
                      4.580044419E-03, 1.047977675E-6 , 4.078415017E-10,
19 
                      5.024239358E-11 };
 19 
                      5.024239358E-11 };
20 
 
 20 
 
- 
 
 21 
// Global Defines
- 
 
 22 
#define ARRAY_SIZE 21
- 
 
 23 
 
- 
 
 24 
// Function Prototypes
- 
 
 25 
float Bessel_BottomUp (int n, float x, float *constants);
- 
 
 26 
float Bessel_TopDown (int n, float x, float *constants);
- 
 
 27 
void print_table_bottomup (float x);
- 
 
 28 
void print_table_topdown (float x);
- 
 
 29 
void print_equals (int num);
- 
 
 30 
 
21 
/* Find the Bessel Function.
 31 
/* Find the Bessel Function.
22 
 * Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)
 32 
 * Algorithm: I[n+1](x) = I[n-1](x) - (2*n/x) * I[n](x)
23 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
 33 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
24 
 * by simple substitution.
 34 
 * by simple substitution.
25 
 *
 35 
 *
Line 34... Line 44...
34 
    // Prime the array
 44 
    // Prime the array
35 
    storage[0] = constants[IDX0];
 45 
    storage[0] = constants[IDX0];
36 
    storage[1] = constants[IDX1];
 46 
    storage[1] = constants[IDX1];
37 
 
 47 
 
38 
    // Calculate each value in turn
 48 
    // Calculate each value in turn
39 
    for (i=2; i<=20; i++)
 49 
    for (i=2; i<=n; i++)
40 
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];
 50 
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];
41 
 
 51 
 
42 
    // Store the answer, delete the memory
 52 
    // Store the answer, delete the memory
43 
    answer = storage[n];
 53 
    answer = storage[n];
44 
    delete[] storage;
 54 
    delete[] storage;
Line 63... Line 73...
63 
    // Prime the array
 73 
    // Prime the array
64 
    storage[20] = constants[IDX20];
 74 
    storage[20] = constants[IDX20];
65 
    storage[19] = constants[IDX19];
 75 
    storage[19] = constants[IDX19];
66 
 
 76 
 
67 
    // Calculate each value in turn
 77 
    // Calculate each value in turn
68 
    for (i=18; i>=0; i--)
 78 
    for (i=18; i>=n; i--)
69 
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];
 79 
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];
70 
 
 80 
 
71 
    // Store the answer, free the memory
 81 
    // Store the answer, free the memory
72 
    answer = storage[n];
 82 
    answer = storage[n];
73 
    delete[] storage;
 83 
    delete[] storage;
Line 79... Line 89...
79 
/* Print out a table for a given x value, and function.
 89 
/* Print out a table for a given x value, and function.
80 
 * Call like this: print_table (1.0, &BesselFunc);
 90 
 * Call like this: print_table (1.0, &BesselFunc);
81 
 *
 91 
 *
82 
 * It chooses the correct constants appropriately.
 92 
 * It chooses the correct constants appropriately.
83 
 */
 93 
 */
84 
void print_table (float x, float (*bessel_func)(int, float, float*))
 94 
void print_table_bottomup (float x)
85 
{
 95 
{
86 
    int n;
 96 
    int n, tblwidth = 6+16+16+16+16+16+20;
87 
    float true_val, calc_val, abs_err, rel_err;
 97 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
88 
    float *constants;
 98 
    float *constants;
89 
 
 99 
 
90 
    if (x == 1.0) constants = x1_const;
 100 
    if (x == 1.0) constants = x1_const;
91 
    if (x == 2.0) constants = x2_const;
 101 
    if (x == 2.0) constants = x2_const;
92 
    if (x == 5.0) constants = x5_const;
 102 
    if (x == 5.0) constants = x5_const;
93 
 
 103 
 
94 
    // Set the output flags
 104 
    // Reset the output flags
95 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
 105 
    cout << resetiosflags (ios_base::scientific | ios_base::showpoint);
96 
 
 106 
 
97 
    // A line of equal signs
 107 
    // A line of equal signs
98 
    for (n=0; n<(20+20+20+20+10); n++)
 - 
 
99 
        cout << "=";
 108 
    print_equals (tblwidth);
100 
 
 - 
 
101 
    cout << endl;
 - 
 
102 
 
 109 
 
103 
    // Print the x value we're using
 110 
    // Print the x value we're using
104 
    cout << "Using x = " << x << endl;
 111 
    cout << "Using x = " << x
- 
 
 112 
         << "          "
- 
 
 113 
         << "Working Bottom-up" << endl;
- 
 
 114 
 
- 
 
 115 
    // Set the output flags we're using
- 
 
 116 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
105 
 
 117 
 
106 
    // Print the table header
 118 
    // Print the table header
107 
    cout << setw(10) << "n Val."
 119 
    cout << setw(6) << "n Val."
108 
         << setw(20) << "TRUE"
 120 
         << setw(16) << "TRUE VAL"
109 
         << setw(20) << "COMPUTED"
 121 
         << setw(16) << "COMPUTED"
110 
         << setw(20) << "ABS ERROR"
 122 
         << setw(16) << "ABS ERROR"
111 
         << setw(20) << "REL ERROR"
 123 
         << setw(16) << "REL ERROR"
- 
 
 124 
         << setw(16) << "I[n-2](x)"
- 
 
 125 
         << setw(20) << "(2(n-1)/x)I[n-1](x)"
- 
 
 126 
         << setw(16) << "IraVal"
112 
         << endl;
 127 
         << endl;
113 
 
 128 
 
114 
    // A line of equal signs
 129 
    // A line of equal signs
115 
    for (n=0; n<(20+20+20+20+10); n++)
 - 
 
116 
        cout << "=";
 130 
    print_equals (tblwidth);
117 
 
 - 
 
118 
    cout << endl;
 - 
 
119 
 
 131 
 
120 
    // Run once for n = 5,10,15,20
 132 
    // Run once for n = 5,10,15,20
121 
    for (n=5; n<=20; n+=5)
 133 
    for (n=5; n<=20; n+=5)
122 
    {
 134 
    {
123 
        if (n==5 ) true_val = constants[IDX5];
 135 
        if (n==5 ) true_val = constants[IDX5];
124 
        if (n==10) true_val = constants[IDX10];
 136 
        if (n==10) true_val = constants[IDX10];
125 
        if (n==15) true_val = constants[IDX15];
 137 
        if (n==15) true_val = constants[IDX15];
126 
        if (n==20) true_val = constants[IDX20];
 138 
        if (n==20) true_val = constants[IDX20];
127 
 
 139 
 
128 
        calc_val = bessel_func(n, x, constants);
 140 
        calc_val = Bessel_BottomUp(n, x, constants);
129 
        abs_err  = abs (true_val - calc_val);
 141 
        abs_err  = abs (true_val - calc_val);
130 
        rel_err  = abs_err / abs (true_val);
 142 
        rel_err  = abs_err / abs (true_val);
- 
 
 143 
        part1    = Bessel_BottomUp(n-2, x, constants);
- 
 
 144 
        part2    = (2.0 * (n-1) / x) * Bessel_BottomUp(n-1, x, constants);
131 
 
 145 
 
132 
        cout << setw(10) << n
 146 
        cout << setw(6) << n
133 
             << setw(20) << true_val
 147 
             << setw(16) << true_val
134 
             << setw(20) << calc_val
 148 
             << setw(16) << calc_val
135 
             << setw(20) << abs_err
 149 
             << setw(16) << abs_err
136 
             << setw(20) << rel_err
 150 
             << setw(16) << rel_err
- 
 
 151 
             << setw(16) << part1
- 
 
 152 
             << setw(20) << part2
- 
 
 153 
             << setw(16) << part1 - part2
137 
             << endl;
 154 
             << endl;
138 
    }
 155 
    }
139 
 
 156 
 
140 
    cout << endl;
 157 
    cout << endl;
141 
}
 158 
}
142 
 
 159 
 
- 
 
 160 
void print_table_topdown (float x)
- 
 
 161 
{
- 
 
 162 
    int n, tblwidth = 6+16+16+16+16+16+20;
- 
 
 163 
    float true_val, calc_val, abs_err, rel_err, part1, part2;
- 
 
 164 
    float *constants;
- 
 
 165 
    bool done = false;
- 
 
 166 
 
- 
 
 167 
    if (x == 1.0) constants = x1_const;
- 
 
 168 
    if (x == 2.0) constants = x2_const;
- 
 
 169 
    if (x == 5.0) constants = x5_const;
- 
 
 170 
 
- 
 
 171 
    // Reset the output flags
- 
 
 172 
    cout << resetiosflags (ios_base::scientific | ios_base::showpoint);
- 
 
 173 
 
- 
 
 174 
    // A line of equal signs
- 
 
 175 
    print_equals (tblwidth);
- 
 
 176 
 
- 
 
 177 
    // Print the x value we're using
- 
 
 178 
    cout << "Using x = " << x
- 
 
 179 
         << "          "
- 
 
 180 
         << "Working Top-down" << endl;
- 
 
 181 
 
- 
 
 182 
    // Set the output flags we're using
- 
 
 183 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
- 
 
 184 
 
- 
 
 185 
    // Print the table header
- 
 
 186 
    cout << setw(6) << "n Val."
- 
 
 187 
         << setw(16) << "TRUE VAL"
- 
 
 188 
         << setw(16) << "COMPUTED"
- 
 
 189 
         << setw(16) << "ABS ERROR"
- 
 
 190 
         << setw(16) << "REL ERROR"
- 
 
 191 
         << setw(16) << "I[n+2](x)"
- 
 
 192 
         << setw(20) << "(2(n+1)/x)I[n+1](x)"
- 
 
 193 
         << setw(16) << "IraVal"
- 
 
 194 
         << endl;
- 
 
 195 
 
- 
 
 196 
    // A line of equal signs
- 
 
 197 
    print_equals (tblwidth);
- 
 
 198 
 
- 
 
 199 
    n = 0;
- 
 
 200 
    // Run once for n = 0,1,5,10,15
- 
 
 201 
    while (!done)
- 
 
 202 
    {
- 
 
 203 
        switch (n)
- 
 
 204 
        {
- 
 
 205 
            case 0:  true_val = constants[IDX0];  break;
- 
 
 206 
            case 1:  true_val = constants[IDX1];  break;
- 
 
 207 
            case 5:  true_val = constants[IDX5];  break;
- 
 
 208 
            case 10: true_val = constants[IDX10]; break;
- 
 
 209 
            case 15: true_val = constants[IDX15]; break;
- 
 
 210 
        }
- 
 
 211 
 
- 
 
 212 
        calc_val = Bessel_TopDown(n, x, constants);
- 
 
 213 
        abs_err  = abs (true_val - calc_val);
- 
 
 214 
        rel_err  = abs_err / abs (true_val);
- 
 
 215 
        part1    = Bessel_TopDown(n+2, x, constants);
- 
 
 216 
        part2    = (2.0 * (n+1) / x) * Bessel_TopDown(n+1, x, constants);
- 
 
 217 
 
- 
 
 218 
        cout << setw(6) << n
- 
 
 219 
             << setw(16) << true_val
- 
 
 220 
             << setw(16) << calc_val
- 
 
 221 
             << setw(16) << abs_err
- 
 
 222 
             << setw(16) << rel_err
- 
 
 223 
             << setw(16) << part1
- 
 
 224 
             << setw(20) << part2
- 
 
 225 
             << setw(16) << (part1 + part2)
- 
 
 226 
             << endl;
- 
 
 227 
 
- 
 
 228 
        switch (n) // Set n for the next loop
- 
 
 229 
        {
- 
 
 230 
            case 0:  n = 1;  break;
- 
 
 231 
            case 1:  n = 5;  break;
- 
 
 232 
            case 5:  n = 10; break;
- 
 
 233 
            case 10: n = 15; break;
- 
 
 234 
            case 15: done = true; break;
- 
 
 235 
        }
- 
 
 236 
    }
- 
 
 237 
 
- 
 
 238 
    cout << endl;
- 
 
 239 
}
- 
 
 240 
 
- 
 
 241 
void print_equals (int num)
- 
 
 242 
{
- 
 
 243 
    int i;
- 
 
 244 
 
- 
 
 245 
    for (i=0; i<num; i++)
- 
 
 246 
        cout << "=";
- 
 
 247 
 
- 
 
 248 
    cout << endl;
- 
 
 249 
}
- 
 
 250 
 
143 
int main (void)
 251 
int main (void)
144 
{
 252 
{
145 
    print_table (1.0, &Bessel_BottomUp);
 253 
    print_table_bottomup (1.0);
146 
    print_table (2.0, &Bessel_BottomUp);
 254 
    print_table_bottomup (2.0);
147 
    print_table (5.0, &Bessel_BottomUp);
 255 
    print_table_bottomup (5.0);
148 
    cout << endl << endl;
 - 
 
149 
    print_table (1.0, &Bessel_TopDown);
 256 
    print_table_topdown  (1.0);
150 
    print_table (2.0, &Bessel_TopDown);
 257 
    print_table_topdown  (2.0);
151 
    print_table (5.0, &Bessel_TopDown);
 258 
    print_table_topdown  (5.0);
152 
 
 259 
 
153 
    return 0;
 260 
    return 0;
154 
}
 261 
}
155 
 
 262