0
点赞
收藏
分享

微信扫一扫

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


#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;
}


举报

相关推荐

0 条评论