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().
 * - Implemented brute_search_lagrange() and brute_search_newton() to find the
 *   worst values. This is taken from the Project #3 Handout.
 * - Added each needed call to main().
 *
 * 2005-11-17
 * - Changed stupid global variables used to return values from a function to
 *   a new type (called bsf_type) and use that instead. What was I thinking???
 *
 */

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

#define MAX_POINTS 26

// A convienent brute-force-search return type
typedef struct { double absmax, xsave; } bsf_type;

/**
 * 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).
 *
 * @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
 * @return a struct of the max abs error and the x value where it occurred
 */
bsf_type brute_search_newton (double *x, double *y, int n, double(*func)(double))
{
    bsf_type answer;
    answer.absmax = 0.0;
    answer.xsave = 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 > answer.absmax)
        {
            answer.absmax = error;
            answer.xsave = xval;
        }

        xval += h;
    }

    return answer;
}

/**
 * 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).
 *
 * @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
 * @return a struct of the max abs error and the x value where it occurred
 */
bsf_type brute_search_lagrange (double *x, double *y, int n, double(*func)(double))
{
    bsf_type answer;
    answer.absmax = 0.0;
    answer.xsave = 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 > answer.absmax)
        {
            answer.absmax = error;
            answer.xsave = xval;
        }

        xval += h;
    }

    return answer;
}

/**
 * Print a nice error line in the format we need.
 *
 * @param points the number of points
 * @param error a struct of the max abs error and the x value where it occurred
 */
void print_error_line (int points, bsf_type error)
{
    printf ("For %-2d Points, MAX ERROR is: %e"
            " -- OCCURS at x = %e\n", points, error.absmax, error.xsave);
}

int main (void)
{
    int i; // number of points
    double x[MAX_POINTS];
    double y[MAX_POINTS];
    bsf_type error;

    // 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)
    {
        GenEvenPts (i, &func1, x, y);
        error = brute_search_newton (x, y, i, &func1);
        print_error_line (i, error);

    }
    
    printf ("\n\n");
    printf ("Newton Polynomial, Equal Points, Function #2\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenEvenPts (i, &func2, x, y);
        error = brute_search_newton (x, y, i, &func2);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("Newton Polynomial, Equal Points, Function #3\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenEvenPts (i, &func3, x, y);
        error = brute_search_newton (x, y, i, &func3);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Equal Points, Function #1\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenEvenPts (i, &func1, x, y);
        error = brute_search_lagrange (x, y, i, &func1);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Equal Points, Function #2\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenEvenPts (i, &func2, x, y);
        error = brute_search_lagrange (x, y, i, &func2);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Equal Points, Function #3\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenEvenPts (i, &func3, x, y);
        error = brute_search_lagrange (x, y, i, &func3);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("Newton Polynomial, Chebychev Points, Function #1\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func1, x, y);
        error = brute_search_newton (x, y, i, &func1);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("Newton Polynomial, Chebychev Points, Function #2\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func2, x, y);
        error = brute_search_newton (x, y, i, &func2);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("Newton Polynomial, Chebychev Points, Function #3\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func3, x, y);
        error = brute_search_newton (x, y, i, &func3);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Chebychev Points, Function #1\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func1, x, y);
        error = brute_search_lagrange (x, y, i, &func1);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Chebychev Points, Function #2\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func2, x, y);
        error = brute_search_lagrange (x, y, i, &func2);
        print_error_line (i, error);
    }
    
    printf ("\n\n");
    printf ("LaGrange Polynomial, Chebychev Points, Function #3\n");
    printf ("============================================================\n");
    for (i=6; i<=26; i+=5)
    {
        GenChebychevPts (i, &func3, x, y);
        error = brute_search_lagrange (x, y, i, &func3);
        print_error_line (i, error);
    }
    
    return 0;
}