用最小二乘法進行任意階數的曲線拟合——C語言實作

/*********************************************
*        利用公式 ax=y--A'Ax=A'y 拟合曲線     *
*            用最小二乘法進行曲線拟合          *
*********************************************/

#include <math.h>
#include <stdio.h>
#define N 30
double power(double a[], double x, int i, int j);
void arryans(double a[][N], double b[N], double x[N], int n);
void arrytrans(double a[][N], double b[][N], int n, int m);
void arrymul(double a[][N], double b[][N], double c[][N], int m, int l, int n);
void arrymuls(double a[][N], double b[], double c[], int m, int n);
void arrycon(double a[], double b[][N], int n, int m);
void arryzeng(double a[][N], double b[], double c[][N], int m, int n);
int main(void)
{
    int i, n, m;
    double a[N][N], b[N][N], c[N][N], d[N], k[N], Sum = 0.0, Sum1;
    double x[N] = { 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0 };
    double y[N] = { 33.40, 79.50, 122.65, 159.05, 189.15,
        214.15, 238.65, 252.50, 267.55, 280.50,
        296.65, 301.40, 310.40, 318.15, 325.15 }; //初始資料

    printf("輸入要拟合的次數:");
    scanf("%d", &n);
    n++;                          //n是拟合的階數
    m = 15;                       //m是初始資料的長度
    arrycon(x, a, m, n);          //a[][]是A
    arrytrans(a, b, m, n);        //b[][]是A'
    arrymul(b, a, c, n, m, n);    //c[][]是A'*A
    arrymuls(b, y, d, n, m);      //d[][]是A'*y
    arryans(c, d, k, n);          //k[]是解，即最後的系數，按升幂排序
    for (i = 0; i < n; i++)
        printf("k[%d]=%f\n", i, k[i]);
    for (i = 0; i < m; i++) {
        Sum1 = power(k, x[i], 0, n);
        Sum += (y[i] - Sum1) * (y[i] - Sum1);
    }
    printf("拟合誤差：%lf\n", sqrt(Sum));
    return 0;
}
void arryans(double a[][N], double b[N], double x[N], int n) //解一個n階的線性方程組
{
    int i, j, k, mi;
    double mx, tmp;

    for (i = 0; i < n - 1; i++) {
        for (j = i + 1, mi = i, mx = fabs(a[i][i]); j < n; j++) //找主元素
            if (fabs(a[j][i]) > mx) {
                mi = j;
                mx = fabs(a[j][i]);
            }
        if (i < mi) { //交換兩行
            tmp = b[i];
            b[i] = b[mi];
            b[mi] = tmp;
            for (j = i; j < n; j++) {
                tmp = a[i][j];
                a[i][j] = a[mi][j];
                a[mi][j] = tmp;
            }
        }
        for (j = i + 1; j < n; j++) { //高斯消元
            tmp = -a[j][i] / a[i][i];
            b[j] += b[i] * tmp;
            for (k = i; k < n; k++)
                a[j][k] += a[i][k] * tmp;
        }
    }
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]; //求解方程
    for (i = n - 2; i >= 0; i--) {
        x[i] = b[i];
        for (j = i + 1; j < n; j++)
            x[i] -= a[i][j] * x[j];
        x[i] /= a[i][i];
    }
}
void arrytrans(double a[][N], double b[][N], int m, int n) //矩陣轉置,m行n列
{
    int i, j;

    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++)
            b[j][i] = a[i][j];
}
void arrymul(double a[][N], double b[][N], double c[][N], int m, int l, int n) //兩個矩陣的乘法
{
    int i, j, k; //c[m][n]=a[m][l]*b[l][n]

    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++) {
            c[i][j] = 0;
            for (k = 0; k < l; k++)
                c[i][j] += a[i][k] * b[k][j];
        }
}
void arrycon(double a[], double b[][N], int m, int n) //矩陣一維變二維
{
    int i, j, k;

    for (i = 0; i < m; i++)
        for (j = 0; j < n; j++) {
            b[i][j] = 1;
            for (k = 0; k < j; k++)
                b[i][j] *= a[i];
        }
}
void arryzeng(double a[][N], double b[], double c[][N], int m, int n) //增廣矩陣
{
    int i, j;

    for (i = 0; i < m; i++)
        for (j = 0; j <= n; j++) {
            if (j == n)
                c[i][j] = b[i];
            else
                c[i][j] = a[i][j];
        }
}
void arrymuls(double a[][N], double b[], double c[], int m, int n)
{ //由于一維數組按行存放，故乘列向量單獨寫
    int i, j;

    for (i = 0; i < m; i++) {
        c[i] = 0;
        for (j = 0; j < n; j++)
            c[i] += a[i][j] * b[j];
    }
}
double power(double a[], double x, int i, int j) //用秦九韶算法計算多項式
{
    if (i >= j)
        return a[i] * x;
    else
        return x * power(a, x, i + 1, j) + a[i];
}
           

該程式可以進行任意多項式拟合。

運作結果