最小二乘多項式擬合僅適用於偶數個座標

Question

我已經編寫了一個程序，該程序通過給出多個座標來計算拋物線的方程，但它僅適用於偶數個座標。如果我輸入一個不規則的數字，它會顯示一些無意義的bs。

下面是代碼（我不能在這裏複製，因爲一些格式問題的代碼）：

#include<iostream> 
#include<iomanip> 
#include<cmath> 
using namespace std; 
int main() 
{ 
    int i, j, k, n, N; 
    n = 2; 
    cout << "Number of data pairs:" << endl; 
    cin >> N; 
    double * x = new double[N]; 
    double * y = new double[N]; 
    cout << endl << "Enter the x-axis values:" << endl; 
    for (i = 0; i<N; i++) 
     cin >> x[i]; 
    cout << endl << "Enter the y-axis values:" << endl; 
    for (i = 0; i<N; i++) 
     cin >> y[i]; 
    double X[5]; 
    for (i = 0; i<2 * n + 1; i++) 
    { 
     X[i] = 0; 
     for (j = 0; j<N; j++) 
      X[i] = X[i] + pow(x[j], i); 
    } 
    double B[3][4], a[3]; //B is the Normal matrix(augmented) that will store the equations, 'a' is for value of the final coefficients 
    for (i = 0; i <= n; i++) 
     for (j = 0; j <= n; j++) 
      B[i][j] = X[i + j]; 
    double Y[3]; 
    for (i = 0; i<n + 1; i++) 
    { 
     Y[i] = 0; 
     for (j = 0; j<N; j++) 
      Y[i] = Y[i] + pow(x[j], i)*y[j]; 
    } 
    for (i = 0; i <= n; i++) 
     B[i][n + 1] = Y[i]; 
    n = n + 1; 
    for (i = 0; i < n; i++) 
    { 
     for (k = i + 1; k < n; k++) 
      if (B[i][i] < B[k][i]) 
       for (j = 0; j <= n; j++) 
       { 
        double temp = B[i][j]; 
        B[i][j] = B[k][j]; 
        B[k][j] = temp; 
       } 
    } 
     for (i = 0; i<n - 1; i++) //loop to perform the gauss elimination 
      for (k = i + 1; k<n; k++) 
      { 
       double t = B[k][i]/B[i][i]; 
       for (j = 0; j <= n; j++) 
        B[k][j] = B[k][j] - t*B[i][j]; //make the elements below the pivot elements equal to zero or elimnate the variables 
      } 
     for (i = n - 1; i >= 0; i--) //back-substitution 
     { 
      a[i] = B[i][n]; 
      for (j = 0; j < n; j++) 
       if (j != i) 
        a[i] = a[i] - B[i][j] * a[j]; 
      a[i] = a[i]/B[i][i]; 
     } 
     cout << endl << "The equation is the following:" << endl;; 
     cout << a[2] << "x^2 + " << a[1] << "x + " << a[0]; 
     cout << endl; 
     delete[]x; 
     delete[]y; 
     system("pause"); 
     return 0; 
} 

我從3得到的輸出和4個座標：

3座標：

數據點數：

輸入的x軸的值：

輸入y軸值

方程式如下：

1.02762e + 47X^2 + -4.316e + 47X + 3.90495e + 47個

4座標：

數據點數：

輸入的x軸的值：

輸入y軸值

方程式如下：

1X^2 + -0x + 0

任何意見或建議嗎？

在此先感謝

Answer 1

的意見已經發布請注意，某些索引超出範圍（超出矩陣的大小）。我寧願輸入x，y數據作爲每行一對值，一個x值和一個y值，這是一個很容易的變化（cin >> x [i] >> y [i]），但是以下示例使用原始問題的順序。用於二次方程擬合常規方法

示例代碼：

