Subversion Repositories programming

/* Copyright: Ira W. Snyder
 * Start Date: 2005-11-28
 * End Date: 2005-11-29
 * License: Public Domain
 *
 * Changelog Follows:
 * 2005-11-28
 * - Implement eval_euler() and eval_trapezoid().
 *
 * 2005-11-29
 * - Implement eval_romberg() and eval_simpson().
 * - eval_trapezoid() was giving horrible results, fix it.
 *
 */

#include <cmath>
#include <cstdio>
using namespace std;

/**
 * The function given in the Project #5 Handout.
 * This is 1 + x^28
 *
 * @param x the point at which to evaluate the function
 */
float func (float x)
{
    return 1.0 + pow (x, 28);
}

/**
 * Evaluate the integral of the function given, over the interval given,
 * using the Euler method.
 *
 * @param a the lower bound of the integral
 * @param b the upper bound of the integral
 * @param n the number of intervals to use
 * @param f the function to evaluate
 */
float eval_euler (const float a, const float b, const int n, float(*f)(float))
{
    const float h = (b - a) / n;
    float answer = 0.0, xi = 0.0;
    int i = 0;

    for (i=0; i<=n; i++)
    {
        xi = a + (i * h);
        answer += f(xi);
    }

    return h * answer;
}

/**
 * Evaluate the integral of the function given, over the interval given,
 * using the trapezoid method.
 *
 * This is derived from the book's pseudocode on Pg 207.
 *
 * @param a the lower bound of the integral
 * @param b the upper bound of the integral
 * @param n the number of intervals to use
 * @param f the function to evaluate
 */
float eval_trapezoid (const float a, const float b, const int n, float(*f)(float))
{
    const float h = (b - a) / n;
    float sum = 0.0;
    int i = 0;

    sum = (1.0/2.0) * (f(a) + f(b));

    for (i=1; i<=n-1; i++)
        sum += f(a + (i * h));

    return h * sum;
}

/**
 * Evaluate the integral of the function given, over the interval given,
 * using the Romberg method.
 *
 * This is directly based on the pseudocode on Pg 223-224 of the textbook.
 *
 * @param a the lower bound of the integral
 * @param b the upper bound of the integral
 * @param n the number of intervals to use
 * @param f the function to evaluate
 */
float eval_romberg (const float a, const float b, const int n, float(*f)(float))
{
    float R[n+1][n+1];
    int i, j, k;
    float h, sum;

    /* Check for values that make us run wild, and refuse to run */
    if (n > 25)
    {
        printf ("Cowardly refusing to run... You shouldn't need to use a\n"
                "value of n > 25. Remember, the inner loop runs (2^n) - 1.\n");
        return 0.0;
    }

    h = b - a;
    R[0][0] = (h/2.0) * (f(a) + f(b));

    for (i=1; i<=n; i++)
    {
        h = h/2.0;
        sum = 0.0;

        for (k=1; k<=(pow(2.0, i) - 1); k+=2)
            sum = sum + f(a + k * h);

        R[i][0] = (1.0/2.0) * R[i-1][0] + sum * h;

        for (j=1; j<=i; j++)
            R[i][j] = R[i][j-1] + (R[i][j-1] - R[i-1][j-1]) / (pow(4.0, j) - 1);
    }

    return R[n][n];
}

/**
 * Evaluate the integral of the function given, over the interval given,
 * using the Simpson method.
 *
 * This is derived from the formula on Pg 237-238
 *
 * @param a the lower bound of the integral
 * @param b the upper bound of the integral
 * @param n the number of intervals to use
 * @param f the function to evaluate
 */
float eval_simpson (const float a, const float b, const int n, float(*f)(float))
{
    const float h = (b - a) / n;
    float sum1=0.0, sum2=0.0;
    int i;

    /* Get the first sum in the formula on Pg 238 */
    for (i=1, sum1=0.0; i<=n/2; i++)
        sum1 += f (a + (2 * i - 1) * h);

    sum1 *= 4;

    /* Get the second sum in the formula on Pg 238 */
    for (i=1, sum2=0.0; i<=(n-2)/2; i++)
        sum2 += f (a + 2 * i * h);

    sum2 *= 2;

    /* Add it all together, and return it */
    return (h / 3.0) * ((f(a) + f(b)) + sum1 + sum2);
}

int main (void)
{
    const float a = 0.0;
    const float b = 3.0;
    const float real = 3.0 + (pow(3.0, 29.0) / 29.0);
    float h, comp, abs, rel;
    int i, n;

    /* Euler Method Header */
    printf ("Euler's Method\n");
    printf ("i       h               Computed        Abs Error       Rel Error\n");
    printf ("====================================================================\n");

    /* Euler Method Table */
    for (i=1; i<=25; i++)
    {
        h = pow(2.0, -i);
        n = (int)((b - a) / h);
        comp = eval_euler (a, b, n, &func);
        abs = fabs(real - comp);
        rel = fabs(real - comp) / fabs(real);

        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
    }

    /* Trapezoid Method Header */
    printf ("\n\nTrapezoid Method\n");
    printf ("i       h               Computed        Abs Error       Rel Error\n");
    printf ("====================================================================\n");

    /* Trapezoid Method Table */
    for (i=1; i<=25; i++)
    {
        h = pow(2.0, -i);
        n = (int)((b - a) / h);
        comp = eval_trapezoid (a, b, n, &func);
        abs = fabs(real - comp);
        rel = fabs(real - comp) / fabs(real);

        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
    }

    /* Romberg Method Header */
    printf ("\n\nRomberg's Method\n");
    printf ("i       h               Computed        Abs Error       Rel Error\n");
    printf ("====================================================================\n");

    /* Romberg Method Table
     *
     * NOTE: Don't go higher than 3 for this one, since the Romberg method's inner
     *       loop runs (2^24) - 1 = 16777216 times at i=3 and
     *       (2^48) - 1 = 2.81474 * 10^14 times at i = 4.
     */
    for (i=1; i<=3; i++)
    {
        h = pow(2.0, -i);
        n = (int)((b - a) / h);
        comp = eval_romberg (a, b, n, &func);
        abs = fabs(real - comp);
        rel = fabs(real - comp) / fabs(real);

        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
    }

    /* Simpson Method Header */
    printf ("\n\nSimpson's Method\n");
    printf ("i       h               Computed        Abs Error       Rel Error\n");
    printf ("====================================================================\n");

    /* Simpson Method Table */
    for (i=1; i<=25; i++)
    {
        h = pow(2.0, -i);
        n = (int)((b - a) / h);
        comp = eval_simpson (a, b, n, &func);
        abs = fabs(real - comp);
        rel = fabs(real - comp) / fabs(real);

        printf ("%.2d\t%e\t%e\t%e\t%e\n", i, h, comp, abs, rel);
    }

    return 0;
}