Subversion Repositories programming

Rev

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

Rev 144 Rev 145
Line 50... Line 50...
50 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
 50 
 * Which is equvalent to: I[n](x) = I[n-2](x) - (2 * (n-1)/x) * I[n-1](x)
51 
 * by simple substitution.
 51 
 * by simple substitution.
52 
 *
 52 
 *
53 
 * This is the recurrence relation for PART #1.
 53 
 * This is the recurrence relation for PART #1.
54 
 */
 54 
 */
55 
float Bessel_BottomUp (int n, float x, float base0, float base1)
 55 
float Bessel_BottomUp (int n, float x, float *constants)
56 
{
 56 
{
57 
    int i = 0;
 57 
    int i = 0;
58 
    float  answer;
 58 
    float  answer;
59 
    float *storage = new float[ARRAY_SIZE];
 59 
    float *storage = new float[ARRAY_SIZE];
60 
 
 60 
 
61 
    // Prime the array
 61 
    // Prime the array
62 
    storage[0] = base0;
 62 
    storage[0] = constants[IDX0];
63 
    storage[1] = base1;
 63 
    storage[1] = constants[IDX1];
64 
 
 64 
 
65 
    // Calculate each value in turn
 65 
    // Calculate each value in turn
66 
    for (i=2; i<=20; i++)
 66 
    for (i=2; i<=20; i++)
67 
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];
 67 
        storage[i] = storage[i-2] - (2.0 * (n-1) / x) * storage[i-1];
68 
 
 68 
 
Line 79... Line 79...
79 
 * Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)
 79 
 * Which is equivalent to: I[n](x) = (2 * (n+1) / x) * I[n+1](x) + I[n+2](x)
80 
 * by simple substitution.
 80 
 * by simple substitution.
81 
 *
 81 
 *
82 
 * This is the recurrence relation for PART #2.
 82 
 * This is the recurrence relation for PART #2.
83 
 */
 83 
 */
84 
float Bessel_TopDown (int n, float x, float base19, float base20)
 84 