#include<iostream> 
#include<iomanip> 
#include<cmath> 
using namespace std; 
int main() 
{ 
    int i, j, k, N; 
    cout << "Number of data pairs:" << endl; 
    cin >> N; 
    double * x = new double[N]; 
    double * y = new double[N]; 
    cout << endl << "Enter the x-axis values:" << endl; 
    for (i = 0; i<N; i++) 
     cin >> x[i]; 
    cout << endl << "Enter the y-axis values:" << endl; 
    for (i = 0; i<N; i++) 
     cin >> y[i]; 
    double B[3][4] = {0.0};   // generate augmented matrix 
    for(k = 0; k < 3; k++){ 
     for (i = 0; i < N; i++) { 
      for (j = 0; j < 3; j++) { 
       B[k][j] += pow(x[i], j + k);} 
      B[k][3] += y[i]*pow(x[i], k);}} 
    for(k = 0; k < 3; k++){   // invert matrix 
     double q = B[k][k];   // divide row by B[k][k] 
     for(i = 0; i < 4; i++){ 
      B[k][i] /= q;} 
     for(j = 0; j < 3; j++){  // zero out column B[][k] 
      if(j == k) 
       continue; 
      double m = B[j][k];   
      for(i = 0; i < 4; i++){ 
       B[j][i] -= m*B[k][i];}}} 
    cout << endl << "The equation is the following:" << endl;; 
    cout << B[2][3] << " x^2 + " << B[1][3] << " x + " << B[0][3]; 
    cout << endl; 
    delete[]x; 
    delete[]y; 
    system("pause"); 
    return 0; 
} 

實施例用於通用度方程常規代碼：

#include<iostream> 
#include<iomanip> 
#include<cmath> 
using namespace std; 
int main() 
{ 
    int d, i, j, k, n; 
    cout << "Degree of equation:" << endl; 
    cin >> d; 
    double **A = new double *[d+1]; 
    for (k = 0; k < d+1; k++) { 
     A[k] = new double[d+2]; 
     for (i = 0; i < d+2; i++) { 
      A[k][i] = 0.0;}} 
    cout << "Number of data pairs:" << endl; 
    cin >> n; 
    double * x = new double[n]; 
    double * y = new double[n]; 
    cout << endl << "Enter the x-axis values:" << endl; 
    for (i = 0; i < n; i++) 
     cin >> x[i]; 
    cout << endl << "Enter the y-axis values:" << endl; 
    for (i = 0; i < n; i++) 
     cin >> y[i]; 
    for(k = 0; k < d+1; k++){ 
     for (i = 0; i < n; i++) { 
      for (j = 0; j < d+1; j++) { 
       A[k][j] += pow(x[i], j + k);} 
      A[k][d+1] += y[i]*pow(x[i], k);}} 
    for(k = 0; k < d+1; k++){  // invert matrix 
     double q = A[k][k];   // divide A[k][] by A[k][k] 
            // if q == 0, would need to swap rows         
     for(i = 0; i < d+2; i++){ 
      A[k][i] /= q;} 
     for(j = 0; j < d+1; j++){ // zero out column A[][k] 
      if(j == k) 
       continue; 
      double m = A[j][k];   
      for(i = 0; i < d+2; i++){ 
       A[j][i] -= m*A[k][i];}}} 
    cout << endl << "The equation is the following:" << endl; 
    for(k = d; k >= 2; k--) 
     cout << A[k][d+1] << " x^" << k << " + "; 
    cout << A[1][d+1] << " x" << " + " << A[0][d+1] << endl; 
    for (k = 0; k < d+1; k++) 
     delete[] A[k]; 
    delete[]A; 
    delete[]y; 
    delete[]x; 
    system("pause"); 
    return 0; 
} 

樣品測試輸入：

3 
4 
1 
2 
3 
4 
10 
49 
142 
313 

我通常使用一種替代算法，避免必須反轉矩陣，這是更好的更高階的多項式，但不需要二次方程。

如果對替代算法感興趣，下面是一個pdf文件的鏈接，該文件描述了算法幷包含僞代碼。我有舊的工作代碼，但它需要轉換（它使用cgets（），這是很少支持了））。

