二次関数で近似

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
でるた

二次関数で近似

#1

投稿記事 by でるた » 6年前

昨日も似たような質問をしましたが、今度は二次関数で近似するプログラムです
計算結果の数値を正常に直したいのです
問題があれば指摘お願いします

ちなみに、int PIVOTとint GAUSSは、予め課題と一緒に用意されたプログラムです

コード:

#include <stdio.h>
#include <math.h>
#define EPS 1.0e-5
#define N 2 /*配列の行数と列数の-1*/
#define n 5 /*入力するデータの数*/

int main(void){
	int i,k;
	double x[n],y[n];
	double sum_x=0,sum_x2=0,sum_x3=0,sum_x4=0,sum_y=0,sum_xy=0,sum_x2y=0;
	double ma[N][N],mb[N],mx[N];

	for(i=0;i<n;i++){
		printf("x座標を入力[%d]:",i);
		scanf("%lf",&x[i]);
		printf("y座標を入力[%d]:",i);
		scanf("%lf",&y[i]);
	}

	for(k=0;k<n;k++){
		sum_x+=x[n];
			ma[0][1]=ma[1][0]=sum_x;
		sum_x2+=pow(x[n],2);
			ma[0][2]=ma[2][0]=sum_x2;
		sum_x3+=pow(x[n],3);
			ma[1][2]=ma[2][1]=sum_x3;
		sum_x4+=pow(x[n],4);
			ma[2][2]=sum_x4;
		sum_y=y[n];
			mb[0]=sum_y;
		sum_xy+=x[n]*y[n];
			mb[1]=sum_xy;
		sum_x2y+=pow(x[n],2)*y[n];
			mb[2]=sum_x2y;
	}

	printf("f(x)=%g+%gx+%gx^2",mx[0],mx[1],mx[2]);

	return 0;
}
	
int PIVOT(int *num,int nr,int k,double ma[])
/*nr:行数  ma:配列ma  mb:配列mb  mx:配列mx*/
{
	int i;
	double aa,bb;
	*num=k; 
	aa=fabs(ma[nr*k+k]);
		for(i=k+1;i<nr;i++)
		{
			if(fabs(ma[nr*i+k])>aa)
				{
					*num=i;
					aa=fabs(ma[nr*i+k]);
				}
		}
		if(fabs(aa)<=EPS) return (1);
		if(*num==k)       return (0);
		for(i=k;i<nr;i++)
			{
				bb=ma[nr*k+i];
				ma[nr*k+i]=ma[nr*(*num)+i];
				ma[nr*(*num)+i]=bb;
			}
			return 0;
}

int GAUSS(int nr,double ma[],double mb[],double mx[])
{
	int i,j,k;
	int num; /*交換行番号*/
	double cc; /*置換変数*/
	for(k=0;k<nr-1;k++)
		{
			if(PIVOT(&num,nr,k,ma)!=0) return(1);
			if(num!=k)
				{
					cc=mb[num];mb[num]=mb[k];mb[k]=cc;
				}
			for(i=k+1;i<nr;i++)
				{
					cc=ma[nr*i+k]/ma[nr*k+k];
					for(j=k+1;j<nr;i++)
						ma[nr*i+j]=ma[nr*i+j]-cc*ma[nr*k+j];
					mb[i]=mb[i]-cc*mb[k];
				}
		}
	for(k=nr-1;k>=0;k--)
		{
			if(fabs(ma[nr*k+k])<=EPS) return(1);
			cc=0.0;
			for(j=k+1;j<nr;j++)
				cc+=ma[nr*k+j]*mx[j];
			mx[k]=(mb[k]-cc)/ma[nr*k+k];
		}
	return 0;
}
ーーーコンパイル結果ーーー
x={1.0,2.0,3.0,4.0,5.0}
y={7.987,2.986,1.986,2.224,5.678}
を入力した場合

f(x)=14.86+8.3157x+1.2963x^2 ・・・と表示されるはずですが

f(x)=6.37431e-314+9.88131e-323+1^2 ・・・と表示してしまいます

アバター
usao
記事: 1566
登録日時: 6年前

Re: 二次関数で近似

