おしえて

シンプソン公式の数値計算

締切済

質問者：アヒルのももちゃん
質問日時：2017/03/26 21:41
回答数：2件

以下のプログラムを実行しました。
補間点のファイルは以下のとおりです。
0
0.1963495408
0.3926990817
0.5890486225
0.7853981634
0.9817477042
1.1780972451
1.3744467859
1.5707963268

プログラム＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
/***
simpson公式の検証プログラム
教科書　：　ｐ46
***/
#include<stdio.h>
#include<math.h>

#define DIVISION 100
#define INPUT 100
#define FILE_NAME 1024

/***値の交換を行うプログラム***/
void swap(double *a, double *b)
{
double temp = *a;
*a = *b;
*b = temp;
}

/***バブルソートを行う関数***/
void bsort(double a[], int n)
{
int i, j;

for(i=0; i<n-1; i++) {
for(j=n-1; j>i; j--) {
if(a[j-1]>a[j]) {
swap(&a[j-1], &a[j]);
}
}
}
}

/***多項式の値を返す関数 ***/
double fi_x(double x, double xi[], int i)
{
double f = (x - xi[i+1])*(x - xi[i+2])/(xi[i] - xi[i+1])/(xi[i] - xi[i+2])*sin(xi[i])
+ (x - xi[i]) *(x - xi[i+2])/(xi[i+1] - xi[i]) /(xi[i+1] - xi[i+2])*sin(xi[i+1])
+ (x - xi[i]) *(x - xi[i+1])/(xi[i+2] - xi[i]) /(xi[i+2] - xi[i+1])*sin(xi[i+2]);
return f;
}

/***積分の値を返す関数 :積分区間[x[i], x[i+2]]***/
double sekibun(int i, double xi[])
{
double temp;

temp = 1.0/6*(
(xi[i+2] - xi[i])*(2*xi[i] + xi[i+2] - 3*xi[i+1])*sin(xi[i])
+ pow((xi[i] - xi[i+2]),3)*sin(xi[i+1])/(xi[i+1] - xi[i])/(xi[i+1] - xi[i+2])
+ (xi[i+2] - xi[i])*(2*xi[i+2] + xi[i] - 3*xi[i+1])*sin(xi[i+2])/(xi[i+2] - xi[i+1])
);
return temp;
}

/***メインプログラム***/
int main(void)
{
/*変数とファイル名の定義*/
FILE *fp, *gp;
double xi[INPUT+1]; //補間点の座標
double X[DIVISION+1], Y[DIVISION+1]; //近似関数の座標
double I = 0.; //積分結果を格納する変数
int N = 0, cnt = 0, i, j;
char frame[FILE_NAME];

/*補間点の読み取り*/
printf("ファイル名：");
scanf("%s", frame);
if( (fp = fopen(frame, "r")) == NULL )
printf("ファイルをオープンできません\n");
else {
while (fscanf(fp, "%lf", &xi[cnt++]) == 1) {
}
fclose(fp);
}

N = cnt - 2;

/*補間点は偶数であるという条件*/
if(!N%2)
N -= 1;

/*読み取った値をバブルソート*/
bsort(xi, N);

/*補間関数のx座標*/
X[0] = xi[0];
for(i=1; i<=DIVISION; i++)
X[i] = X[i-1] + (xi[N] - xi[0]) / DIVISION;

/*補間関数のy座標*/
for(i=0; i<=N-2; i+2) {
for(j=0; j<=DIVISION; j++) {
if(xi[i] <= X[j] && X[j] <= xi[i+2])
Y[j] = fi_x(X[j], xi, i);
}
}

/*積分の実行*/
for(i=0; i<=N-2; i+2)
I += sekibun(i,xi);

printf("[%lf,%lf]での積分の結果は%lfです。\n", X[0], X[DIVISION], I);

/*グラフの作成*/
gp = popen("gnuplot4 -persist", "w");
fprintf(gp, "set multiplot\n");
fprintf(gp, "set xrange[%lf : %lf]\n", X[0], X[DIVISION]);
fprintf(gp, "set yrange[%lf : %lf]\n", -1.5, 1.5);
fprintf(gp, "plot sin(x) with points\n");
fprintf(gp, "plot '-' with lines linetype 1 title \"sin\" \n"); //データの送信

for(i=0; i<DIVISION; i++)
{
fprintf(gp, "%lf\t%lf\n", X[i], Y[i]);
}

fprintf(gp, "e\n");

pclose(gp);

return 0;
}

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行結果
[]$gcc simpson_rule-t1.c -lm
[]$./a.out
ファイル名：d1.dat
^Z
[6]+ 停止 ./a.out
プログラムが終了しません。