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()
 *
 */

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

static const float a = 0.0;
static const float b = 3.0;

float f (float x)
{
    return 1.0 + pow (x, 28);
}

float eval_euler (const int n)
{
    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;
}

float eval_trapezoid (const int n)
{
    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;
}

float eval_romberg (const int n)
{
    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];
}

float eval_simpson (const int n)
{

}

int main (void)
{
    int i = 5;
    printf ("i: %d -- euler: %e\n", i, eval_euler (i));
    printf ("i: %d -- trap: %e\n", i, eval_trapezoid (i));
    printf ("i: %d -- romberg: %e\n", i, eval_romberg (i));
    return 0;
}