Subversion Repositories programming

Rev

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

Rev 168 Rev 169
Line 6... Line 6...
6 
 * Changelog Follows:
 6 
 * Changelog Follows:
7 
 * 2005-11-28
 7 
 * 2005-11-28
8 
 * - Implement eval_euler() and eval_trapezoid().
 8 
 * - Implement eval_euler() and eval_trapezoid().
9 
 *
 9 
 *
10 
 * 2005-11-29
 10 
 * 2005-11-29
11 
 * - Implement eval_romberg()
 11 
 * - Implement eval_romberg() and eval_simpson().
12 
 *
 12 
 *
13 
 */
 13 
 */
14 
 
 14 
 
15 
#include <cmath>
 15 
#include <cmath>
16 
#include <cstdlib>
 16 
#include <cstdlib>
17 
#include <cstdio>
 17 
#include <cstdio>
18 
#include <iostream>
 18 
#include <iostream>
19 
using namespace std;
 19 
using namespace std;
20 
 
 20 
 
- 
 
 21 
/**
21 
static const float a = 0.0;
 22 
 * The function given in the Project #5 Handout.
22 
static const float b = 3.0;
 23 
 * This is 1 + x^28
23 
 
 24 
 *
- 
 
 25 
 * @param x the point at which to evaluate the function
- 
 
 26 
 */
24 
float f (float x)
 27 
float func (float x)
25 
{
 28 
{
26 
    return 1.0 + pow (x, 28);
 29 
    return 1.0 + pow (x, 28);
27 
}
 30 
}
28 
 
 31 
 
- 
 
 32 
/**
- 
 
 33 
 * Evaluate the integral of the function given, over the interval given,
- 
 
 34 
 * using the Euler method.
- 
 
 35 
 *
29 
float eval_euler (const int n)
 36 
 * @param a the lower bound of the integral
- 
 
 37 
 * @param b the upper bound of the integral
- 
 
 38 
 * @param n the number of intervals to use
- 
 
 39 
 * @param f the function to evaluate
- 
 
 40 
 */
- 
 
 41 
float eval_euler (const float a, const float b, const int n, float(*f)(float))
30 
{
 42 
{
31 
    const float h = (b - a) / n;
 43 
    const float h = (b - a) / n;
32 
    float answer = 0.0, xi = 0.0;
 44 
    float answer = 0.0, xi = 0.0;
33 
    int i = 0;
 45 
    int i = 0;
34 
 
 46 
 
Line 39... Line 51...
39 
    }
 51 
    }
40 
 
 52 
 
41 
    return h * answer;
 53 
    return h * answer;
42 
}
 54 
}
43 
 
 55 
 
- 
 
 56 
/**
- 
 
 57 
 * Evaluate the integral of the function given, over the interval given,
44 
float eval_trapezoid (const int n)
 58 
 * using the trapezoid method.
- 
 
 59 
 *
- 
 
 60 
 * @param a the lower bound of the integral
- 
 
 61 
 * @param b the upper bound of the integral
- 
 
 62 
 * @param n the number of intervals to use
- 
 
 63 
 * @param f the function to evaluate
- 
 
 64 
 */
- 
 
 65 
float eval_trapezoid (const float a, const float b, const int n, float(*f)(float))
45 
{
 66 
{
46 
    const float h = (b - a) / n;
 67 
    const float h = (b - a) / n;
47 
    float answer = 0.0, xi_last = 0.0, xi_now = 0.0;
 68 
    float answer = 0.0, xi_last = 0.0, xi_now = 0.0;
48 
    int i = 0;
 69 
    int i = 0;
49 
 
 70 
 
Line 57... Line 78...
57 
    }
 78 
    }
58 
 
 79 
 
59 
    return answer;
 80 
    return answer;
60 
}
 81 
}
61 
 
 82 
 
- 
 
 83 
/**
- 
 
 84 
 * Evaluate the integral of the function given, over the interval given,
62 
float eval_romberg (const int n)
 85 
 * using the Romberg method.
- 
 
 86 
 *
- 
 
 87 
 * This is directly based on the pseudocode on Pg 223-224 of the textbook.
- 
 
 88 
 *
- 
 
 89 
 * @param a the lower bound of the integral
- 
 
 90 
 * @param b the upper bound of the integral
- 
 
 91 
 * @param n the number of intervals to use
- 
 
 92 
 * @param f the function to evaluate
- 
 
 93 
 */
- 
 
 94 
float eval_romberg (const float a, const float b, const int n, float(*f)(float))
63 
{
 95 
{
64 
    float R[n+1][n+1];
 96 
    float R[n+1][n+1];
65 
    int i, j, k;
 97 
    int i, j, k;
66 
    float h, sum;
 98 
    float h, sum;
67 
 
 99 
 
Line 83... Line 115...
83 
    }
 115 
    }
84 
 
 116 
 
85 
    return R[n][n];
 117 
    return R[n][n];
86 
}
 118 
}
87 
 
 119 
 
- 
 
 120 
/**
- 
 
 121 
 * Evaluate the integral of the function given, over the interval given,
88 
float eval_simpson (const int n)
 122 
 * using the Simpson method.
- 
 
 123 
 *
- 
 
 124 
 * This is derived from the formula on Pg 237-238
- 
 
 125 
 *
- 
 
 126 
 * @param a the lower bound of the integral
- 
 
 127 
 * @param b the upper bound of the integral
- 
 
 128 
 * @param n the number of intervals to use
- 
 
 129 
 * @param f the function to evaluate
- 
 
 130 
 */
- 
 
 131 
float eval_simpson (const float a, const float b, const int n, float(*f)(float))
89 
{
 132 
{
- 
 
 133 
    const float h = (b - a) / n;
- 
 
 134 
    float sum1=0.0, sum2=0.0;
- 
 
 135 
    int i;
- 
 
 136 
 
- 
 
 137 
    /* Get the first sum in the formula on Pg 238 */
- 
 
 138 
    for (i=1, sum1=0.0; i<=n/2; i++)
- 
 
 139 
        sum1 += f (a + (2 * i - 1) * h);
90 
 
 140 
 
- 
 
 141 
    sum1 *= 4;
- 
 
 142 
 
- 
 
 143 
    /* Get the second sum in the formula on Pg 238 */
- 
 
 144 
    for (i=1, sum2=0.0; i<=(n-2)/2; i++)
- 
 
 145 
        sum2 += f (a + 2 * i * h);
- 
 
 146 
 
- 
 
 147 
    sum2 *= 2;
- 
 
 148 
 
- 
 
 149 
    /* Add it all together, and return it */
- 
 
 150 
    return (h / 3.0) * ((f(a) + f(b)) + sum1 + sum2);
91 
}
 151 
}
92 
 
 152 
 
93 
int main (void)
 153 
int main (void)
94 
{
 154 
{
95 
    int i = 5;
 155 
    const float a = 0.0;
96 
    printf ("i: %d -- euler: %e\n", i, eval_euler (i));
 156 
    const float b = 3.0;
97 
    printf ("i: %d -- trap: %e\n", i, eval_trapezoid (i));
 157 
    int n = 0;
98 
    printf ("i: %d -- romberg: %e\n", i, eval_romberg (i));
 - 
 
- 
 
 158 
 
99 
    return 0;
 159 
    return 0;
100 
}
 160 
}
101 
 
 161