liezhuyuan.h

曲线拟合及最小二乘法的相关程序

/**********************************************************************
 * 预处理定义
**********************************************************************/
#include <stdio.h>
#include <math.h>

/**********************************************************************
 * 函数声明部分
**********************************************************************/
int ChangeMaxToTop(double **a,int row,int line,int m,int n);
int LieZhuYuan(double **a,int m,int n);
double QiuRoot(int num,double **a,int m,int n);
double Array(double **a,int i,int j,int n);
double DoubleAbs(double x);

/**********************************************************************
 * 函 数 名:getRoot
 * 功能描述:求根
 * 输    入:a -- 待求根方程系数矩阵
 *			 m -- 二维数组行数
 *			 n -- 二维数组列数
 * 输    出:x -- 根数组
 * 返    回:无
***********************************************************************/
void getRoot(double x[],double **a,int m,int n)
{
	int i;
	double **p;
	p = a;
	LieZhuYuan(a,m,n);

	for(i = 0;i < m;i ++)
	{
		x[i] = QiuRoot(i + 1,a,m,n);
		printf("x%d = %f\n",i + 1,x[i]);
	}
	/*
	x[1] = QiuRoot(2,p,3,4);
	printf("x2 = %f\n",x[1]);

	x[0] = QiuRoot(1,p,3,4);
	printf("x1 = %f\n",x[0]);
	*/
}
/**********************************************************************
 * 函 数 名:Array
 * 功能描述:计算在一维数组中的二维数据
 * 输    入:a -- 一维数组 
 *			 i -- 二维数组行下标
 *           j -- 二维数组列下标
 *			 n -- 二维数组列数
 * 输    出:无
 * 返    回:a[i][j]的值
**********************************************************************/
double Array(double **a,int i,int j,int n)
{
	double temp = *((double*)a + n * i + j);
	return temp;
}

/**********************************************************************
 * 函 数 名:ChangeMaxToTop
 * 功能描述:求第row行第line列中的最大的行，并与第row行交换位置
 * 输    入:a -- 待交换数组
 *			 row -- 需要交换的数组行下标
 *           line -- 需要交换的数组列下标
 *			 m -- 二维数组行数
 *			 n -- 二维数组列数
 * 输    出:交换后的数组
 * 返    回:1 -- 成功
 *			 0 -- 失败
**********************************************************************/
int ChangeMaxToTop(double **a,int row,int line,int m,int n)
{
	int i,j;
	double temp;
	temp = Array(a,row,line,n);
	j = row;
	for(i = row + 1;i < m;i++) //找到绝对值比当前主元绝对值更大的数
		if(DoubleAbs(temp) < DoubleAbs(Array(a,i,line,n)))	
		{
			temp = Array(a,i,line,n);
			j = i;
		}
	if(j < row || j > m)
		return 0;
	else if(j == row)	//当前主元最大，不需要调换
		return 1;
	for(i = 0;i < n;i++)
	{
		temp = Array(a,j,i,n);
		*((double*)a + n * j + i) = Array(a,row,i,n);
		*((double*)a + n * row + i) = temp;
	}
	return 1;
}

/**********************************************************************
 * 函 数 名:LieZhuYuan
 * 功能描述:列主元法消元
 * 输    入:a -- 系数矩阵
 *			 m -- 矩阵行数
 *			 n -- 矩阵列数
 * 输    出:a -- 销元后系数矩阵
 * 返    回:1 -- 成功
 *			 0 -- 失败
**********************************************************************/
int LieZhuYuan(double **a,int m,int n)
{
	int i,j,k,flag;
	double temp;
	for(i = 0;i < m - 1;i++)
	{
		flag = ChangeMaxToTop(a,i,i,m,n);	//交换最大主元
		if (flag == 0)
			return 0;
		for(k = i + 1;k < m;k++)
		{
			temp = Array(a,k,i,n) / (0 - Array(a,i,i,n));
			for(j = 0;j < n;j++)	//置主元下面的系数为0
				*((double*)a + n * k + j) = Array(a,k,j,n) + Array(a,i,j,n) * temp;
		}
	}
	return 1;
}

/**********************************************************************
 * 函 数 名:QiuRoot
 * 功能描述:求根
 * 输    入:num -- 要求的根
 *			 a -- 系数矩阵
 *			 m -- 矩阵行数
 *			 n -- 矩阵列数
 * 输    出:无
 * 返    回:第num个根值
**********************************************************************/
double QiuRoot(int num,double **a,int m,int n)
{
	int j;
	double temp = 0;		//临时变量，用来求和
	if(num == m)
		return Array(a,m - 1,n - 1,n) / Array(a,m - 1,m - 1,n);		//第m个根
	else
	{
		for(j = num;j < m;j++)
			temp = Array(a,num - 1,j,n) * QiuRoot(j + 1,a,m,n) + temp;	//递归调用
		return (Array(a,num - 1,n - 1,n) - temp) / Array(a,num - 1,num - 1,n);
	}
}

/**********************************************************************
 * 函 数 名:DoubleAbs
 * 功能描述:求double型数据的绝对值
 * 输    入:x -- double型数据
 * 输    出:无
 * 返    回:x的绝对值
**********************************************************************/
double DoubleAbs(double x)
{
	if (x < 0)
		return 0 - x;
	else
		return x;
}