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.

 * - Split the Newton algorithm into two parts: NewtonMethod() and NewtonCoef().

 *

 */

#include <cstdio>

#include <cmath>

using namespace std;

#define MAX_POINTS 26

// Global Variables for easier return values

double absmax, xsave;

/**

 * 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 Newton Coefficients.

 * This only needs to be called once per set of x,y values.

 *

 * @param x the x[i] values

 * @param y the y[i] values

 * @param a the a[i] values will be returned here

 * @param n the number of points

 */

void NewtonCoef (double *x, double *y, double *a, int n)

{

    int i, j;

    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]);

}

/**

 * 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 a the a[i] values (Newton Coefficients)

 * @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, double *a, int n, double point)

{

    int i;

    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]);

    }

}

/**

 * Brute-force search the Newton Interpolating Poly for the x and y values.

 * The error calculations will be for the given function (func1, func2, or func3).

 *

 * NOTE: The "return values" are stored in two global variables: absmax and xsave.

 *

 * @param x the x[i] values

 * @param y the y[i] values

 * @param n the number of points for this polynomial

 * @param func the function to use for the error calculation

 */

void brute_search_newton (double *x, double *y, int n, double(*func)(double))

{

    absmax = 0.0;

    double xval = -5.0;

    double h = 10.0/500.0;

    double temp, error;

    double a[n];

    int i;

    NewtonCoef(x, y, a, n);

    for (i=0; i<=500; i++)

    {

        temp = NewtonMethod(x, y, a, n, xval);

        error = abs(temp - func(xval));

        if (error > absmax)

        {

            absmax = error;

            xsave = xval;

        }

        xval += h;

    }

}

/**

 * Brute-force search the LaGrange Interpolating Poly for the x and y values.

 * The error calculations will be for the given function (func1, func2, or func3).

 *

 * NOTE: The "return values" are stored in two global variables: absmax and xsave.

 *

 * @param x the x[i] values

 * @param y the y[i] values

 * @param n the number of points for this polynomial

 * @param func the function to use for the error calculation

 */

void brute_search_lagrange (double *x, double *y, int n, double(*func)(double))

{

    absmax = 0.0;

    double xval = -5.0;

    double h = 10.0/500.0;

    double temp, error;

    int i;

    for (i=0; i<=500; i++)

    {

        temp = LaGrangeMethod(x, y, n, xval);

        error = abs(temp - func(xval));

        if (error > absmax)

        {

            absmax = error;

            xsave = xval;

        }

        xval += h;

    }

}

int main (void)

{

    int size, i;

    double x[MAX_POINTS];

    double y[MAX_POINTS];

    // In the following code, the correct functions are run for each of the

    // 12 times this needs to be run.

    printf ("Newton Polynomial, Equal Points, Function #1\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

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

        brute_search_newton (x, y, size, &func1);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("Newton Polynomial, Equal Points, Function #2\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenEvenPts (size, &func2, x, y);

        brute_search_newton (x, y, size, &func2);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("Newton Polynomial, Equal Points, Function #3\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenEvenPts (size, &func3, x, y);

        brute_search_newton (x, y, size, &func3);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Equal Points, Function #1\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

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

        brute_search_lagrange (x, y, size, &func1);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Equal Points, Function #2\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenEvenPts (size, &func2, x, y);

        brute_search_lagrange (x, y, size, &func2);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Equal Points, Function #3\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenEvenPts (size, &func3, x, y);

        brute_search_lagrange (x, y, size, &func3);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("Newton Polynomial, Chebychev Points, Function #1\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

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

        brute_search_newton (x, y, size, &func1);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("Newton Polynomial, Chebychev Points, Function #2\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenChebychevPts (size, &func2, x, y);

        brute_search_newton (x, y, size, &func2);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("Newton Polynomial, Chebychev Points, Function #3\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenChebychevPts (size, &func3, x, y);

        brute_search_newton (x, y, size, &func3);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Chebychev Points, Function #1\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

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

        brute_search_lagrange (x, y, size, &func1);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Chebychev Points, Function #2\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenChebychevPts (size, &func2, x, y);

        brute_search_lagrange (x, y, size, &func2);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    printf ("\n\n");

    printf ("LaGrange Polynomial, Chebychev Points, Function #3\n");

    printf ("============================================================\n");

    for (i=6; i<=26; i+=5)

    {

        size = i;

        GenChebychevPts (size, &func3, x, y);

        brute_search_lagrange (x, y, size, &func3);

        printf ("For %-2d Points, MAX ERROR is: %e"

                " -- OCCURS at x = %e\n", size, absmax, xsave);

    }

    return 0;

}