（C言語）プログラムが動きません。

質問

下のC言語のプログラムを書いたのですが、動きません。上手く作動していれば画面に「結果をrunge1.csvに書き込みました」と表示されてテキストファイルが作成されるはずですが、プログラムを実行しても何のメッセージも表示されずファイルも作成されません。ですが、bcpadのエラーメッセージの欄には何も表示されず、どこを直せばいいのか分かりません。アドバイスを下さい。お願いします。
/*ルンゲ・クッタ法による微分方程式y'+y*y=xの解*/
#include<stdio.h>
#include<math.h>
#define N 100
int main(void)
{
char *file="runge1.csv";
double x0=0.0,x1=5.0,y0=0.0,x,y,h,k1,k2,k3,k4;
int i;
FILE *fp;

h=(x0-x1)/N;
x=x0; y=y0;
fp=fopen(file,"w");
fprintf(fp,"x,Runge-y\n");
fprintf(fp,"%3.2f,%3.5f\n",x,y);

for(i=1;i<=N;i++){
x+=h;
k1=(x-y*y)*h;
k2=(x-y*y+k1/2.0)*h;
k3=(x-y*y+k2/2.0)*h;
k4=(x-y*y+k3)*h;
y=(x-y*y)+(k1+2*k2+2*k3+k4)/6.0;
fprintf(fp,"%3.2f,%3.5f\n",x,y);
}
fclose(fp);
printf("結果を%sに書き込みました。\n",file);

return(0);
}

回答No.5

和泉 博（@hiroshi09s） · ベストアンサー

#4さんに追加補足すると、４次のルンゲ・クッタ法を最初からファイルに落とさず、結果を見るようにします。それでOKなら、ファイルに落とされることをお勧めします。↓のプログラムでは動作を確認後、 fp を戻せば出来上がりです。
　なお、UNIX系でしたら、「 ./a.out>runge1.csv 」のようにリダイレクトして終わりです。

#include<stdio.h>
#include<math.h>
#define N 100

double f(double x, double y) {
return x - y * y;
}

int main(void) {
char *file="runge1.csv";
double x0=0.0,x1=5.0,y0=0.0,x,y,h,k1,k2,k3,k4;
int i;
FILE *fp;

h=(x1-x0)/N;
x=x0; y=y0;
fp = stdout; //fp=fopen(file,"w");
fprintf(fp,"x,Runge-y\n");
fprintf(fp,"%3.2f,%3.5f\n",x,y);

for(i=1;i<=N;i++){
k1=f(x,y)*h;
k2=f(x+0.5*h, y+0.5*k1)*h;
k3=f(x+0.5*h, y+0.5*k2)*h;
k4=f(x+h, y+k4)*h;
y += (k1+2.0*(k2+k3)+k4)/6.0;
x += h;
fprintf(fp,"%3.2f,%3.5f\n",x,y);
}
//fclose(fp);
printf("結果を%sに書き込みました。\n",file);

return(0);
}

参考URL：

http://www.me.tokushima-u.ac.jp/~miw/wadai/IP3/06/10/10.html

お礼

アドバイスありがとうございます。私はWindowsを使用していますので、ご指南いただいた点を直したら無事動きました。オイラー法のときは,
x+=h;
y+=(x-y*y)*h;
で上手くいったのですが・・・。まだプログラミング初心者なので難しかったです。

回答No.6

和泉 博（@hiroshi09s）

#5です。プログラム転記ミスです。失礼しました。 m(_ _)m

k4=f(x+h, y+k4)*h;

　↓↓↓ 訂正 ↓↓↓

k4=f(x+h, y+k3)*h;

お礼

ありがとうございます。理解の助けになりました。

回答No.4

kmee

x0→x1へ向かって近似していくので
>h=(x0-x1)/N;
ではなく
h=(x1-x0)/N
ですね。

また、現在のx,yを使って次のx,yを近似するのですから
>x+=h;
と先に足してはいけません。

