Transcript ニュートン法A(収束過程がわかる

ニュートン法の解の計算
情報電子工学系学科
電気電子工学コース・情報通信システム工学コース
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);}
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()
｛ 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,"%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>>
z1.dat
:データファイル
<<executable>>
re2.exe
:可視化プログラム
<<file>>
z2.csv
:表計算用可視化ファイル
サンプル点数
グラフ作成者
<<file>>
<<executable>>
z1.dat
re2.exe
:データファイル
:可視化プログラム
サンプル点数
グラフ作成者
<<file>>
z2.txt
:gnuplot用可視化ファイル
フローチャート(生成部）
開始
(main関数)
aの入力
開始(save関数)
ループ
=0;i<count;i++
(fp=fopen(“z1.dat”,”
wb”))=NULL
ループ
Count++
Save(count,b)
終了
exit(1)
i=0;i<c;i++
Newa=af(a)/de(a)
b[count]=newa
ループ終了
Newa,count
の出力
h[i]=f(i)
fwrite(&h[i],sizeo
f(double),1,fp)
fwrite(&d[i],sizeo
f(double),1,fp)
fclose(fp)
Count==max
exit(1)
i=1;i<d[0]+1;i+
+
i=0;i<d[0]+1;i+
+
h[i]=f(-[0]+i)
fwrite(&h[i],sizeo
f(double),1,fp)
終了
フローチャート（可視化部①）
開始(main関数)
i=0;i<k+1;i++
fprintf(fpw,"%f %d\n",a[i],b[i]);
if(x==1)
fprintf(fpw,"%f %f\n",a[i],e[i]);
fprintf(fpw,"%s",v)
(fpr=fopen(“z1.dat”,”rb”
))==NULL
exit(1)
1
fread(&a[i],sizeof(dou
ble),1,fpr)
fread(&b[i],sizeof(i
nt),1,fpr)
j=i
fprintf(fpw,"%f %f\n",a[i],e[
i]);
exit(1);
i=0;i<2*l+1;i
++
g=c[i]
fprintf(fpw,"%f",i-l);
fprintf(fpw," %f\n",c[i]);
i=0;i<2*a[0]+1;i+
+
fprintf(fpw,"%f",i-l);
fprintf(fpw," %f\n",c[i])
ｌの入力
kの入力
i=0;i<k;i++
fprintf(fpw,"%f %d\n",a[i+1
],b[i]);
fprintf(fpw,"%s",v)
Fread(&c[i],sizeof(do
uble),1,fp)
break
(fpw=fopen(
data_file,"w"
))==NULL
i=0;i<2*a[0]+1;i+
+
Xの入力
X==1||x==2
data_file="z2.txt";
i=0;i<k+1;i+
+
Fread(&e[i],sizeof(do
uble),1,fpr)
x==2
(fpw=fopen
( "z2.csv","
w"))==NU
LL
exit(1)
i=0;i<2*a[0]+1;i++)
fprintf(fpw,"%f",i-l);
fprintf(fpw," , %f\n",c[i
])
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の
初期値が２０に設定されているためであると考えられる。
