C/C++语言曲线的拟合的最小二乘方法

应用计算方法C语言程序：03

接 应用计算方法C语言程序：02

编写C语言曲线的拟合的最小二乘方法，以计算方法课本例题8-3为例。

1、获取矩阵A和b

#define POW_N 3 //表示非线性函数阶次
double A[POW_N][POW_N] = { 0 };
double b[POW_N] = { 0 };
double data1[12] = { 256,201,159,61,77,40,17,25,103,156,222,345 };
void coefficient_matrix(double *data, int data_len)//求解系数矩阵
{
	for (int i = 0; i < POW_N; i++)
	{
		for (int j = 0; j < POW_N; j++)
		{
			for (int k = 1; k <= data_len; k++)
			{
				A[i][j] += pow(k, i + j);
			}
		}
	}

	for (int i = 0; i < POW_N; i++)
	{
		for (int k = 1; k <= data_len; k++)
		{
			b[i] += pow(k, i) * (data[k-1]);
		}
	}
}

2、使用LU分解求解拟合曲线系数

double L[3][3] = { 0 }, U[3][3] = { 0 };
#define POW_N 3
double A[POW_N][POW_N] = { 0 };
double b[POW_N] = { 0 };
double data1[12] = { 256,201,159,61,77,40,17,25,103,156,222,345 };
void coefficient_matrix(double *data, int data_len)；//求解系数矩阵
void Doolittle(double a[3][3]);LU分解
int main()
{
	//获取系数矩阵
	coefficient_matrix(data1,12);

	Doolittle(A);

	double b1[3][3] = { 0 };
	double b2[3][3] = { 0 };
	Matrix_inverse(U, 3, b1);//求矩阵U的逆
	Matrix_inverse(L, 3, b2);//求矩阵L的逆
	
	
	//x=U^(-1) * L^(-1) * b
	double b3[3][3] = { 0 };
	//U^(-1) * L^(-1)
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				b3[i][j] += b1[i][k] * b2[k][j];
			}
		}
	}
	//x=U^(-1) * L^(-1) * b
	double b4[3] = { 0 };
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			b4[i]+= b3[i][j] * b[j];
		}
	}

	//输出方程的解
	cout << endl;
	for (int i = 0; i < 3; i++)
	{
		cout << "x" << i + 1 << " = " << b4[i] << " ";
	}

	cout << endl << endl << endl;

	return 0;
}

3、运行结果