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;

}