ガウスのの単純消去法のプログラムです。
ガウスのの単純消去法のプログラムです。
前進消去の第k段階が終わった段階でaij,biが表示されるようにしたいんですがどうすればいいでしょうか↓
よろしくお願いします。
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define ERROR -1
#define EPS 1.0e-15
int gausssimp(double *, double *, double *, int);
int main(void)
{
double *a,*b,*x,val;
char s[32];
int i,j,n,result;
printf("Input size n=");
gets(s);
sscanf(s,"%d",&n);
if((a=(double *)calloc(n*n,sizeof(double)) ) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}
if((b=(double *)calloc(n,sizeof(double))) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}
if((x=(double *)calloc(n,sizeof(double))) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}
/* A,b の成分を入力 */
printf("----- A -----\n");
for(i=0; i<n; i++){
for(j=0; j<n; j++){
printf("a(%2d,%2d)=",i+1,j+1);
gets(s);
sscanf(s,"%lf",&val);
a[n*i+j]=val;
}
}
printf("\n----- b -----\n");
for(i=0; i<n; i++){
printf("b(%2d)=",i+1);
gets(s);
sscanf(s,"%lf",&val);
b[i]=val;
}
/* Gaussの単純消去法による求解 */
result = gausssimp(x,a,b,n);
/* 解の表示 */
if(result == ERROR){
printf("ERROR occurs. pivot 0\n");
} else {
printf("\n----- solution -----\n");
for(i=0; i<n; i++){
printf("x(%2d)=%.8e\n",i+1, x[i]);
}
}
free(a);
free(b);
free(x);
return 0;
}
int gausssimp(double *x, double *a, double *b, int n)
{
int i,j,k;
double tmp,p,sum;
/* step 1: 前進消去 */
/**** 追加 ****/
/* 前進消去の各段階を終えるごとに，式がどのように変化しているかわかるように表示する */
/**************/
for(k=0; k<n-1;k++){
if(a[n*k+k] == 0.0) { /* ピボットの値が0．割り算でエラーが起きる．*/
return ERROR;
} else {
/* k+1番目以降の式から x[k] の項を消去 */
for(i=k+1; i<n; i++){
p=a[n*i+k]/a[n*k+k];
for(j=0; j<n; j++){
a[n*i+j]=a[n*i+j]-p*a[n*k+j];
}
b[i]=b[i]-p*b[k];
printf("a[%d %d]=%d b[%d]=%d",i,j,a[n*i+j]-p*a[n*k+j],i,b[i]-p*a[n*k+j]);
↑これではできませんでした。。。
}
/**********************************/
}
printf("k=%d\n",k);
}
/* step 2: 後退代入 */
for(k=n-1; k>=0; k--){
if(fabs(a[n*k+k]) < EPS) return(ERROR);
sum=0.0;
for(j=k+1; j<n; j++) sum+=a[n*k+j]*x[j];
x[k]=(b[k] - sum)/a[k*n+k];
}
return 0;
}
補足
回答ありがとうございます。 ご指摘のとおりでした。 for(i=N;i>=0;i--) と変更したところエラーも出ず、正確に演算できました。 これから関数化してみたいと思います。