float Bessel_TopDown (int n, float x, float *constants)
85 
{
 85 
{
86 
    int i = 18;
 86 
    int i = 18;
87 
    float answer;
 87 
    float answer;
88 
    float *storage = new float[ARRAY_SIZE];
 88 
    float *storage = new float[ARRAY_SIZE];
89 
 
 89 
 
90 
    // Prime the array
 90 
    // Prime the array
91 
    storage[20] = base20;
 91 
    storage[20] = constants[IDX20];
92 
    storage[19] = base19;
 92 
    storage[19] = constants[IDX19];
93 
 
 93 
 
94 
    // Calculate each value in turn
 94 
    // Calculate each value in turn
95 
    for (i=18; i>=0; i--)
 95 
    for (i=18; i>=0; i--)
96 
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];
 96 
        storage[i] = (2.0 * (n+1) / x) * storage[i+1] + storage[i+2];
97 
 
 97 
 
Line 106... Line 106...
106 
/* Print out a table for a given x value, and function.
 106 
/* Print out a table for a given x value, and function.
107 
 * Call like this: print_table (1.0, &BesselFunc);
 107 
 * Call like this: print_table (1.0, &BesselFunc);
108 
 *
 108 
 *
109 
 * It chooses the correct constants appropriately.
 109 
 * It chooses the correct constants appropriately.
110 
 */
 110 
 */
111 
void print_table (float x, float (*bessel_func)(int, float, float, float))
 111 
void print_table (float x, float (*bessel_func)(int, float, float*))
112 
{
 112 
{
113 
    int n;
 113 
    int n;
114 
    float base0, base1, base19, base20;
 - 
 
115 
    float true_val, calc_val, abs_err, rel_err;
 114 
    float true_val, calc_val, abs_err, rel_err;
116 
    float *constants;
 115 
    float *constants;
117 
 
 116 
 
118 
    if (x == 1.0)
 - 
 
119 
        constants = x1_const;
 117 
    if (x == 1.0) constants = x1_const;
120 
 
 - 
 
121 
    if (x == 2.0)
 - 
 
122 
        constants = x2_const;
 118 
    if (x == 2.0) constants = x2_const;
123 
 
 - 
 
124 
    if (x == 5.0)
 - 
 
125 
        constants = x5_const;
 119 
    if (x == 5.0) constants = x5_const;
126 
 
 120 
 
127 
    // Set the output flags
 121 
    // Set the output flags
128 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
 122 
    cout << setiosflags (ios_base::scientific | ios_base::showpoint);
129 
 
 123 
 
- 
 
 124 
    // A line of equal signs
- 
 
 125 
    for (n=0; n<(20+20+20+20+10); n++)
- 
 
 126 
        cout << "=";
- 
 
 127 
 
- 
 
 128 
    cout << endl;
- 
 
 129 
 
- 
 
 130 
    // Print the x value we're using
- 
 
 131 
    cout << "Using x = " << x << endl;
- 
 
 132 
 
130 
    // Print the table header
 133 
    // Print the table header
131 
    cout << setw(10) << "n Val."
 134 
    cout << setw(10) << "n Val."
132 
         << setw(20) << "TRUE"
 135 
         << setw(20) << "TRUE"
133 
         << setw(20) << "COMPUTED"
 136 
         << setw(20) << "COMPUTED"
134 
         << setw(20) << "ABS ERROR"
 137 
         << setw(20) << "ABS ERROR"
Line 139... Line 142...
139 
    for (n=0; n<(20+20+20+20+10); n++)
 142 
    for (n=0; n<(20+20+20+20+10); n++)
140 
        cout << "=";
 143 
        cout << "=";
141 
 
 144 
 
142 
    cout << endl;
 145 
    cout << endl;
143 
 
 146 
 
144 
    // Print out the line for n = 5
 147 
    // Run once for n = 5,10,15,20
145 
    n = 5;
 - 
 
146 
    true_val = constants[IDX5];
 - 
 
147 
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
 - 
 
148 
    abs_err = abs (true_val - calc_val);
 - 
 
149 
    rel_err = abs_err / abs (true_val);
 - 
 
150 
 
 - 
 
151 
    cout << setw(10) << n
 148 
    for (n=5; n<=20; n+=5)
152 
         << setw(20) << true_val
 - 
 
153 
         << setw(20) << calc_val
 - 
 
154 
         << setw(20) << abs_err
 - 
 
155 
         << setw(20) << rel_err
 - 
 
156 
         << endl;
 - 
 
157 
 
 - 
 
158 
    // Print out the line for n = 10
 - 
 
159 
    n = 10;
 149 
    {
160 
    true_val = constants[IDX10];
 150 
        if (n==5 ) true_val = constants[IDX5];
161 
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
 151 
        if (n==10) true_val = constants[IDX10];
162 
    abs_err = abs (true_val - calc_val);
 152 
        if (n==15) true_val = constants[IDX15];
163 
    rel_err = abs_err / abs (true_val);
 - 
 
164 
 
 - 
 
165 
    cout << setw(10) << n
 - 
 
166 
         << setw(20) << true_val
 153 
        if (n==20) true_val = constants[IDX20];
167 
         << setw(20) << calc_val
 - 
 
168 
         << setw(20) << abs_err
 - 
 
169 
         << setw(20) << rel_err
 - 
 
170 
         << endl;
 - 
 
171 
 
 154 
 
172 
    // Print out the line for n = 15
 - 
 
173 
    n = 15;
 - 
 
174 
    true_val = constants[IDX15];
 - 
 
175 
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
 155 
        calc_val = bessel_func(n, x, constants);
176 
    abs_err = abs (true_val - calc_val);
 156 
        abs_err  = abs (true_val - calc_val);
177 
    rel_err = abs_err / abs (true_val);
 157 
        rel_err  = abs_err / abs (true_val);
178 
 
 158 
 
179 
    cout << setw(10) << n
 159 
        cout << setw(10) << n
180 
         << setw(20) << true_val
 160 
             << setw(20) << true_val
181 
         << setw(20) << calc_val
 161 
             << setw(20) << calc_val
182 
         << setw(20) << abs_err
 162 
             << setw(20) << abs_err
183 
         << setw(20) << rel_err
 163 
             << setw(20) << rel_err
184 
         << endl;
 164 
             << endl;
185 
 
 - 
 
186 
    // Print out the line for n = 20
 - 
 
187 
    n = 20;
 165 
    }
188 
    true_val = constants[IDX20];
 - 
 
189 
    calc_val = bessel_func(n, x, constants[IDX0], constants[IDX1]);
 - 
 
190 
    abs_err = abs (true_val - calc_val);
 - 
 
191 
    rel_err = abs_err / abs (true_val);
 - 
 
192 
 
 - 
 
193 
    cout << setw(10) << n
 - 
 
194 
         << setw(20) << true_val
 - 
 
195 
         << setw(20) << calc_val
 - 
 
196 
         << setw(20) << abs_err
 - 
 
197 
         << setw(20) << rel_err
 - 
 
198 
         << endl;
 - 
 
199 
}
 166 
}
200 
 
 167 
 
201 
int main (void)
 168 
int main (void)
202 
{
 169 
{
203 
    print_table (1.0, &Bessel_BottomUp);
 170 
    print_table (1.0, &Bessel_BottomUp);