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 5double 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]?x_i^k$，之后根據傳入的flag確定到底應該選擇f(xi)還是$\sum _{k=0}^{max\_n}a[k]?x_i^k$

2.ppt上計算正交多項式的時候用到了$?(x)(x?B_1)$，很明顯應該拆成$?(x)x?B_1?(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; }

總結

以上是生活随笔為你收集整理的6-7 Orthogonal Polynomials Approximation (40分)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。