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;

}