1. 程式人生 > >最小二乘法C語言的實現

最小二乘法C語言的實現

1.實驗目的:

進一步熟悉曲線擬合的最小二乘法。

掌握程式語言字元處理程式的設計和除錯技術。

2.實驗要求:

輸入:已知點的數目以及各點座標 。

輸出:根據最小二乘法原理以及各點座標求出擬合曲線 。

3.程式流程:

(1)輸入已知點的個數;

(2)分別輸入已知點的X座標;

(3)分別輸入已知點的Y座標;

(4)通過呼叫函式,求出擬合曲線。

最小二乘法原理如下:

根據一組給定的實驗資料,求出自變數x與因變數y的函


4.難點,要根據已知座標點求出多項式的係數。

5.程式流程圖

6.原始碼(本程式已通過測試)

#include<stdio.h>
#include<math.h>
#define X 50
#define Y 50
float x[X],y[Y];
int n;//輸入的資料總組數即座標的總個數
void init();//初始化並輸入相關資料
void confrim();//確認輸入的資料
void deal();//根據輸入的座標點計算出擬合曲線
void modify();//用於修改輸錯的相應座標這樣可以避免一些資料重新輸入
void main()
{
	
	int select;
	system("color f1");//dos命令使介面變顏色
	init();//
	confrim();
	printf("請選擇要擬合成幾次多項式(提示:如果是一次函式就輸入1二次函式就輸入2):");
	scanf("%d",&select);//輸入你要選擇擬合的函式的次數
	deal(select);
}

void init()//初始化並輸入相關資料
{
	int i;
	printf ("\n*********************************************************\n");
	printf ("\n歡迎使用最小二乘法資料處理程式\n");
	printf ("\n請輸入您要處理的資料的組數(提示:程式定義一對x,y值為一組資料):");
	
	while(1)
	{
		scanf("%d",&n);
		if(n<=1)
		{
			printf("\n理的資料的組數不能小於或等於1");
			printf ("\n請重新輸入您要處理的資料的組數:");
		}
		
		else if(n>50)
		{
			printf ("\n對不起,本程式暫時無法處理50組以上的資料");
			printf ("\n請重新輸入您要處理的資料的組數:");
		} 
		else break;
	}
	
	for (i=0;i<n;i++)//輸入相應座標點將其存到數組裡
	{
		printf ("\n請輸入第%d個x的值x%d=",i+1,i+1);
		scanf ("%f",&x[i]);
		printf ("\n請輸入對應的y的值:y%d=",i+1);
		scanf ("%f",&y[i]);
	}
	system("color f2");//
	system("cls");//清屏
}

void deal(int select)//採用克萊默法則求解方程
{
	int i;
	float a0,a1,a2,temp,temp0,temp1,temp2;
	float sy=0,sx=0,sxx=0,syy=0,sxy=0,sxxy=0,sxxx=0,sxxxx=0;//定義相關變數
	for(i=0;i<n;i++)
	{
		sx+=x[i];//計算xi的和
		sy+=y[i];//計算yi的和
		sxx+=x[i]*x[i];//計算xi的平方的和
		sxxx+=pow(x[i],3);//計算xi的立方的和
		sxxxx+=pow(x[i],4);//計算xi的4次方的和
		sxy+=x[i]*y[i];//計算xi乘yi的的和
		sxxy+=x[i]*x[i]*y[i];//計算xi平方乘yi的和
	}
	temp=n*sxx-sx*sx;//方程的係數行列式
	temp0=sy*sxx-sx*sxy;
	temp1=n*sxy-sy*sx;
	a0=temp0/temp;
	a1=temp1/temp;
	if(select==1)
	{
		printf("經最小二乘法擬合得到的一元線性方程為:\n"); 
		printf("f(x)=%3.3fx+%3.3f\n",a1,a0); 
				system("pause");
	}
	temp=n*(sxx*sxxxx-sxxx*sxxx)-sx*(sx*sxxxx-sxx*sxxx)//方程的係數行列式
		+sxx*(sx*sxxx-sxx*sxx);
	temp0=sy*(sxx*sxxxx-sxxx*sxxx)-sxy*(sx*sxxxx-sxx*sxxx)
		+sxxy*(sx*sxxx-sxx*sxx);
	temp1=n*(sxy*sxxxx-sxxy*sxxx)-sx*(sy*sxxxx-sxx*sxxy)
		+sxx*(sy*sxxx-sxy*sxx);
	temp2=n*(sxx*sxxy-sxy*sxxx)-sx*(sx*sxxy-sy*sxxx)
		+sxx*(sx*sxy-sy*sxx);
	a0=temp0/temp;
	a1=temp1/temp;
	a2=temp2/temp;
	if(select==2)
	{
		printf("經最小二乘法擬合得到的二次近似方程為:\n"); 
		printf("f(x)=%3.3fx2+%3.3fx+%3.3f\n",a2,a1,a0); 
		system("pause");
	}
	
}
void modify()//修改輸錯的相應座標

{  
	int z;
	char flag;
	while(1)
	{
		printf("請輸入你要修改的是第幾組資料:");
		scanf("%d",&z);
		printf ("\n請輸入你要修改的第%d個x的值x%d=",z,z);
		scanf ("%f",&x[z-1]);
		printf ("\n請輸入你要修改的對應的y的值:y%d=",z);
		scanf ("%f",&y[z-1]);
		printf("是否繼續修改資料是Y否N:");
		getchar();
		scanf("%c",&flag);
		if(flag=='N'||flag=='n')
			break;
	}
	system("cls");//清屏
	confrim();
}
void confrim()
{
	char flag;
	int i;
	while(1)
	{
		for(i=0;i<n;i++)
		{
			printf ("請輸入第%d個x的值x%d=",i+1,i+1);
			printf ("%f",x[i]);
			printf ("  輸入對應的y的值:y%d=",i+1);
			
			printf("%f",y[i]);
			printf("\n");
		}
		printf("確認你輸入的資料是Y否N(即重新輸入)修改M:");
		getchar();
		scanf("%c",&flag);
		if(flag=='y'||flag=='Y')
			break;
		else if(flag=='n'||flag=='N')
			init();
		else 
		{
			modify();
			break;
		}
	}	
}

測試

1擬合二次曲線結果(書上的)

2.擬合一次曲線結果