實驗要求：

給定函數f(x)在 m 個采樣點處的值f()以及每個點的權重（1<= i <= m），求曲線擬合的正交多項式，滿足最小二乘誤差

。

函數接口定義：

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

在接口定義中：f是給定函數f(x)，m為采樣點個數，x傳入m個采樣點，w傳入每個點的權重；c傳回求得的曲線擬合的正交多項式的系數，eps傳入最小二乘誤差上限TOL，并且傳回擬合的實際最小二乘誤差

。

另外，在裁判程序中還定義了常數MAXN，當擬合多項式的階數達到MAXN卻仍未達到精度要求時，即返回 n = MAXN時的結果。

偽代碼：

---

?

? Algorithm: Orthogonal Polynomials Approximation

To approximate a given function by a polynomial with error bounded

by a given tolerance.

Input: number of data m; x[m]; y[m]; weight w[m]; tolerance TOL;

??????????? maximum degree of polynomial Max_n.

Output: coefficients of the approximating polynomial.

1? Set? j0(x) o 1;? a0 = (j0, y)/(j0, j0);? P(x) = a0 j0(x);? err = (y, y) - a0 (j0, y);

2? Set? a1= (xj0, j0)/(j0, j0);? j1(x) = (x - a1) j0(x);

1 = (j1, y)/(j1, j1);? P(x) += a1 j1(x);? err? -= a1 (j1, y);

Step 3? Set? k = 1;

4? While (( k < Max_n)&&(|err|3TOL))? do? steps 5-7

??????????????? Step 5? k ++;

6? ak= (xj1, j1)/(j1, j1);? bk-1 = (j1, j1)/(j0, j0);

2(x) = (x - ak) j1(x) - bk-1 j0(x);? ak = (j2, y)/(j2, j2);

(x) += ak j2(x);? err? -= ak (j2, y);

7? Set j0(x) = j1(x); j1(x) = j2(x);

Step 8? Output ( );? STOP.

---

c語言實現：#include<stdio.h> #include<math.h>#define MAXM 200 /* 采樣點個數上限*/ #define MAXN 5 /* 擬合多項式階數上限*//* 這里僅給出2個測試用例，正式裁判程序可能包含更多測試用例 */ double qi_x[MAXN+1][MAXN+1];//存儲基底函數 double ak[MAXN+1];//存儲擬合函數對應基地的系數double f1(double x) {return sin(x); }double f2(double x) {return exp(x); }//用于求x的n次方 double power(double x,int n) {int i;double y=1;for(i=0;i<n;i++){y=y*x;}return y; }//用于求解兩個基底函數的內積 double get_Inner_product(double x[],double y1[],double y2[] ,double w[],int m) {int i,j;double sum=0;double sum1,sum2;for(i=0;i<m;i++){sum1=0;sum2=0;for(j=0;j<=MAXN;j++){sum1=sum1+y1[j]*power(x[i],j);sum2=sum2+y2[j]*power(x[i],j);}sum=sum+w[i]*sum1*sum2;}return sum; }//用于完成x*a(x)的功能void carry_x(double x[],double *temp) {int i;for(i=MAXN+1;i>0;i--)temp[i]=x[i-1];temp[0]=0; }//用于完成擬合函數和基底函數的內積 double get_function_inner_product(double (*f)(double t),double x[],double y1[],double w[],int m) {int i,j;double sum=0;double sum1;for(i=0;i<m;i++){sum1=0;for(j=0;j<=MAXN;j++){sum1=sum1+y1[j]*power(x[i],j);}sum=sum+w[i]*f(x[i])*sum1;}// printf("%f\n",sum);return sum; }//求解正交基底 void Array_op(int n,double fg,double fg1) {int i;carry_x(qi_x[n-1],qi_x[n]);if(n==1){for(i=0;i<=n;i++)qi_x[n][i]=qi_x[n][i]-fg*qi_x[n-1][i];}else{for(i=0;i<=n;i++){qi_x[n][i]=qi_x[n][i]-fg*qi_x[n-1][i]-fg1*qi_x[n-2][i];}//qi_x[n][0]=-fg*qi_x[n-1][0]-fg1*qi_x[n-2][0];} }//完成函數擬合并返回誤差和擬合次數 int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps ) {double err;double sum=0;double m1,n;double temp[MAXN+1]={0};int i;qi_x[0][0]=1;ak[0]=get_function_inner_product(f,x,qi_x[0],w,m)/get_Inner_product(x,qi_x[0],qi_x[0],w,m);for(i=0;i<m;i++){sum=sum+w[i]*f(x[i])*f(x[i]);}err=sum-ak[0]*get_function_inner_product(f,x,qi_x[0],w,m);// printf("%f\n",err);carry_x(qi_x[0],temp);m1=get_Inner_product(x,temp,qi_x[0],w,m)/get_Inner_product(x,qi_x[0],qi_x[0],w,m);Array_op(1,m1,0.0);ak[1]=get_function_inner_product(f,x,qi_x[1],w,m)/get_Inner_product(x,qi_x[1],qi_x[1],w,m);err=err-ak[1]*get_function_inner_product(f,x,qi_x[1],w,m);int k=1;while((k<MAXN)&&(fabs(err)>=*eps)){k=k+1;double temp1[MAXN+1]={0};carry_x(qi_x[k-1],temp1);m1=get_Inner_product(x,temp1,qi_x[k-1],w,m)/get_Inner_product(x,qi_x[k-1],qi_x[k-1],w,m);n=get_Inner_product(x,qi_x[k-1],qi_x[k-1],w,m)/get_Inner_product(x,qi_x[k-2],qi_x[k-2],w,m);Array_op(k,m1,n);// printf("%f\n",get_function_inner_product(f,x,qi_x[k],w,m));ak[k]=get_function_inner_product(f,x,qi_x[k],w,m)/get_Inner_product(x,qi_x[k],qi_x[k],w,m);// printf("%f\n",ak[k]);err=err-ak[k]*get_function_inner_product(f,x,qi_x[k],w,m);}*eps=err;return k; }//輸出實驗結果 void print_results( int n, double c[], double eps) { /* 用于輸出結果 */int i;int m,k;for(m=0;m<MAXN+1;m++){c[m]=0;for(k=0;k<MAXN+1;k++)c[m]=c[m]+ak[k]*qi_x[k][m];}printf("%d\n", n);for (i=0; i<=n; i++)printf("%8.4e ", c[i]);printf("\n");printf("error = %6.2e\n", eps);printf("\n"); }int main() {int m, i, n;double x[MAXM], w[MAXM], c[MAXN+1], eps;/* 測試1 */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);/* 測試2 */memset(qi_x,0,sizeof(qi_x));memset(ak,0,sizeof(ak));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; } --------------------- 作者：couchy-wu 來源：CSDN 原文：https://blog.csdn.net/qq_41989372/article/details/84638914 版權聲明：本文為博主原創文章，轉載請附上博文鏈接！

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

?

總結

以上是生活随笔為你收集整理的正交多项式拟合函数的全部內容，希望文章能夠幫你解決所遇到的問題。

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