http://rcgldr.net/misc/opls.pdf

實施例的代碼。它真的很舊，最初是在Atari ST上運行的。我將它轉換爲與Visual Studio一起工作。所有靜態的原因是爲了減少鏈接器必須處理的符號數量。

/*------------------------------------------------------*/ 
/*  fit4.c   polynomial fit program   */ 
/*  originally written in 1990, minimal updates  */ 
/*------------------------------------------------------*/ 
/* disable Visual Studio warnings for old functions */ 
#define _CRT_SECURE_NO_WARNINGS 1 
#include <stdio.h> 

/*        mmax = max # coefficients */ 
/*        nmax = max # points */ 
/*        bfrsz = bfr size */ 
/*        linsz = size of line bfr */ 
#define mmax 11 
#define nmax 300 
#define bfrsz 0x2000 
#define linsz 64 

static void polyf(); 
static void gcoef(); 
static void gvar(); 
static void calc(); 
static void calcx(); 
static void rdata(); 
static void gdata(); 
static int gtlin(); 
static char gtchr(); 
static int conrs(); 

static char cbfr[64];   /* console response bfr */ 
static char line[linsz]; 
static int lineno; 

static char sbfr[bfrsz];  /* file params */ 
static FILE *sfp; 
static char *sptr, *send; 
static int gteof; 

static int wf; 

static int m, n;    /* input values */ 
static double x[nmax]; 
static double y[nmax]; 
static double w[nmax]; 
static double b[mmax];   /* generated values */ 
static double A[mmax]; 
static double B[mmax]; 
static double L[mmax]; 
static double W[mmax]; 
static double p2[nmax]; 
static double p1[nmax]; 
static double p0[nmax]; 
static double c[mmax];   /* coefficients for y(x) */ 
static double z[nmax];   /* calculated y[] */ 
static double xi, zi, vr; 
static double D0, D1;   /* constants */ 

static double *pp2;    /* pointers to logical p2, p1, p0 */ 
static double *pp1; 
static double *pp0; 
static double *ppx; 

static double *px, *pf, *pw; /* for gdata */ 

main() 
{ 
int i; 

name0: 
    printf("\nEnter name of data file: "); /* get file name */ 
    if(!conrs()) 
     return(0); 
    sfp = fopen(&cbfr[2], "rb"); /* open file */ 
    if(sfp == (FILE *)0){ 
     printf("\nfile not found"); 
     goto name0;} 

    wf = 0;      /* ask for weighting */ 
    printf("\nUsing weights (Y/N)? "); 
    if(!conrs()) 
     return(0); 
    if('Y' == (cbfr[2]&0x5f)) 
     wf = 1; 

deg0: 
    printf("\nEnter degree of equation (1-n): "); /* get # terms */ 
    if(!conrs()) 
     return(0); 
    sscanf(&cbfr[2], "%d", &m); 
    if(m >= mmax) 
     goto deg0; 

    sptr = send = (char *)0; 
    gteof = 0; 
    lineno = 0; 

    gdata();     /* get data */ 

    printf("\n%5d points found ", n); 

    polyf();     /* generate b[], A[], B[] */ 
    gcoef();     /* generate coefficients */ 
    gvar();      /* generate variance */ 

    for(i = 0; i <= m; i++){ 
     printf("\n%3d b%12.4le A%12.4le B%12.4le", 
       i, b[i], A[i], B[i]); 
     printf(" L%12.4le W%12.4le", L[i], W[i]);} 

    printf("\nvariance = %12.4le\n", vr); 

    for(i = m; i; i--) 
     printf("%12.4le X**%1d + ", c[i], i); 
    printf("%12.4le\n", c[0]); 

    for(i = 0; i < n; i++){  /* calculate results */ 
     xi = x[i]; 
     calc(); 
     z[i] = zi;} 

    for(i = 0; i < n; i += 1) /* display results */ 
     printf("\n%14.6le %14.6le %14.6le %14.6le", 
       x[i], y[i], z[i], y[i]-z[i]); 
    printf("\n"); 
    return(0); 
} 