y'=f(x,y)=x-y*y について
k2=f(x+h/2 , y + k1/2) * h
ですから、
>k2=(x-y*y+k1/2.0)*h;
という式にはなりません。k3,k4も同様

> y=(x-y*y)+(k1+2*k2+2*k3+k4)/6.0;
= f(x,y) +(k1+2*k2+2*k3+k4)/6.0;
ですから、これも計算式が違います。

コツとして、
double f(double x, double y) {
return x - y * y ;
}
と、関数を定義すると、計算部分は公式のまま書けます。

お礼

アドバイスをありがとうございます。結構うっかりミスがあったのですね。何度もチェックしたつもりだったのですが・・・。指南いただいた点を修正したら無事動きました。

回答No.3

noname#144013

こんにちは。

「ルンゲ・クッタ法」云々などの数学的部分については解りませんので、それ以外
の部分についてコメントさせて頂きます。

bcpadをご使用されているようですので、コンパイラは、「Borland C++ Compiler」
をご使用でしょうか？（OSは、Windows系？）

当方で、ご提示のプログラムを、BCPad ＋ Borland C++ Compiler 5.5 の環境で
コンパイル＆実行してみました。

その結果、実行時に例外エラーが発生しているようでした。

そこで、プログラムに例外処理を追加して、エラー要因を確認してみました。
注）この例外処理は、OSがWindows系＋処理系が対応している場合、を前提と
　　しています。

以下は、その例外処理を追加したソースです。

■例外処理追加ソース
注１）ファイルオープン時のエラーチェック追加、CSVファイルのフォーマット変更
　　　（ループ回数の値出力の追加）などの変更も行っています。
注２）インデント等のため全角スペースを入れています。
　　　ご使用の際は半角スペースorタブに置換して下さい。

/////↓ここから///////////////
/*
* ルンゲ・クッタ法による微分方程式[ y'+y*y=x ]の解
*/
#include <windows.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <process.h>

#define N 100

int Except(LPEXCEPTION_POINTERS pExp);

int main(void)
{
　　char *file="runge1.csv";
　　double x0=0.0;
　　double x1=5.0;
　　double y0=0.0;
　　double x, y, h, k1, k2, k3, k4;
　　int i;
　　FILE *fp;

　　h=(x0-x1)/N;
　　x=x0; y=y0;

　　fp=fopen(file,"w");
　　if(fp==NULL){
　　　　printf("file open error.\n");
　　　　exit(1);
　　}

　　i=0;
　　fprintf(fp, "i,x,Runge-y\n");
　　fprintf(fp, "%3d,%3.2f,%3.5f\n", i, x, y);

　　for(i=1;i<=N;i++){
　　　　__try {
　　　　　　x+=h;
　　　　　　k1=(x-y*y)*h;
　　　　　　k2=(x-y*y+k1/2.0)*h;
　　　　　　k3=(x-y*y+k2/2.0)*h;
　　　　　　k4=(x-y*y+k3)*h;
　　　　　　y=(x-y*y)+(k1+2*k2+2*k3+k4)/6.0;
　　　　　　fprintf(fp,"%3d,%3.2f,%3.5f\n", i, x, y);
　　　　}
　　　　__except( Except(GetExceptionInformation()),
　　　　　　　　　　EXCEPTION_EXECUTE_HANDLER ){
　　　　　　//例外処理
　　　　　　printf("\n[Variable Infomation]:\n");
　　　　　　printf(" i = %d\n", i);
　　　　　　printf(" h = %f\n", h);
　　　　　　printf(" x = %f\n", x);
　　　　　　printf(" k1= %f\n", k1);
　　　　　　printf(" k2= %f\n", k2);
　　　　　　printf(" k3= %f\n", k3);
　　　　　　printf(" k4= %f\n", k4);
　　　　　　printf(" y = %f\n", y);
　　　　　　fclose(fp);
　　　　　　exit(2);
　　　　}
　　}

　　fclose(fp);
　　printf("結果を%sに書き込みました。\n",file);

　　return(0);
}

