6-7 Orthogonal Polynomials Approximation (40分)

Given a function f and a set of m>0 distinct points x​1​​ <x​2​​ <⋯<x​m​​ . You are supposed to write a function to approximate f by an orthogonal polynomial using the exact function values at the given m points with a weight w(x​i​​ ) assigned to each point x​i​​ . The total error err=∑​i=1​m​​ w(x​i​​ )[7f(x​i​​ )−P​n​​ (x​i​​ )]​2​​ must be no larger than a given tolerance.

Format of function:

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

where the function pointer double (*f)(double t) defines the function f; int m is the number of points; double x[] contains points x
​1​​ <x​2​​ <⋯<x​m​​ ; double w[] contains the values of a weight function at the given points x[]; double c[] contains the coefficients c​0​​ ,c​1​​ ,⋯,c​n​​ of the approximation polynomial P​n​​ (x)=c​0​​ +c
​1​​ x+⋯+c​n​​ x​n​​ ; double *eps is passed into the function as the tolerance for the error, and is supposed to be returned as the value of error. The function OPA is supposed to return the degree of the approximation polynomial.

Note: a constant Max_n is defined so that if the total error is still not small enough when n = Max_n, simply output the result obtained when n = Max_n.

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

思路

1.需要利用一个结构实现多项式相加，多项式内积。

最简单的方法就是使用数组，并且不需要记录最高次数，只需要将前面的项系数都变为0即可。相当于一个次数为MAX_n的多项式，前面某些项系数为0。

多项式加法：系数直接相加。

多项式乘x:将系数后移一位

多项式乘常数：系数直接乘

多项式内积：比较麻烦的是，变量可能是y也可能是x。一种思路是两个都先按照多项式算出$\sum _{k=0}^{max\_n}a[k]\cdot x_i^k$，之后根据传入的flag确定到底应该选择f(xi)还是$\sum _{k=0}^{max\_n}a[k]\cdot x_i^k$

2.ppt上计算正交多项式的时候用到了$\varphi (x)(x-B_1)$，很明显应该拆成$\varphi (x)x-B_1\varphi (x)$

3.内积别忘了乘w[i]

4.如果内积没算错，这题大概率就不会错。

上代码

double gx[MAX_m];
double gy[MAX_m];
double gw[MAX_m];
int gm;

double Polynomials_Sum(double* f,int flagf, double* g,int flagg)/*计算f(x)，g(y)的求和,flg表示后一个变量*/
{
    double ret = 0;
    for (int i = 0; i < gm; i++)
    {   
        double sumf = 0;
        double sumg = 0;
        for (int j = 0; j <= MAX_n; j++)
        {
            sumf +=  f[j] * pow(gx[i], j);
            sumg +=  g[j] * pow(gx[i], j);
        }   
        sumf = (flagf == 0) ? sumf : gy[i];
        sumg = (flagg == 0) ? sumg : gy[i];
        ret += gw[i] * sumg * sumf;
    }
    return ret;
}

int OPA(double (*f)(double t), int m, double x[], double w[], double c[], double* eps)
{
    for (int i = 0; i < m; i++)
    {
        gx[i] = x[i];
        gw[i] = w[i];
        gy[i] = f(x[i]);
        gm = m;
    }
    
    double phi0[MAX_n+1] = { 0 };
    double phi1[MAX_n+1] = { 0 };
    double phi2[MAX_n+1] = { 0 };
    phi0[0] = 1;

    double y[MAX_n+1] = {0};
    y[1] = 1;

    double a0=Polynomials_Sum(phi0, 0,y,1) / Polynomials_Sum(phi0,0, phi0, 0);
    
    for (int i = 0; i <= MAX_n; i++)
    {
        c[i] = phi0[i] * a0;
    }
    
    double error = Polynomials_Sum(y,1, y, 1) - a0 * Polynomials_Sum(phi0,0, y,1 );
    
    double xphi0[MAX_n+1] = {0};
    for (int i = 0; i < MAX_n; i++)
        xphi0[i+1] = phi0[i];

    double B1 = Polynomials_Sum(xphi0,0, phi0, 0) / Polynomials_Sum(phi0,0, phi0, 0);
    
    phi1[0] = -B1;
    phi1[1] = 1;

    a0 = Polynomials_Sum(phi1, 0, y, 1) / Polynomials_Sum(phi1,0, phi1, 0);
    for (int i = 0; i <= MAX_n; i++)
    {
        c[i] += a0 * phi1[i];
    }
    
    error -= a0 * Polynomials_Sum(phi1,0, y, 1);
    int k = 1;
    
    while (k < MAX_n && fabs(error) >= *eps)
    {
        k++;
        xphi0[0] = 0;
        for (int i = 0; i < MAX_n; i++)
            xphi0[i+1] = phi1[i];
        
        B1 = Polynomials_Sum(xphi0,0, phi1,0) / Polynomials_Sum(phi1,0, phi1, 0);
        
        double C1 = Polynomials_Sum(xphi0,0, phi0, 0) / Polynomials_Sum(phi0,0, phi0, 0);
        
        phi2[0] = 0;
        if (phi1[MAX_n] != 0) {
            *eps = error;
            return k;
        }
        for (int i = 0; i < MAX_n; i++)
            phi2[i + 1] = phi1[i];

        for (int i = 0; i <= MAX_n; i++)
            phi2[i] = phi2[i] - B1 * phi1[i] - C1 * phi0[i];
       
        double t1 = Polynomials_Sum(phi2, 0, y, 1);
        double t2 = Polynomials_Sum(phi2, 0, phi2, 0);
        a0 =  t1/t2 ;
        
        for (int i = 0; i <= MAX_n; i++)
            c[i] += a0 * phi2[i];
        
        error -= a0 * Polynomials_Sum(phi2,0, y, 1);
        for (int i = 0; i <= MAX_n; i++)
        {
            phi0[i] = phi1[i];
            phi1[i] = phi2[i];
        }    
    }
    *eps = error;
    return k;
}