/*------------------------------------------------------*/ 
/*  polyf   poly fit      */ 
/*  in:    x[], y[], w[], n, m    */ 
/*  out:   b[], A[], B[], L[], W[]   */ 
/*------------------------------------------------------*/ 
static void polyf() 
{ 
int i, j; 

    D0 = (double)0.;   /* init */ 
    D1 = (double)1.; 
    pp2 = p2; 
    pp1 = p1; 
    pp0 = p0; 

    j = 0; 
    A[j] = D0;     /* calc A, p[j], p[j-1], L, W */ 
    L[j] = D0;     /* note A[0] not used */ 
    W[j] = D0; 
    for(i = 0; i < n; i++){ 
     pp0[i] = D1; 
     pp1[i] = D0; 
     L[j] += w[i]; 
     W[j] += w[i]*y[i];} 
    B[0] = D0; 
    b[j] = W[j]/L[j]; 

    for(j = 1; j <= m; j++){ 
     ppx = pp2;    /* save old p[j], p[j-1] */ 
     pp2 = pp1; 
     pp1 = pp0; 
     pp0 = ppx; 
     A[j] = D0;    /* calc A */ 
     for(i = 0; i < n; i++){ 
      A[j] += w[i]*x[i]*pp1[i]*pp1[i]/L[j-1];} 
     L[j] = D0;    /* calc p[j], L, W */ 
     W[j] = D0; 
     for(i = 0; i < n; i++){ 
      pp0[i] = (x[i]-A[j])*pp1[i]-B[j-1]*pp2[i]; 
      L[j] += w[i]*pp0[i]*pp0[i]; 
      W[j] += w[i]*y[i]*pp0[i];} 
     B[j] = L[j]/L[j-1];  /* calc B[], b[] */ 
     b[j] = W[j]/L[j];} 
} 

/*------------------------------------------------------*/ 
/*  gcoef generate coefficients     */ 
/*  in:  b[], A[], B[]       */ 
/*  out: c[]          */ 
/*  uses: p0[], p1[], p2[]      */ 
/*------------------------------------------------------*/ 
static void gcoef() 
{ 
int i, j; 
    for(i = 0; i <= m; i++){   /* init */ 
     c[i] = p2[i] = p1[i] = p0[i] = 0.;} 
    p0[0] = D1; 
    c[0] += b[0]*p0[0]; 
    for(j = 1; j <= m; j++){   /* generate coefs */ 
     p2[0] = p1[0]; 
     p1[0] = p0[0]; 
     p0[0] = -A[j]*p1[0]-B[j-1]*p2[0]; 
     c[0] += b[j]*p0[0]; 
     for(i = 1; i <= j; i++){ 
      p2[i] = p1[i]; 
      p1[i] = p0[i]; 
      p0[i] = p1[i-1]-A[j]*p1[i]-B[j-1]*p2[i]; 
      c[i] += b[j]*p0[i];}} 
} 

/*------------------------------------------------------*/ 
/*  gvar generate variance      */ 
/*------------------------------------------------------*/ 
static void gvar() 
{ 
int i; 
double tt; 
    vr = 0.; 
    for(i = 0; i < n; i++){ 
     xi = x[i]; 
     calc(); 
     tt = y[i]-zi; 
     vr += tt*tt;} 
    vr /= n-m-1; 
} 

/*------------------------------------------------------*/ 
/*  calc   calc zi, given xi    */ 
/*  in:    c[]        */ 
/*------------------------------------------------------*/ 
static void calc() 
{ 
int i; 
    zi = c[m]; 
    for(i = m-1; i >= 0; i--) 
     zi = zi*xi + c[i]; 
} 

