「最小二乗法」の編集履歴（バックアップ）一覧はこちら

「最小二乗法」(2007/01/23 (火) 11:58:29) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

-課題1
課題2、3は値を変化

 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
 #define N 6	//データ数
 #define M 4 //次数
 #define EPS 10e-5
  
 double a[M+1][M+2];
 
 int jordan(void);
 
 double expor(double x , double n){
 	int i;
	double bar = 1.0;
	for( i = 0 ; i < n ; i ++ )
		bar *= x;
	return bar;
 }
 
 void main(){
	double x[N] = {0.0,1.0,2.0,3.0,3.1,5.0};
	double y[N] = {0.0,1.1,2.5,4.0,4.1,5.0};
	int i , j , k;
	
	for( i = 0 ; i < M+1 ; i++ )
		for( j = 0 ; j < M+2 ; j++ )
			a[i][j] = 0.0;
	
	for( i = 0 ; i <= M ; i++ ){
		for( j = 0 ; j <= M ; j++ ){
			for( k = 0 ; k < N ; k++ )
				a[i][j] += expor( x[k] , i + j);
 //			printf("%10.2lf ",a[i][j]);
		}
 //		printf("\n");
	}
	
	for( j = 0 ; j <= M ; j++ )
		for(k = 0 ; k < N ; k++ )
			a[j][M+1] += y[k] * expor(x[k] , j);
		
	if(jordan() == 1)
		return;
		
	for( i = 0 ; i <= M ; i++ )
		printf("A%d=%7.3lf\n" , i , a[i][M+1]);
	return;
 }
 
 int jordan(){
	double   s, q , pivot , del;
	int      i, j, k, l;    
	for( i = 0 ; i < M+1 ; i++ ){
		pivot = a[i][i];
		if( fabs(pivot) < EPS ){
			printf("ピボット数が許容誤差以下\n");
			return 1;
		}
		for( j = i ; j < M+2 ; j++ )
			a[i][j] = a[i][j]/pivot;
		for( k = 0 ; k < M+1 ; k++ ){
			del = a[k][i];
			for( j = i ; j < M+2 ; j++ )
				if( k != i )
				a[k][j] -= del*a[i][j];
		}
	}
	return 0;
 }