ニュートン法の解の計算 情報電子工学系学科 電気電子工学コース・情報通信システム工学コース 10024061 小林 祐樹 10024074 佐藤 景治 10024075 佐藤 建甫 10024077 佐藤 宣
はじめに 今回の課題内容 プログラムの製作①(ソースコード<生成部>) プログラムの製作②(ソースコード<可視部>) 計算結果 アクティビティ図(データフロー) フローチャート グラフへの可視化 今回の課題の考察
課題内容 出力 数値&表示用ファイルを表示 出力に変換 ニュートン法による数値計算 プログラム入力 方程式を入力する
プログラムの製作①ソースコード<生成部> #include <stdio.h> #include <stdlib.h> #include <math.h> #define max 1000 #define eps 1.0e-5 #define F x*x-10*x+25 #define DF 2*x-10 double f(double x); double df(double x); void save(int c,double d[max]); int main() { int count,i,j=0; double a,newa,b[max]; FILE *fp; count=0; printf("ニュートン法を用いてf(x)=0の解を求めます。\n"); printf("初期値を入力してください。\n"); scanf("%lf",&a); b[0]=a; for(;;) { count++; newa=a-f(a)/df(a); b[count]=newa; if(fabs(newa-a)<eps) break; a=newa; if(count==max) { printf("収束しませんでした\n"); exit(1);} }
ソースコード②<生成部> printf("解の値は%f\n収束するのに%d回かかりました。\n", newa,count); printf("計算過程を以下に示します。\n"); for(i=0;i<count+1;i++) { printf("x[%d] %f\n",i,b[i]);} save(count,b); return 0;} double f(double x) {return F;} double df(double x) { return DF;} void save(int c, double d[max]) { int i,j=0,k=0; double h[max],l[max]; char a[3]="\n" ; FILE *fp;
ソースコード③<生成部> if( (fp=fopen( "z1.dat","wb"))==NULL) { fprintf(stderr,"ファイルを作成できません。\n"); exit(1);} for(i=0;i<d[0]+1;i++) {fwrite(&d[i],sizeof(double),1,fp); fwrite( &j,sizeof(int),1,fp);} { h[i]=f(-d[0]+i); fwrite(&h[i],sizeof(double),1,fp);} for(i=1;i<d[0]+1;i++) {h[i]=f(i); fwrite(&h[i],sizeof(double),1,fp); } for(i=0;i<c+1;i++) { l[i]=f(d[i]); fwrite(&l[i],sizeof(double),1,fp);} fclose(fp);}
出力結果 ニュートン法を用いてf(x)=x*x-10*x+25=0の解を求めます。 初期値を入力してください。 20 解の値は5.000007 収束するのに21回かかりました。 計算過程を以下に示します。 x[0] 20.000000 x[11] 5.007324 x[1] 12.500000 x[12] 5.003662 x[2] 8.750000 x[13] 5.001831 x[3] 6.875000 x[14] 5.000916 x[4] 5.937500 x[15] 5.000458 x[5] 5.468750 x[16] 5.000229 x[6] 5.234375 x[17] 5.000114 x[7] 5.117188 x[18] 5.000057 x[8] 5.058594 x[19] 5.000029 x[9] 5.029297 x[20] 5.000014 x[10] 5.014648 x[21] 5.000007
プログラム製作<可視部> #include <stdio.h> #include <stdlib.h> #include <math.h> #define max 1000 #define eps 1.0e-5 #define F "x*x-10*x+25" void main() { int count,x,b[max],i,j,k; double newa,a[max],c[max],e[max],f[max],g,l; char v[3]="\n", *data_file; FILE *fpr,*fpw,*gp; count=0; if( (fpr=fopen( "z1.dat","rb"))==NULL) { fprintf(stderr,"ファイルを開けません1。\n"); exit(1)}
プログラム製作<可視部> printf("どちらのファイルを作成しますか?当てはまる数字を入力してください。\n"); printf("txt形式ファイル--->1\n"); printf("csv形式ファイル--->2\n"); while(1){ scanf("%d",&x); if(x==1 || x==2) break; printf("入力し直してください。"); } printf("初期値を入力してください。"); scanf("%lf",&l); printf("何回で収束しましたか?"); scanf("%d",&k); for(i=0;i<l+1;i++) { fread(&a[i],sizeof(double),1,fpr); fread(&b[i],sizeof(int),1,fpr); j=i;} for(i=0;i<2*l+1;i++) {fread( &c[i],sizeof(double),1, fpr); g=c[i];} for(i=0;i<k+1;i++) { fread(&e[i],sizeof(double),1,fpr);}
プログラム製作<可視部> if(x==1) { /*データファイルの作成*/ data_file="z2.txt"; if( (fpw=fopen(data_file,"w"))==NULL) { fprintf(stderr,"ファイルを開けません2。\n"); exit(1);} for(i=0;i<2*a[0]+1;i++) {fprintf(fpw,"%f",i-l); fprintf(fpw," %f\n",c[i]);} fprintf(fpw,"%s",v); for(i=0;i<k;i++) { fprintf(fpw,"%f %d\n",a[i],b[i]); fprintf(fpw,"%f %f\n",a[i],e[i]); fprintf(fpw,"%s",v);} { fprintf(fpw,"%f %f\n",a[i],e[i]); fprintf(fpw,"%f %d\n",a[i+1],b[i]); fprintf(fpw,"%s",v);}}
プログラム製作<可視部> if(x==2) { if( (fpw=fopen( "z2.csv","w"))==NULL) { fprintf(stderr,"ファイルを開けません2。\n"); exit(1);} for(i=0;i<2*a[0]+1;i++) {fprintf(fpw,"%f",i-l); fprintf(fpw," , %f\n",c[i]);} fprintf(fpw,"%s",v); for(i=0;i<k;i++) { fprintf(fpw,"%f , %d\n",a[i],b[i]); fprintf(fpw,"%f , %f\n",a[i],e[i]); fprintf(fpw,"%s",v);}} fclose(fpr); fclose(fpw); if(x==1) {gp = popen("gnuplot -persist","w"); fprintf(gp,"set grid\n"); fprintf(gp, "set xrange [-%f:%f]\n",l,l); fprintf(gp, "set yrange [-%f:%f]\n",g,g); fprintf(gp, "plot \"%s\" with lines linetype 1 title \"x*x-2*x+1\"\n",data_file); fprintf(gp,"replot %s \n",F); pclose(gp);}}
アクティビティ図(データフロー) <<file>> <<executable>> z1.dat z2.csv re2.exe :データファイル :可視化プログラム :表計算用可視化ファイル サンプル点数 グラフ作成者 <<file>> <<executable>> <<file>> z2.txt re2.exe z1.dat :gnuplot用可視化ファイル :データファイル :可視化プログラム サンプル点数 グラフ作成者
フローチャート(生成部) h[i]=f(i) exit(1) i=0;i<c;i++ fclose(fp) 終了 開始 (main関数) 開始(save関数) ループ =0;i<count;i++ aの入力 (fp=fopen(“z1.dat”,”wb”))=NULL i=1;i<d[0]+1;i++ ループ Save(count,b) h[i]=f(i) exit(1) 終了 i=0;i<c;i++ Count++ fwrite(&h[i],sizeof(double),1,fp) Newa=a-f(a)/de(a) fwrite(&d[i],sizeof(double),1,fp) b[count]=newa fclose(fp) Count==max exit(1) i=0;i<d[0]+1;i++ 終了 ループ終了 h[i]=f(-[0]+i) Newa,countの出力 fwrite(&h[i],sizeof(double),1,fp)
フローチャート(可視化部①) i=0;i<k+1;i++ i=0;i<k;i++ exit(1) exit(1); j=i 1 fprintf(fpw,"%f %d\n",a[i],b[i]); fprintf(fpw,"%f %f\n",a[i],e[i]); fprintf(fpw,"%s",v) i=0;i<k+1;i++ if(x==1) 開始(main関数) data_file="z2.txt"; fread(&a[i],sizeof(double),1,fpr) (fpr=fopen(“z1.dat”,”rb”))==NULL i=0;i<k;i++ (fpw=fopen(data_file,"w"))==NULL fread(&b[i],sizeof(int),1,fpr) exit(1) fprintf(fpw,"%f %f\n",a[i],e[i]); fprintf(fpw,"%f %d\n",a[i+1],b[i]); fprintf(fpw,"%s",v) exit(1); j=i 1 i=0;i<2*a[0]+1;i++ x==2 Xの入力 fprintf(fpw,"%f",i-l); fprintf(fpw," %f\n",c[i]); i=0;i<2*l+1;i++ (fpw=fopen( "z2.csv","w"))==NULL X==1||x==2 Fread(&c[i],sizeof(double),1,fp) exit(1) i=0;i<2*a[0]+1;i++ break g=c[i] i=0;i<2*a[0]+1;i++) fprintf(fpw,"%f",i-l); fprintf(fpw," %f\n",c[i]) fprintf(fpw,"%f",i-l); fprintf(fpw," , %f\n",c[i]) lの入力 i=0;i<k+1;i++ kの入力 Fread(&e[i],sizeof(double),1,fpr) i=0;i<k;i++
フローチャート(可視化部②) i=0;i<k;i++ fclose(fpr); fclose(fpw); 終了 fprintf(fpw,"%f , %d\n",a[i],b[i]); fprintf(fpw,"%f , %f\n",a[i],e[i]); fprintf(fpw,"%s",v) fclose(fpr); fclose(fpw); 終了
グラフへの可視化 図 1 図 2 図 1;gnuplotによるグラフ可視化 図 2; excel によるグラフ可視化
課題に対する考察 生成部でforを駆使していくことにより、出力結果では2 0回の計算結果を算出しているが、今回製作したプログ ラムでは、どんな回数でも表示できるようになっているの ではないかと考えている。 グラフへの可視化において、図 2のexcelを用いて表 示したグラフでは、x>0の部分で途切れているが、ここ では今回自ら設定した方程式f(x)=x^2-10*x+25=0の 初期値が20に設定されているためであると考えられる。