有志による某国立高専某学科某クラスのページ。夏休み中の連絡とかにも使ってください。課題に迷ったらココ!
カテゴリ
全体
初めによくお読みください
おしらせ
英訳
代数幾何
解析
応用物理
熱力学
流れ
情報処理
設計
電子情報
そのほか
未分類
以前の記事
フォロー中のブログ
メモ帳
最新のトラックバック
ライフログ
検索
その他のジャンル
ファン
記事ランキング
ブログジャンル
画像一覧
情報処理課題L
SREです。最小自乗法のプログラム載せます。
ちなみに見やすいように若干プログラムに訂正を加えてあります。

#include
#include

#define N 1
#define N1 N+1
#define M 5

double x[M]={1.0,2.0,3.0,4.0,5.0};
double f[M]={2.0,2.5,2.9,3.5,4.4};
double A[N1][N1],y[N1];

int leastsq(void);
int sweep(void);

int main(void)
{
int i,j,sw;
double xx,wk,wg;

sw=leastsq();

for(i=0;i printf("y[%d]= %f\n",i,y[i]);

if(sw)
{
printf("x=");
scanf("%lf",&xx);

wk=0.0;
for(i=0;i {
wg=1.0;
for(j=0;j {
wg=wg*xx;
}

wk=wk+y[i]*wg;
}
printf("y=c+b+x^2…+a*x^N=%f\n",wk);
}
return 0;
}

int sweep(void)
{
int i,j,k,pivot_row;
double p,q,big,temp;
for(k=0;k {
for(i=0;i<(N1-k);i++)
{
big=fabs(A[k][k]);
if(big {
big=fabs(A[k+i][k]);
pivot_row=k+i;
}
}
if(big==0.0)
{
printf("このプログラムにおいて入力された配列は適用範囲外と認識されました。\n");
return 0;
}
if(k!=pivot_row)
{
for(j=0;j {
temp=A[k][j];A[k][j]=A[pivot_row][j];A[pivot_row][j]=temp;
}
temp=y[k];y[k]=y[pivot_row];y[pivot_row]=temp;
}
p=A[k][k];
for(i=0;i {
A[k][i]=A[k][i]/p;
}
y[k]=y[k]/p;
for(i=0;i {
q=A[i][k];
if(i!=k)
{
for(j=0;j {
A[i][j]=A[i][j]-q*A[k][j];
}
y[i]=y[i]-q*y[k];
}
}
}
return 1;
}

int leastsq()
{
double wk,wg1,wg2;
int i,j,k,m;

for(i=0;i {
for(j=i;j {
wk=0.0;
for(k=0;k {
wg1=1.0;
for(m=0;m {
wg1=wg1*x[k];
}

wg2=1.0;
for(m=0;m {
wg2=wg2*x[k];
}

wk=wk+wg1*wg2;
}

A[i][j]=wk;
A[j][i]=wk;
}
}
for(i=0;i {
wk=0.0;
for(k=0;k {
wg1=1.0;
for(m=0;m {
wg1=wg1*x[k];
}

wk=wk+wg1*f[k];
}
y[i]=wk;
}

sweep();

return 1;
}
[PR]
by m0511xx | 2009-02-05 10:07