back.cpp

计算方法中最小二乘法程序设计的源代码,具有很强的通用性

/*
#include "stdio.h"
#include "math.h"

#define M 11	//采样个数
#define N 4		//未知数个数

double	Y[M]={	3.6599, 2.6261, 2.0555, 1.6241, 1.3675, 1.1514, 0.9599, 0.8654,
0.7766, 0.7277, 0.6331};
double	X[M][N+1]={	{1, 25, 8.4634, 6.3548, 1.5561},
{1, 35, 5.8744, 4.6514, 1.2122},
{1, 45, 4.6508, 3.7096, 1.0355}, 
{1, 55, 3.4163, 2.8078, 0.9106}, 
{1, 65, 2.8506, 2.2576, 0.8349}, 
{1, 75, 2.2308, 1.9725, 0.7325}, 
{1, 85, 1.7921, 1.6598, 0.6401}, 
{1, 95, 1.5319, 1.4081, 0.6160}, 
{1, 105,1.3758, 1.1834, 0.5707}, 
{1, 115,1.1659, 1.1148, 0.5270},
{1, 125,0.9435, 0.9263, 0.4800}};
double	B[N+1];
double	YJ[M];
double	A[N+1][N+2];

void main()
{	
	//计算增广矩阵元素
	int	i, j, k;
	
	for(k=0; k<=N; k++)
	{
		for(j=0; j<=N; j++)
		{
			A[k][j] = 0;
			
			for(i=0; i<M; i++)
				A[k][j] = A[k][j]+X[i][k]*X[i][j];
		}
		
		A[k][N+1] = 0;
		
		for(i=0; i<M; i++)
			A[k][N+1] = A[k][N+1]+Y[i]*X[i][k];
	}
	
	//列主元高斯－约当消去法解方程组
	double	c, ftemp;
	int		p;
	
	for(k=0;k<=N;k++)
	{
		c=A[k][k];
		p=k;
		
		for(i=k+1; i<=N; i++)
		{
			if(fabs(A[i][k]) > c)
			{
				c=A[i][k];
				p=i;
			}
		}
		if (fabs(c) < 1.0e-9) 
			return;
		
		for(j=k; j<=N+1; j++){
			ftemp=A[k][j];
			A[k][j]=A[p][j];
			A[p][j]=ftemp;
		}
		
		for(i=0; i<=N; i++){
			if (i==k) continue;
			
			for(j=k+1; j<=N+1; j++)
				A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k];
		}
	}
	
	for(i=0; i<=N; i++)
		B[i]=A[i][N+1]/A[i][i];
	
	
	for(i=0; i<M; i++)
	{
		YJ[i]=0;
		for(j=0; j<=N; j++)
			YJ[i]+=B[j]*X[i][j];
	}
	
	p=0;
	return;
}
*/