NA-PTA-NA Lab7 正交多项式拟合的免病态方法

Orthogonal Polynomials Approximation

Format of function:

int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps );

Sample program of judge:

#include <stdio.h>
#include <math.h>

#define MAX_m 200
#define MAX_n 5

double f1(double x)
{
    return sin(x);
}

double f2(double x)
{
    return exp(x);
}

int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps );

void print_results( int n, double c[], double eps)
{    
    int i;

    printf("%d\n", n);
    for (i=0; i<=n; i++)
        printf("%12.4e ", c[i]);
    printf("\n");
    printf("error = %9.2e\n", eps);
    printf("\n");
}

int main()
{
    int m, i, n;
    double x[MAX_m], w[MAX_m], c[MAX_n+1], eps;

    m = 90;
    for (i=0; i<m; i++) {
        x[i] = 3.1415926535897932 * (double)(i+1) / 180.0;
        w[i] = 1.0;
    }
    eps = 0.001;
    n = OPA(f1, m, x, w, c, &eps);
    print_results(n, c, eps);

    m = 200;
    for (i=0; i<m; i++) {
        x[i] = 0.01*(double)i;
        w[i] = 1.0;
    }
    eps = 0.001;
    n = OPA(f2, m, x, w, c, &eps);
    print_results(n, c, eps);

    return 0;
}

/* Your function will be put here */

Sample Output:

3
 -2.5301e-03   1.0287e+00  -7.2279e-02  -1.1287e-01 
error =  6.33e-05

4
  1.0025e+00   9.6180e-01   6.2900e-01   7.0907e-03   1.1792e-01 
error =  1.62e-04

My Answer:
这个代码还有两个点没过，但是不知道为什么。

#include <string.h>
#include <stdlib.h>
int M;
double W[MAX_m];

double myPow(double x, int n)
{
	double res = 1.0;
	while (n--) res *= x;
	return res;
}
double CalculatePoly(double poly[], double x)
{
	double res = 0;
	for (int i = 0; i < MAX_n + 1; i++)
		res += poly[i] * myPow(x, i);
	return res;
}

double InnerProd_yy(double y[])
{	// y直接提供第一二个基函数值
	double sum = 0;
	for (int i = 0; i < M; i++) 
		sum += W[i] * y[i] * y[i];
	return sum;
}
double InnerProd_y(double poly1[], double x[], double y[])
{	// f的系数在x[]求值提供第一个基函数值  y直接提供第二个基函数值
	double sum = 0;
	for (int i = 0; i < M; i++) 
		sum += W[i] * CalculatePoly(poly1, x[i]) * y[i];
	return sum;
}
double InnerProd(double poly1[], double poly2[], double x[])
{	// f的系数在x[]求值提供第一个基函数值  g的系数在x[]求值提供第二个基函数值
	double sum = 0;
	for (int i = 0; i < M; i++) 
		sum += W[i] * CalculatePoly(poly1, x[i]) * CalculatePoly(poly2, x[i]);
	return sum;
}
 poly deal /
double* x_Mul_Poly(double poly[])
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));

	for (int i = MAX_n; i >= 1; i--) 
		new_poly[i] = poly[i - 1];
	new_poly[0] = 0.0;

	return new_poly;
}
double* Poly_Sub_Val(double poly[], double val)
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));
	
	memcpy(new_poly, poly, (MAX_n + 1) * sizeof(double));
	new_poly[0] -= val;
	
	return new_poly;
}
double* Val_Mul_Poly(double val, double poly[])
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));

	for (int i = 0; i <= MAX_n; i++)
		new_poly[i] = val * poly[i];

	return new_poly;
}
double* Get_x_poly()
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));

	memset(new_poly, 0, (MAX_n + 1) * sizeof(double));
	new_poly[1] = 1.0;

	return new_poly;
}
double* Poly_Add_Poly(double* poly1, double* poly2)
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));

	for (int i = 0; i <= MAX_n; i++)
		new_poly[i] = poly1[i] + poly2[i];

	return new_poly;
}
double* Poly_Sub_Poly(double* poly1, double* poly2)
{
	double* new_poly = (double*)malloc((MAX_n + 1) * sizeof(double));

	for (int i = 0; i <= MAX_n; i++)
		new_poly[i] = poly1[i] - poly2[i];

	return new_poly;
}

 opa ///

int OPA(double(*f)(double t), int m, double x[], double w[], double c[], double *eps)
{
		M = m;
		memcpy(W, w, MAX_m * sizeof(double));

		double y[MAX_m];
		for (int i = 0; i < m; i++) y[i] = f(x[i]);

		double phi[MAX_n + 1][MAX_n + 1] = { 0 };		// 多项式升序排列 项次数从0到MAX_n  phi[MAX_n]也是需要的 因为phi_k也要计算
		double a[MAX_n + 3] = { 0 };					// 插值多项式phi(x)的系数
		double B[MAX_n + 3] = { 0 };
		double C[MAX_n + 3] = { 0 };
		double err = InnerProd_yy(y);

		// initialize   以B C phi a err的顺序
		phi[0][0] = 1;
		a[0] = InnerProd_y(phi[0], x, y) / InnerProd(phi[0], phi[0], x);
		err -= a[0] * InnerProd_y(phi[0], x, y);

		B[1] = InnerProd(x_Mul_Poly(phi[0]), phi[0], x) / InnerProd(phi[0], phi[0], x);
		memcpy(phi[1], Poly_Sub_Val(Get_x_poly(), B[1]), (MAX_n + 1) * sizeof(double));
		a[1] = InnerProd_y(phi[1], x, y) / InnerProd(phi[1], phi[1], x);
		err -= a[1] * InnerProd_y(phi[1], x, y);

		int k;
		for (k = 2; k <= MAX_n; k++) {
			B[k] = InnerProd(x_Mul_Poly(phi[k - 1]), phi[k - 1], x) / InnerProd(phi[k - 1], phi[k - 1], x);
			C[k] = InnerProd(x_Mul_Poly(phi[k - 1]), phi[k - 2], x) / InnerProd(phi[k - 2], phi[k - 2], x);
			memcpy(phi[k], Poly_Sub_Poly(
				Poly_Sub_Poly(x_Mul_Poly(phi[k - 1]), Val_Mul_Poly(B[k], phi[k - 1])),
				Val_Mul_Poly(C[k], phi[k - 2])),
				(MAX_n + 1) * sizeof(double));
			a[k] = InnerProd_y(phi[k], x, y) / InnerProd(phi[k], phi[k], x);
			err -= a[k] * InnerProd_y(phi[k], x, y);

			if (fabs(err) < *eps) break;
		}
		int deg = k;	// res_poly的次数是0-k

		double* res_poly = (double*)malloc((MAX_n + 1) * sizeof(double));
		memset(res_poly, 0, (MAX_n + 1) * sizeof(double));
		for (int k = 0; k <= deg; k++) {
			res_poly = Poly_Add_Poly(res_poly, Val_Mul_Poly(a[k], phi[k]));
		}

		// store the res_poly[] in c[]
		memcpy(c, res_poly, (MAX_n + 1) * sizeof(double));

		*eps = err;

		return deg;
}