/*
* Except : 例外情報の取得＆表示
*/
int Except(LPEXCEPTION_POINTERS pExp)
{
　　char strExpCode[128];
　　if(pExp){
　　　　strcpy(strExpCode, "UNKNOW");
　　　　switch(pExp->ExceptionRecord->ExceptionCode){
　　　　case EXCEPTION_ACCESS_VIOLATION:
　　　　　　strcpy(strExpCode, "EXCEPTION_ACCESS_VIOLATION");
　　　　　　break;
　　　　case EXCEPTION_ARRAY_BOUNDS_EXCEEDED:
　　　　　　strcpy(strExpCode, "EXCEPTION_ARRAY_BOUNDS_EXCEEDED");
　　　　　　break;
　　　　case EXCEPTION_BREAKPOINT:
　　　　　　strcpy(strExpCode, "EXCEPTION_BREAKPOINT");
　　　　　　break;
　　　　case EXCEPTION_DATATYPE_MISALIGNMENT:
　　　　　　strcpy(strExpCode, "EXCEPTION_DATATYPE_MISALIGNMENT");
　　　　　　break;
　　　　case EXCEPTION_FLT_DENORMAL_OPERAND:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_DENORMAL_OPERAND");
　　　　　　break;
　　　　case EXCEPTION_FLT_DIVIDE_BY_ZERO:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_DIVIDE_BY_ZERO");
　　　　　　break;
　　　　case EXCEPTION_FLT_INEXACT_RESULT:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_INEXACT_RESULT");
　　　　　　break;
　　　　case EXCEPTION_FLT_INVALID_OPERATION:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_INVALID_OPERATION");
　　　　　　break;
　　　　case EXCEPTION_FLT_OVERFLOW:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_OVERFLOW");
　　　　　　break;
　　　　case EXCEPTION_FLT_STACK_CHECK:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_STACK_CHECK");
　　　　　　break;
　　　　case EXCEPTION_FLT_UNDERFLOW:
　　　　　　strcpy(strExpCode, "EXCEPTION_FLT_UNDERFLOW");
　　　　　　break;
　　　　case EXCEPTION_ILLEGAL_INSTRUCTION:
　　　　　　strcpy(strExpCode, "EXCEPTION_ILLEGAL_INSTRUCTION");
　　　　　　break;
　　　　case EXCEPTION_IN_PAGE_ERROR:
　　　　　　strcpy(strExpCode, "EXCEPTION_IN_PAGE_ERROR");
　　　　　　break;
　　　　case EXCEPTION_INT_DIVIDE_BY_ZERO:
　　　　　　strcpy(strExpCode, "EXCEPTION_INT_DIVIDE_BY_ZERO");
　　　　　　break;
　　　　case EXCEPTION_INT_OVERFLOW:
　　　　　　strcpy(strExpCode, "EXCEPTION_INT_OVERFLOW");
　　　　　　break;
　　　　case EXCEPTION_INVALID_DISPOSITION:
　　　　　　strcpy(strExpCode, "EXCEPTION_INVALID_DISPOSITION");
　　　　　　break;
　　　　case EXCEPTION_NONCONTINUABLE_EXCEPTION:
　　　　　　strcpy(strExpCode, "EXCEPTION_NONCONTINUABLE_EXCEPTION");
　　　　　　break;
　　　　case EXCEPTION_PRIV_INSTRUCTION:
　　　　　　strcpy(strExpCode, "EXCEPTION_PRIV_INSTRUCTION");
　　　　　　break;
　　　　case EXCEPTION_SINGLE_STEP :
　　　　　　strcpy(strExpCode, "EXCEPTION_SINGLE_STEP");
　　　　　　break;
　　　　case EXCEPTION_STACK_OVERFLOW:
　　　　　　strcpy(strExpCode, "EXCEPTION_STACK_OVERFLOW");
　　　　　　break;
　　　　default:
　　　　　　break;
　　　　}
　　　　printf("[Exception Infomation]:\n");
　　　　printf(" Code\t= %d[0x%08X] : %s\r\n"
　　　　　　" Flags\t= 0x%08X\r\n"
　　　　　　" Address= 0x%08X\r\n"
　　　　　　,pExp->ExceptionRecord->ExceptionCode
　　　　　　,pExp->ExceptionRecord->ExceptionCode
　　　　　　,strExpCode
　　　　　　,pExp->ExceptionRecord->ExceptionFlags
　　　　　　,pExp->ExceptionRecord->ExceptionAddress);
　　}
　　return 0;
}
/////↑ここまで///////////////

