nihe.c

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

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

/**********************************************************************
 * 函数声明部分
**********************************************************************/
void multiply(double **a,double **b,double **c,int row,int colum);
void bMultiply(double **a,double **b,double **c,int row,int colum);
int maxTwo(int x,int y);
int minTwo(int x,int y);
void transpose(double **a,double **b,int row,int colum);
void junct(double **a,double **b,double **c,int row,int colum);


/**********************************************************************
 * 函 数 名:multiply
 * 功能描述:求两个矩阵相乘
 * 输    入:a -- 左矩阵
 *			 b -- 右矩阵
 *			 row -- 矩阵行数
 *			 colum -- 矩阵列数
 * 输    出:c -- 结果矩阵
 * 返    回:无
***********************************************************************/
void multiply(double **a,double **b,double **c,int row,int colum)
{
	int i,j,k;
	for(i = 0;i < colum;i++)
	{
		for(k = 0;k < colum;k++)
		{
			double sum = 0;
			for(j = 0;j < row;j++)
			{
				sum = Array(b,j,i,colum) * Array(a,k,j,row) + sum;
			}
			*((double*)c + minTwo(row,colum) * i + k) = sum;
		}
	}
}

/**********************************************************************
 * 函 数 名:bMultiply
 * 功能描述:求根
 * 输    入:a -- 左矩阵
 *			 b -- 右矩阵
 *			 row -- 矩阵行数
 *			 colum -- 矩阵列数
 * 输    出:c -- 结果矩阵
 * 返    回:无
***********************************************************************/
void bMultiply(double **a,double **b,double **c,int row,int colum)
{
	int i,j,k;
	for(i = 0;i < row;i++)
	{
		double sum = 0;
		for(j = 0;j < colum;j++)
		{
			sum = Array(a,i,j,colum) * Array(b,j,0,1) + sum;
		}
		*((double*)c + 1 * i + 0) = sum;
	}
}

/**********************************************************************
 * 函 数 名:maxTwo
 * 功能描述:求两个数中较大的一个
 * 输    入:x -- 待比较数1
 *			 y -- 待比较数2
 * 输    出:无
 * 返    回:两个数中较大着
***********************************************************************/
int maxTwo(int x,int y)
{
	return x > y ? x : y;
}

/**********************************************************************
 * 函 数 名:minTwo
 * 功能描述:求两个数中较小的一个
 * 输    入:x -- 待比较数1
 *			 y -- 待比较数2
 * 输    出:无
 * 返    回:两个数中较小着
***********************************************************************/
int minTwo(int x,int y)
{
	return x > y ? y : x;
}

/**********************************************************************
 * 函 数 名:transpose
 * 功能描述:转置矩阵
 * 输    入:a -- 左矩阵
 *			 row -- 矩阵行数
 *			 colum -- 矩阵列数
 * 输    出:b -- a矩阵的转置
 * 返    回:无
***********************************************************************/
void transpose(double **a,double **b,int row,int colum)
{
	int i,j;
	for(i = 0;i < row;i++)
		for(j = 0;j < colum;j++)
			*((double*)b + row * j + i) = Array(a,i,j,colum);
}

/**********************************************************************
 * 函 数 名:junct
 * 功能描述:对接两个矩阵
 * 输    入:a -- 左矩阵
 *			 b -- 右矩阵
 *			 row -- 矩阵行数
 *			 colum -- 矩阵列数
 * 输    出:c -- 结果矩阵
 * 返    回:无
***********************************************************************/
void junct(double **a,double **b,double **c,int row,int colum)
{
	int i,j;
	for(i = 0;i < row;i++)
	{
		for(j = 0;j < colum;j++)
		{
			*((double*)c + (row + 1) * i + j) = Array(a,i,j,colum);
		}
		*((double*)c + (row + 1) * i + j) = Array(b,i,0,1);
	}
}

/**********************************************************************
 * 函 数 名:main
 * 功能描述:主函数
 * 输    入:无
 * 输    出:无
 * 返    回:无
***********************************************************************/
int main(void)
{
	
	double a[11][3] = {0,0,1,0.1*0.1,0.1,1,0.2*0.2,0.2,1,
		0.3*0.3,0.3,1,0.4*0.4,0.4,1,0.5*0.5,0.5,1,
		0.6*0.6,0.6,1,0.7*0.7,0.7,1,0.8*0.8,0.8,1,
		0.9*0.9,0.9,1,1,1,1};
	double b[11][1] = {-0.447,1.978,3.28,6.16,7.08,
		7.34,7.66,9.56,9.48,9.30,11.2};
	double a_zh[3][11],xishu_result[3][3],b_result[3][1],x_result[3][4];
	double x[3];
	transpose(a,a_zh,11,3);
	multiply(a_zh,a,xishu_result,11,3);
	bMultiply(a_zh,b,b_result,3,11);
	junct(xishu_result,b_result,x_result,3,3);
	getRoot(x,x_result,3,4);
	printf("a = %f\n",x[0]);
	printf("b = %f\n",x[1]);
	printf("c = %f\n",x[2]);
	
	/*
	double a[16][2] = {1,log(0.15),1,log(0.4),1,log(0.6),1,log(1.01),
		1,log(1.5),1,log(2.2),1,log(2.4),1,log(2.7),
		1,log(2.9),1,log(3.5),1,log(3.8),1,log(4.4),
		1,log(4.6),1,log(5.1),1,log(6.6),1,log(7.6)};
	double b[16][1] = {log(4.4964),log(5.1284),log(5.6931),log(6.2884),
		log(7.0989),log(7.5507),log(7.5106),log(8.0756),log(7.8708),
		log(8.2403),log(8.5303),log(8.7394),log(8.9981),log(9.1450),
		log(9.5070),log(9.9115)};
	double a_zh[2][16],xishu_result[2][2],b_result[2][1],x_result[2][3];
	double A_B[2];
	double a_root,b_root;
	transpose(a,a_zh,16,2);							//求aT
	multiply(a_zh,a,xishu_result,16,2);				//aT*a
	bMultiply(a_zh,b,b_result,2,16);				//aT*b
	junct(xishu_result,b_result,x_result,2,2);		//求系数矩阵
	getRoot(A_B,x_result,2,3);						//求根
	a_root = exp(A_B[0]);							//根a
	b_root = A_B[1];								//根b
	printf("a = %f\n",a_root);						//打印a
	printf("b = %f\n",b_root);						//打印b
	*/
}