#2

投稿記事 by usao » 6年前

(1)そもそも計算処理自体がまともに走っているのかどうかが,読み手には不明です.
 コードを貼ってあとは読め,という状態で投げるのではなく,「どういう方法で近似するのか」くらいは説明するべきかと思います.
(2)関数PIVOTとGAUSSはmain()から全く呼ばれておらず,
 コード内に存在している意図が不明です.用いないのが正しいのでしょうか?
(3)Nの値が2なのに対して,配列アクセスが[2]でされており,領域外参照になっています.

でるた

Re: 二次関数で近似

#3

投稿記事 by でるた » 6年前

説明不足でした、ごめんなさい。

sum_xはΣ[i-1,n]xi、sum_x2はΣ[i-1,n]x^2i
sum_x3はΣ[i-1,n]x^3i、sum_4はΣ[i-1,n]x^4i、の計算をして
配列maに格納します

sum_yはΣ[i-1,n]yi、sum_xyはΣ[i-1,n]xiyi
sum_x2yはΣ[i-1,n]x^2iyiの計算をして
配列mbに格納します

mbとmaの逆行列を掛けることで、配列mxが求めることができ、
mx[0]をa0、mx[1]をa1、mx[2]をa2とすると
f(x)=a0+a1x+a2x^2 ...の式が求められます


不足だったらまた教えて下さい

アバター
usao
記事: 1566
登録日時: 6年前

Re: 二次関数で近似

#4

投稿記事 by usao » 6年前

とりあえずやりたいことは データ点群に対する最小二乗当てはめ ですね.

コード:

//※面倒なのでΣを省略している
[ 1    x    x^2 ][ a0 ]   [ y ]
[ x    x^2  x^3 ][ a1 ] = [ xy ]
[ x^2  x^3  x^4 ][ a2 ]   [ x^2*y ]
から,左辺の行列の逆行列を右辺にかけて,[a0 a1 a2]^T を得ようというわけですね.

…で,現状,この式の各項を計算する箇所しかなく,逆行列云々の話がありませんので,当然,まともな結果は出ないわけですが.
(何の初期化もされていないmx[]の内容を表示しているだけ)
まずは,その実装を追加するか,
あるいは現状使われていない関数がその処理なのだとすれば,main()から適切に呼び出して使うなりする必要がありますね.

その他,
・ma[0][0]に値が入れられていないように見受けます.
・ma, mb は,forが終わってから代入すればよいのでは?
というのも気になりますね.

でるた

Re: 二次関数で近似

#5

投稿記事 by でるた » 6年前

指摘していただいた
・main()にint PIVOT,GAUSSを入れる
・ma[0][0]にnを代入
・ma,mbへの値の代入をforが終わってから
をしてみました。

数値が少し変わりましたが、やはり望んでいる値ではありませんでした。

mx[]が初期化されていないとなると
GAUSS内のmx[k]=(mb[k]-cc)/ma[nr*k+k]が働いていない
ということでしょうか。

アバター
usao
記事: 1566
登録日時: 6年前

Re: 二次関数で近似

#6

投稿記事 by usao » 6年前

>(3)Nの値が2なのに対して,配列アクセスが[2]でされており,領域外参照になっています.
これについては解決していますか?

>mx[]が初期化されていない
これは あなたが最初に提示したコードで,11行目で未初期化状態で定義したmx[]が,その後一切さわられることもなく
37行目でprintf()での表示対象になっている ということです.
あなたの現在のコード状態は不明ですので,現状どうなっているかわかりません.

ma,mbにちゃんと値を設定できたなら,あとは
GAUSS()に正しく引数を私さえすれば答えがでるんじゃないでしょうか?
(関数名等の雰囲気的に,おそらくガウスの消去法だと思う=「maの逆行列を求めてmbに掛ける」という操作ではないと思います.)

>ちなみに、int PIVOTとint GAUSSは、予め課題と一緒に用意されたプログラムです
とのことなので,
・これら関数の実装はバグっていないはず
・使い方ややっていること等については説明をうけているはず
と思うので.

閉鎖

“C言語何でも質問掲示板” へ戻る