数值计算方法-曲线拟合问题-最小二乘法

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1100;
double x[maxn], y[maxn];
double s[maxn], t[maxn];
double a[maxn][maxn], b[maxn];
double l[maxn][maxn], u[maxn][maxn];
void solve(int n)      //Doolitle分解法解线性方程组
{
    for(int i = 1; i <= n; i++)        l[i][i] = 1;
    for(int k = 1; k <= n; k++)
    {
        if(k == 1)
        {
            for(int j = 1; j <= n; j++) u[1][j] = a[1][j];
            for(int i = 2; i <= n; i++) l[i][1] = a[i][1] / u[1][1];
        }
        else
        {
            for(int j = k; j <= n; j++)
            {
                double tmp = 0;
                for(int s = 1; s <= k - 1; s++) tmp += l[k][s] * u[s][j];
                u[k][j] = a[k][j] - tmp;
            }
            for(int i = k + 1; i <= n; i++)
            {
                double tmp = 0;
                for(int s = 1; s <= k - 1; s++)                    tmp += l[i][s] * u[s][k];
                l[i][k] = (a[i][k] - tmp) / u[k][k];
            }
        }
    }
    y[1] = b[1];
    for(int k = 2; k <= n; k++)
    {
        double tmp = 0;
        for(int s = 1; s <= k - 1; s++)            tmp += l[k][s] * y[s];
        y[k] = b[k] - tmp;
    }
    x[n] = y[n] / u[n][n];
    for(int k = n - 1; k; k--)
    {
        double tmp = 0;
        for(int s = k + 1; s <= n; s++)            tmp += u[k][s] * x[s];
        x[k] = (y[k] - tmp) / u[k][k];
    }
    for(int i = 1; i <= n; i++) printf("x[%d]=%.6f\n", i, x[i]);
}
int main()
{
    cout << "请输入点的个数:\n";
    int n;
    cin >> n;
    cout << "请输入" << n << "个点对应的坐标值" << endl;
    for(int i = 1; i <= n; i++)
    {
        cin >> x[i] >> y[i];
    }
    cout << "请输入多项式系数：\n";
    int m;
    cin >> m;
    for(int i = 0; i <= 2 * m; i++)
    {
        double tmp = 0;
        for(int j = 1; j <= n; j++)           tmp = tmp + pow(x[j], i);
        s[i] = tmp;
    }
    for(int i = 0; i <= m; i++)
    {
        double tmp = 0;
        for(int j = 1; j <= n; j++)
        {
            tmp = tmp + pow(x[j], i) * y[j];
        }
        t[i] = tmp;
    }
    int tn = m + 1;
    for(int i = 1; i <= m + 1; i++)
    {
        int pos = i - 1;
        for(int j = 1; j <= m + 1; j++)
        {
            a[i][j] = s[pos + j - 1];
        }
    }
    for(int i = 1; i <= tn; i++)      b[i] = t[i - 1];
    solve(tn);
    return 0;
}/*请输入点的个数:9请输入9个点对应的坐标值0 00.5 3.21.2 4.21.6 4.82.2 6.12.6 6.93 8.23.6 9.94.2 12请输入多项式系数：3x[1]=0.421119x[2]=4.535287x[3]=-1.343178x[4]=0.221377*/经过写检验程序可得：
#include < bits / stdc++.h >
using namespace std;/*x[1]=0.421119x[2]=4.535287x[3]=-1.343178x[4]=0.221377*/
int main()
{
    double a[] = {0.421119, 4.535287, -1.343178, 0.221377};
    double xx[] = {0, 0, 0.5, 1.2, 1.6, 2.2, 2.6, 3, 3.6, 4.2};
    double y[] = {0, 0, 3.2, 4.2, 4.8, 6.1, 6.9, 8.2, 9.9, 12};
    for(int i = 1; i <= 10; i++)
    {
        double x = xx[i];
        double res = 0;
        for(int j = 0; j <= 3; j++)           res = res + a[j] * pow(x, j);
        cout << res << " " << y[i] << endl;
    }
    return 0;
}