Subversion Repositories programming

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

#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <iostream>
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.
 *
 * @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 answer = 0.0, xi_last = 0.0, xi_now = 0.0;
    int i = 0;

    xi_last = a + (i * h); // i = 0 here

    for (i=1; i<=n; i++)
    {
        xi_now = a + (i * h);
        answer += (f (xi_last) + f (xi_now) * (xi_now - xi_last));
        xi_last = xi_now; // shift xi's
    }

    return answer;
}

/**
 * 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;

    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;
    int n = 0;

    return 0;
}