最小二乗法 - パラボラの説明書

• 課題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;
}