Subversion Repositories programming

/* Copyright: Ira W. Snyder
 * Start Date: 2005-11-13
 * End Date: 2005-11-13
 * License: Public Domain
 *
 * Changelog Follows:
 *
 * 2005-11-13
 * - Implemented Functions 1-3 from the Project #3 Handout.
 * - Implemented Algorithm 7 (LaGrange) from Notes Set #3.
 * - Implemented Algorithm 8-9 (Newton Divided Diff) from Notes Set #3.
 * - Implemented GenEvenPts() which will generate evenly spaced points.
 * - Implemented GenChebychevPts() which will generate Chebychev points.
 *
 */

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

#define X_MIN = -5.0
#define X_MAX = 5.0

/**
 * Function #1 from the Project #3 Handout.
 *
 * @param x the place to calculate the value of this function
 * @return the value of the function at x
 */
double func1 (double x)
{
    return 1.0 / (1.0 + pow(x, 2));
}

/**
 * Function #2 from the Project #3 Handout.
 *
 * @param x the place to calculate the value of this function
 * @return the value of the function at x
 */
double func2 (double x)
{
    // NOTE: M_E is the value of e defined by the math.h header
    return (pow(M_E, x) + pow(M_E, -x)) / 2.0;
}

/**
 * Function #3 from the Project #3 Handout.
 *
 * @param x the place to calculate the value of this function
 * @return the value of the function at x
 */
double func3 (double x)
{
    return pow(x, 14) + pow(x, 10) + pow(x, 4) + x + 1;
}

/**
 * Find the value of the LaGrange Interpolating Polynomial at a point.
 *
 * @param x the x[i] values
 * @param y the y[i] values
 * @param n the number of points
 * @param point the point to evaluate at
 * @return the value of the Interpolating Poly at point
 */
double LaGrangeMethod (double *x, double *y, int n, double point)
{
    double value = 0;
    double term  = 0;
    int i, j;

    for (i=0; i<n; i++)
    {
        term = y[i];

        for (j=0; j<n; j++)
            if (i != j)
                term *= (point - x[j])/(x[i] - x[j]);

        value += term;
    }

    return value;
}

/**
 * Find the value of the Newton Interpolating Polynomial at a point.
 * The a[i]'s are found using a divided difference table.
 *
 * @param x the x[i] values
 * @param y the y[i] values
 * @param n the number of points
 * @param point the point to evaluate at
 * @return the value of the Interpolating Poly at point
 */
double NewtonMethod (double *x, double *y, int n, double point)
{
    int i, j;
    double a[n];

    for (i=0; i<n; i++)
        a[i] = y[i];

    for (i=1; i<n; i++)
        for (j=n-1; j>=i; j--)
            a[j] = (a[j] - a[j-1]) / (x[j] - x[j-i]);

    // At this point, all of the a[i]'s have been calculated,
    // using a Divided Difference Table.

    double xterm = 1.0;
    double value = 0.0;

    for (i=0; i<n; i++)
    {
        value += a[i] * xterm;
        xterm *= (point - x[i]);
    }

    return value;
}

/**
 * Generate evenly spaced points for the function given.
 * The algorithm is taken from the Project #3 Handout.
 *
 * @param num_pts the number of points
 * @param *func the function that you want to use to find y[i]
 * @param x[] the array that will hold the x[i] values
 * @param y[] the array that will hold the y[i] values
 */
void GenEvenPts (int num_pts, double(*func)(double), double x[], double y[])
{
    int i;
    double h = 10.0 / (double)(num_pts-1);
    double xtemp = -5.0;

    for (i=0; i<num_pts; i++)
    {
        x[i] = xtemp;
        y[i] = func(x[i]);
        xtemp += h;
    }
}

/**
 * Generate Chebychev points for the function given.
 * The algorithm is taken from the Project #3 Handout.
 *
 * @param num_pts the number of points
 * @param *func the function that you want to use to find y[i]
 * @param x[] the array that will hold the x[i] values
 * @param y[] the array that will hold the y[i] values
 */
void GenChebychevPts (int num_pts, double(*func)(double), double x[], double y[])
{
    int i;

    for (i=0; i<num_pts; i++)
    {
        x[i] = -5.0 * cos((float)i * M_PI / (float)(num_pts-1));
        y[i] = func(x[i]);
    }
}

int main (void)
{
    const int size = 11;
    double x[size];
    double y[size];
    double point = 2.0;

    //GenEvenPts (size, &func1, x, y);
    GenChebychevPts (size, &func1, x, y);

    for (int i=0; i<size; i++)
        printf("x[%d] = %e -- y[%d] = %e\n", i, x[i], i, y[i]);

    printf ("LaGrange = %e\n", LaGrangeMethod(x, y, size, point));
    printf ("Newton = %e\n", NewtonMethod(x, y, size, point));

    return 0;
}