/*------------------------------------------------------*/ 
/*  calcx   calc zi, given xi    */ 
/*  in:    b[], A[], B[]     */ 
/*------------------------------------------------------*/ 
static void calcx() 
{ 
int i; 
double q2, q1, q0; 
    if(m == 0){ 
     zi = b[0]; 
     return;} 
    if(m == 1){ 
     zi = b[0]+(xi-A[1])*b[1]; 
     return;} 
    q1 = b[m]; 
    q0 = b[m-1]+(xi-A[m])*q1; 
    for(i = m-2; i >= 0; i--){ 
     q2 = q1; 
     q1 = q0; 
     q0 = b[i]+(xi-A[i+1])*q1-B[i+1]*q2;} 
    zi = q0; 
} 

/*------------------------------------------------------*/ 
/*  gdata get data        */ 
/*------------------------------------------------------*/ 
static void gdata() 
{ 
    px = &x[0]; 
    pf = &y[0]; 
    pw = &w[0]; 
    while(1){ 
     gtlin();    /* get a line */ 
     if(gteof) 
      break; 
     if(lineno == nmax){ 
      printf("\ntoo many points\n"); 
      break;} 
     if(wf){     /* stuff values */ 
      sscanf(line, "%le%le%le", px, pf, pw);} 
     else{ 
      sscanf(line, "%le%le", px, pf); 
      *pw = 1.0;} 
     px++;     /* bump ptrs */ 
     pf++; 
     pw++;} 
    fclose(sfp);    /* close file */ 
    n = lineno;     /* set # points */ 
} 

/*------------------------------------------------------*/ 
/*  gtlin get a line of data      */ 
/*------------------------------------------------------*/ 
static int gtlin() 
{ 
char chr; 
int col; 
    col = 0; 
    while(1){ 
     chr = gtchr(); 
     switch(chr){ 
      case 0x0a:   /* line feed */ 
      lineno++; 
      return(col); 
      case 0x1a: 
      return(col); 
      default: 
      line[col] = chr; 
      col++; 
      if(col >= linsz){ 
       printf("line # %d too long\n%s",lineno, line); 
       return(col);}}} 
} 

/*------------------------------------------------------*/ 
/*  gtchr get a char        */ 
/*------------------------------------------------------*/ 
static char gtchr() 
{ 
int cnt; 
    if(gteof)     /* check for eof */ 
     return(0x1a); 
    if(sptr == send){ 
     if(!(cnt = (int) fread(sbfr, 1, bfrsz, sfp))){ 
      fclose(sfp); 
      gteof = 1; 
      return(0x1a);} 
     sptr = sbfr; 
     send = sbfr+cnt;} 
    return(*sptr++); 
} 

/*------------------------------------------------------*/ 
/*  conrs   get string from console   */ 
/*------------------------------------------------------*/ 
static int conrs() 
{ 
int i; 
    memset(cbfr, 0, sizeof(cbfr)); /* get a line */ 
    cbfr[0] = sizeof(cbfr)-2; 
    fgets(cbfr+2, sizeof(cbfr)-2, stdin); 
    cbfr[1] = (char)(strlen(&cbfr[2])-1); 
    i = cbfr[1]; 
    cbfr[2+i] = 0; 
    return(i); 
} 

示例數據文件。我將它命名爲fitdat.txt：

1 1 
2 4 
3 9 

編譯並運行該程序。輸入數據文件的名稱，然後在提示使用權重時輸入N，然後輸入2以生成等式的程度。 （如果使用權重，第一列將是權重因子，權重2將是相同的兩個相同數據點的實例，但權重可以是1.5的值）。

Answer 2

在這裏，你正在閱讀的初始化值a[j]：

for (i = n - 1; i >= 0; i--) //back-substitution 
{ 
    a[i] = B[i][n]; 
    for (j = 0; j < n; j++) 
     if (j != i) 
      a[i] = a[i] - B[i][j] * a[j]; 
    a[i] = a[i]/B[i][i]; 
}