人気ブログランキング | 話題のタグを見る


by m0511xx
有志による某国立高専某学科某クラスのページ。夏休み中の連絡とかにも使ってください。課題に迷ったらココ!
情報処理課題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;
}
by m0511xx | 2009-02-05 10:07