手助けお願いします。
問題は x=(-4,-3,-2,-1,0,1,2,3,4,5,6)
y=(13.4,7.2,2.7,1.2,0.6,3.1,6.9,12.5,20.8,31.7,44.0)
を観測点とし、それに最もよく近似すると思われるf(x)を求めることです。
ガウス・ヨルダン法を使って解きます。
ガウス・ヨルダンのプログラムは確実にあっています。
わからないのはそれに渡す配列の作り方。
main関数。
結果をエクセルのグラフで表すこと。
作ったエラーだらけのソース貼ります。
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define ERROR -1 //エラーリターン値
#define NORMAL 0 //ノーマルリターン値
#define mabs(x) ( ( (x) > 0 ) ? (x):-(x)) //絶対値を返すマクロ
#define DBL_EPSILON 0.000000001
//Gauss_Jに渡す配列を作る↓
double setmatrixa(double x,int m,int n){
int i,j,k;
double a[100];
for(i=0;i<m;i++){
for(j=0;j<m;j++){
for(k=0;k<n;k++){
a[(m+2)*i+j]+=Mypow(x[k],(i+j));
}
}
}
return(a);
}
//Gauss_Jに渡す配列を作る↓
double setmatrixb(double y,int m,int n){
int i,k;
double b[100];
for(i=0;i<m;i++){
for(k=0;k<n;k++){
b+=y[k]*Mypow(x[k],i);
}
}
return(b);
}
//POW関数だとなぜかエラーが出るため、代用の関数↓
double Mypow(double x,int s){
double buf_x;
int i;
buf_x=x;
if(s==0){
x=1;
}
else{
for(i=0;i<s-1;i++){
x=x*buf_x;
}
}
return(x);
}
//ここから
int Gauss_J(int n,double *a,double *b);
int partial(int n,int p,double *a,double *b);
void mswap(double *a,double *b);
int Gauss_J(int n,double *a,double *b){
int p,i,j,l; //カウンタ
double pivot,c; //ピボット値
for(p=0;p<n;p++){ //1行目からn行目まで繰り返す
if((partial(n,p,(double *)a,(double *)b))==ERROR){
return(ERROR);
}
pivot=a[p*n+p]; //ピボットを取得
for(i=p;i<n;i++){ //p行目のp列からn列まで
a[p*n+i]/=pivot; //係数行列のp行をpivotで割る
}
b[p]/=pivot; //定数ベクトルのp行をpivotで割る
for(l=0;l<n;l++){ //1行目からn行目まで
if(l!=p){ //p行を除いて
c=a[l*n+p]; //掃きだす
for(j=p;j<n;j++){
a[l*n+j]-=c*a[p*n+j];
}
b[[/url]-=c*b[p];
}
}
}
return(NORMAL);
}
int partial(int n,int p,double *a,double *b){
int l,m,max; //l:カウンタ max:絶対値が最大値の係数を持つ行
double pivot; //ピボット
for(l=max=p,pivot=a[p*n+p];l<n;l++){ //n行目まででp列の絶対値が最大のものを選ぶ
if(mabs(pivot)<mabs(a[l*n+p])){
max=l;
pivot=a[l*n+p];
}
}
if(mabs(pivot)<DBL_EPSILON){ //選ばれたピボットが0ならERROR(-1)を返す
return(ERROR);
}
if(p!=max){ //選ばれた行がp行でなければ、行の交換を行う
for(m=0;m<n;m++){ //係数行列のpとmaxを交換
mswap(&a[p*n+m],&a[max*n+m]);
}
mswap(&b[p],&b[max]); //定数のp行とmax行を交換
}
return (NORMAL);
}
void mswap(double *a,double *b){
double temp;
temp=*a;
*a=*b;
*b=temp;
}
//ここまで教科書丸写し
void main(void){
double x[100],y[100],a[100],b[100];
int i,n,m;
int r;
printf("データ数 = ");
scanf("%d",&n);
printf("%d\n",n);
printf("次元数 = ");
scanf("%d",&m);
printf("%d\n",m);
for(i=0;i<n;i++){
printf("x = ");
scanf("%lf",&x);
printf("%f",x);
printf("y = ");
scanf("%lf",&y);
printf("%f",y);
}
a=setmatrixa(x,m,n);
b=setmatrixb(y,m,n);
r=Gauss_J(m,a,b);
if(r==0){ //エラー判定
printf("isw = %d\n",r);
for(i=0;i<m;i++){
printf("X%d = %lf\n",i+1,b); //結果出力
}
}
else{
printf("isw = %d",r);
}
}
メニュー
カレンダー
Clock
最新メンバー
投稿数
統計データ
Link to us
カウンター