ニュートン法の解の計算

ニュートン法の解の計算 情報電子工学系学科電気電子工学コース・情報通信システム工学コース 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(intc,double d[max]); int main() { intcount,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]) { inti,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);} for(i=0;i<d[0]+1;i++) { 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() ｛intcount,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,"%s",v); for(i=0;i<k;i++) { 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>> <<file>> z1.dat z2.csv re2.exe :データファイル :可視化プログラム :表計算用可視化ファイルサンプル点数グラフ作成者 <<file>> <<executable>> <<file>> z2.txt re2.exe z1.dat :gnuplot用可視化ファイル :データファイル :可視化プログラムサンプル点数グラフ作成者

フローチャート(生成部） 　　　開始　　　　(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)

フローチャート（可視化部①） 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]) ｌの入力 i=0;i<k+1;i++ kの入力 Fread(&e[i],sizeof(double),1,fpr) i=0;i<k;i++

フローチャート（可視化部②） 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); 終了

グラフへの可視化 図　１図　２図　１；gnuplotによるグラフ可視化図　２； excel によるグラフ可視化

課題に対する考察 • 生成部でforを駆使していくことにより、出力結果では２０回の計算結果を算出しているが、今回製作したプログラムでは、どんな回数でも表示できるようになっているのではないかと考えている。 • グラフへの可視化において、図　２のexcelを用いて表示したグラフでは、x>0の部分で途切れているが、ここでは今回自ら設定した方程式f(x)=x＾2-10*x+25=0の初期値が２０に設定されているためであると考えられる。

ニュートン法の解の計算

ニュートン法の解の計算

Presentation Transcript