利用LU分解法的多项式拟合实验

#include<iostream>
#include<vector>
#include<math.h>
using namespace std;


void print(int n, double**A,char *name) {
cout << "接下来打印矩阵" << name<<endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
cout << A[i][j] << "   ";
}
cout << endl;
}
cout << "打印完毕" << endl;
}


void print(int n, double*s, char *name) {
cout << "接下来打印向量" << name<<endl;
for (int i = 0; i < n; i++) {
cout << s[i] << "   ";
}
cout << "打印完毕" << endl;
}


void Init(vector<double>&x,vector<double> &y) {
cout << "有多少组数据要输入：" << endl;
int num;
cin >> num;
if (num == 0) {
cout << "输入了数字0" << endl;
exit(0);
}
printf("请按行并以空格分隔输入这些数据点（例：2.5 65.8）：\n");
for (int i = 0; i < num; i++) {
double a, b;
cin >> a >> b;
x.push_back(a);
y.push_back(b);
}
cout << "当前输入的数据为：" << endl;
for (int i = 0; i < num; i++) {
cout << "(" << x[i] << "," << y[i] << ")";
if ((i+1) % 5 == 0)cout << endl;
}
cout << endl;
}


void Jacobi(const int n, double *s, double *p) {
double **A;
A = new double*[n + 1];


//根据中间结果获取系数矩阵
for (int i = 0; i < n + 1; i++) {
A[i] = new double[n + 1];
for (int j = 0; j < n + 1; j++) {
A[i][j] = s[i + j];
}
}


print(n + 1, A,"A");
double error = 0.001;


//初始化
double **U, **L;
U = new double*[n + 1];
L = new double*[n + 1];
for (int i = 0; i < n + 1; i++) {
U[i] = new double[n + 1];
L[i] = new double[n + 1];
for (int j = 0; j < n + 1; j++) {
U[i][j] = 0.0;
L[i][j] = 0.0;
}
}


//计算LU矩阵
for (int k = 0; k < n + 1; k++) {
for (int j = k; j < n + 1; j++) {
double sum = 0;
for (int m = 0; m < k; m++)sum += L[k][m] * U[m][j];
U[k][j] = A[k][j] - sum;
}
if (abs(U[k][k]) < error) {
cout << "求解失败" << endl;
exit(0);
}
for (int i = k + 1; i < n + 1; i++) {
double sum = 0;
for (int m = 0; m < k; m++)sum += L[i][m] * U[m][k];
L[i][k] = (A[i][k] - sum) / U[k][k];
}
}


print(n + 1, U, "U");
print(n + 1, L, "L");


double *Y = new double[n + 1];
double *X = new double[n + 1];

//求解LY=b
Y[0] = p[0];
for (int i = 1; i < n + 1; i++) {
double sum = 0;
for (int j = 0; j <= i - 1; j++)sum += L[i][j] * Y[j];
Y[i] = p[i] - sum;
}


//求解UX=Y
X[n] = Y[n] / U[n][n];
for (int i = n-1; i >=0; i--) {
double sum = 0;
for (int j = i+1; j < n+1; j++)sum += U[i][j] * X[j];
X[i] = (Y[i] - sum)/U[i][i];
}


//输出参数
cout << "参数为：" << endl;
for (int i = 0; i < n + 1; i++) {
printf("%c=%f\n", 'a' + i, X[i]);
}
cout << endl;


//计算结果
double result = 0;
for (int i = 0; i < n + 1; i++) {
result += X[i] * pow(2010, i);
}
cout << "结果为：" << result << endl;
}


void Calculate(const int n, const vector<double>x, const vector<double> y) {

double *s = new double[2 * n + 1];
double *p = new double[n + 1];


//初始化
for (int i = 0; i < 2 * n + 1; i++)
s[i] = 0.0;
for (int i = 0; i < n + 1; i++)
p[i] = 0.0;


//保留中间结果
for (int i = 0; i < (int)x.size(); i++) {
for (int j = 0; j < n + 1; j++) {
p[j] += y[i] * pow(x[i], j);//p[i]=y*x的i次方
}
for (int j = 0; j < 2 * n + 1; j++) {
s[j] += pow(x[i], j);  //s[i]=X的i次方
}
}


print(n + 1, p,"P");


Jacobi(n,s,p);
}


int main() {
int n;
cout << "请输入需要拟合的多项式次数" << endl;
cin >> n;
vector<double>x,y;
Init(x, y);
Calculate(n, x, y);
system("pause");
}