C++最小二乘法多项式拟合
时间: 2023-09-04 14:15:29 浏览: 51
最小二乘法多项式拟合是一种常见的数据拟合方法。它的基本思路是:假设有一组数据 $(x_1, y_1), (x_2, y_2), \cdots, (x_n, y_n)$,我们要用一个 $m$ 次多项式 $f(x)=a_0+a_1x+a_2x^2+\cdots+a_mx^m$ 来拟合这些数据,使得拟合出的多项式 $f(x)$ 在各个数据点上的误差平方和最小。
假设我们要拟合的多项式是 $f(x)=a_0+a_1x+a_2x^2+\cdots+a_mx^m$,我们可以将每个数据点 $(x_i, y_i)$ 代入多项式中,得到:
$$
f(x_i)=a_0+a_1x_i+a_2x_i^2+\cdots+a_mx_i^m
$$
我们希望拟合出的多项式 $f(x)$ 在各个数据点上的误差平方和最小,因此我们需要最小化下面的误差平方和:
$$
S=\sum_{i=1}^n(f(x_i)-y_i)^2=\sum_{i=1}^n(a_0+a_1x_i+a_2x_i^2+\cdots+a_mx_i^m-y_i)^2
$$
为了求解 $a_0, a_1, \cdots, a_m$,我们需要对 $S$ 求偏导数,并令其为 0:
$$
\begin{aligned}
\frac{\partial S}{\partial a_0}&=2\sum_{i=1}^n(a_0+a_1x_i+a_2x_i^2+\cdots+a_mx_i^m-y_i)=0 \\
\frac{\partial S}{\partial a_1}&=2\sum_{i=1}^n(a_0+a_1x_i+a_2x_i^2+\cdots+a_mx_i^m-y_i)x_i=0 \\
\cdots \\
\frac{\partial S}{\partial a_m}&=2\sum_{i=1}^n(a_0+a_1x_i+a_2x_i^2+\cdots+a_mx_i^m-y_i)x_i^m=0 \\
\end{aligned}
$$
将上面的方程组写成矩阵形式:
$$
\begin{bmatrix}
n & \sum_{i=1}^n x_i & \cdots & \sum_{i=1}^n x_i^m \\
\sum_{i=1}^n x_i & \sum_{i=1}^n x_i^2 & \cdots & \sum_{i=1}^n x_i^{m+1} \\
\vdots & \vdots & \ddots & \vdots \\
\sum_{i=1}^n x_i^m & \sum_{i=1}^n x_i^{m+1} & \cdots & \sum_{i=1}^n x_i^{2m} \\
\end{bmatrix}
\begin{bmatrix}
a_0 \\
a_1 \\
\vdots \\
a_m \\
\end{bmatrix}
=
\begin{bmatrix}
\sum_{i=1}^n y_i \\
\sum_{i=1}^n x_iy_i \\
\vdots \\
\sum_{i=1}^n x_i^my_i \\
\end{bmatrix}
$$
可以用线性方程组的求解方法,求出 $a_0, a_1, \cdots, a_m$ 的值,从而得到最小二乘法拟合的多项式 $f(x)$。
下面是 C 语言实现最小二乘法多项式拟合的代码,假设要拟合的数据存储在数组 x 和 y 中,要拟合的多项式的次数为 m:
```c
#include <stdio.h>
#include <stdlib.h>
#define MAXN 100
int main() {
int n, m, i, j, k;
double x[MAXN], y[MAXN], a[MAXN], s[MAXN][MAXN];
printf("Enter the number of data points: ");
scanf("%d", &n);
printf("Enter the degree of polynomial: ");
scanf("%d", &m);
printf("Enter the data points:\n");
for (i = 0; i < n; i++) {
scanf("%lf%lf", &x[i], &y[i]);
}
// 构造矩阵
for (i = 0; i <= m; i++) {
for (j = 0; j <= m; j++) {
s[i][j] = 0;
for (k = 0; k < n; k++) {
s[i][j] += pow(x[k], i + j);
}
}
}
// 构造向量
for (i = 0; i <= m; i++) {
a[i] = 0;
for (j = 0; j < n; j++) {
a[i] += pow(x[j], i) * y[j];
}
}
// 解线性方程组
for (i = 0; i <= m; i++) {
for (j = i + 1; j <= m; j++) {
double t = s[j][i] / s[i][i];
for (k = i; k <= m; k++) {
s[j][k] -= t * s[i][k];
}
a[j] -= t * a[i];
}
}
for (i = m; i >= 0; i--) {
for (j = i + 1; j <= m; j++) {
a[i] -= s[i][j] * a[j];
}
a[i] /= s[i][i];
}
// 输出拟合多项式
printf("The fitted polynomial is:\n");
printf("f(x) = ");
for (i = m; i >= 0; i--) {
printf("%.2f", a[i]);
if (i > 0) {
printf("x");
if (i > 1) {
printf("^%d", i);
}
printf(" + ");
}
}
printf("\n");
return 0;
}
```