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;

}