＜上記プログラムの実行結果＆出力CSVファイル＞
◎実行結果
[Exception Infomation]:
Code　 = -1073741679[0xC0000091] : EXCEPTION_FLT_OVERFLOW
Flags　 = 0x00000000
Address= 0x00401231

[Variable Infomation]:
i = 22
h = -0.050000
x = -1.100000
k1= 2.368268460679573460000000000000000000000e+216
k2= 2.309061749162584184000000000000000000000e+216
k3= 2.310541916950509031000000000000000000000e+216
k4= 2.252741364832048024000000000000000000000e+216
y = -4.505533302063516528000000000000000000000e+217

◎出力CSVファイル
i,x,Runge-y
　0,0.00,0.00000
　1,-0.05,-0.04756
　2,-0.10,-0.09727
　3,-0.15,-0.15169
　4,-0.20,-0.21213
　5,-0.25,-0.28061
　6,-0.30,-0.36027
　7,-0.35,-0.45640
　8,-0.40,-0.57863
　9,-0.45,-0.74654
10,-0.50,-1.00575
11,-0.55,-1.48538
12,-0.60,-2.66949
13,-0.65,-7.39690
14,-0.70,-52.71162
15,-0.75,-2643.71838
16,-0.80,-6648378.04575
17,-0.85,-42045225928516.06250
18,-0.90,-1681584354667049036000000000.00000
19,-0.95,-2.689816127625114148000000000000000000000e+54
20,-1.00,-6.882250301579525437000000000000000000000e+108
21,-1.05,-4.505533302063516528000000000000000000000e+217

上記の結果から判るように、エラーの要因は、forループ内の22回目のループの
演算での、浮動小数点演算のオーバーフロー発生によるもののようです。
※質問者さんの現象も、同様な要因によるものかもしれません。

以上の内容から、対策としては、

　１．オーバーフローにならないように、各パラメータの設定、ループ回数の設定
　　　などを行う。

　２．演算の前に、オーバーフローになりそうかどうか、値のチェックを行う。

などの処理が必要だと思われます。
以上、見当違いの場合はすみません。

■参考サイト
※例外処理については、下記サイトを参考にさせて頂きました。
デバッグ関連の処理
http://hp.vector.co.jp/authors/VA014436/prg_memo/windows/vctips/046.html

VC++構造化例外メモ
http://www.ne.jp/asahi/hishidama/home/tech/vcpp/seh.html

※例外コード（ExceptionCode）の意味などについては、下記サイトをご覧下さい。
GetExceptionCode 関数
http://msdn.microsoft.com/ja-jp/library/cc428942.aspx

お礼

詳しいアドバイスをありがとうございます。BCPad ＋ Borland C++ Compiler 5.5を私も使用しています。実行結果を見ていて思ったのですが、hが負数になるはずがないのでそこで間違えていたようです。他にも細かい点を修正したら無事動きました。

回答No.2

Trick--o--

デバッガを使わない（使えない）なら、
一行ずつprintfを入れてどこまで実行できているか調べる。

お礼

アドバイスありがとうございます。無事原因が分かりました。「デバッガ」というのを初めて知りました。便利そうですね。今度使ってみようかと思います。

回答No.1

shinh

取り敢えず、　fopen　の　戻り値で　エラーが発生していないか　チェックを入れたほうが良いですね。

お礼

アドバイスをありがとうございました。もう一度チェックしたらどうもforループ内に問題があったようです。そこを修正したら無事動きました。