はさみ撃ちによる方程式の解を求めるプログラミング
ただし、emacsで書いたので、若干違うところがあるので注意！
<stdio.h>
<math.h>
<unistd.h>
extern char *optarg;
extern int optind, opterr, optopt;
double e_pow;
double func( double x )
{
//return sin( x ) - 0.2;
return x*(x*(x-6)+3)+16;
}
void iterate( double a, double b, int n )
{
double ya, yb, c, yc, delta;
ya = func(a);
yb = func(b);
c = a + (a-b)*ya/(yb-ya);
yc = func(c);
n++;
printf( "%03d> %18.14lf %18.14lf %18.14lf %16.6le¥n", n, a, c, b, yc );
if (yc == 0) {
printf( "Ans(%03d): %20.14lf %16.6le¥n", n, c, yc );
} else if (fabs(c-a) < fabs(c)*e_pow||fabs(c-b) < fabs(c)*e_pow) {
printf( "Ans(%03d): %20.14lf %16.6le¥n", n, c, yc );
} else if ( (a == c) || (c == b) ) {
printf( "Ans(%03d): %20.14lf %16.6le¥n", n, c, yc );
} else if (ya * yc < 0) {
iterate(a, c, n);
} else {
iterate(c, b, n);
}
}
int main(int argc,char* argv[]){
double a, b, ya, yb;
double xfrom, xto;
int xsteps, xidx,e_gosa=15;
int option;
opterr = 0;
while( 1 ){
option = getopt(argc, argv, "he:d");
if( option == -1 ) break;
switch( option ){
case 'h' : printf("-%c.¥n-e:significant figure¥n", option ); break;
case 'e' : printf("-%c with arg %s.¥n", option, optarg);e_gosa=atoi(optarg); break;
case '?' : printf("Unknown option '%c'¥n", optopt); break;
}
optarg = NULL;
}
printf("x: from to steps>");
scanf( "%lf %lf %d", &xfrom, &xto, &xsteps);
if(1<=e_gosa&&e_gosa<=15){
e_pow=pow(0.1,e_gosa);
}else{
e_gosa=15;
}
a = xfrom;
if ((ya = func(a)) == 0) {
printf( "Ans(%03d): %20.14lf %16.6le¥n", 0, a, ya );
}
for( xidx=1; xidx<=xsteps; xidx++ ) {
b = (xfrom*(xsteps-xidx)+xto*xidx)/xsteps;
yb = func(b);
if (yb == 0) {
printf( "Ans(%03d): %20.14lf %16.6le¥n", 0, b, yb );
} else if (ya * yb < 0) {
iterate(a, b, 0);
}
a = b; ya = yb;
}
return 0;
}
結果
./a.out -e 10
x: from to steps>-10 10 20
001> -2.00000000000000 -1.21428571428571 -1.00000000000000 1.719752e+00
002> -2.00000000000000 -1.27125232381274 -1.21428571428571 4.352996e-01
003> -2.00000000000000 -1.28539180858941 -1.27125232381274 1.066664e-01
004> -2.00000000000000 -1.28883984775152 -1.28539180858941 2.592916e-02
005> -2.00000000000000 -1.28967703306446 -1.28883984775152 6.290747e-03
006> -2.00000000000000 -1.28988008692532 -1.28967703306446 1.525494e-03
007> -2.00000000000000 -1.28992932366956 -1.28988008692532 3.698867e-04
008> -2.00000000000000 -1.28994126191066 -1.28992932366956 8.968400e-05
009> -2.00000000000000 -1.28994415648563 -1.28994126191066 2.174494e-05
010> -2.00000000000000 -1.28994485830875 -1.28994415648563 5.272309e-06
011> -2.00000000000000 -1.28994502847371 -1.28994485830875 1.278331e-06
012> -2.00000000000000 -1.28994506973212 -1.28994502847371 3.099457e-07
013> -2.00000000000000 -1.28994507973569 -1.28994506973212 7.514983e-08
014> -2.00000000000000 -1.28994508216117 -1.28994507973569 1.822093e-08
015> -2.00000000000000 -1.28994508274925 -1.28994508216117 4.417869e-09
016> -2.00000000000000 -1.28994508289184 -1.28994508274925 1.071156e-09
017> -2.00000000000000 -1.28994508292641 -1.28994508289184 2.597140e-10
Ans(017): -1.28994508292641 2.597140e-10
001> 2.00000000000000 2.75000000000000 3.00000000000000 -3.281250e-01
002> 2.00000000000000 2.71111111111111 2.75000000000000 -4.040604e-02
003> 2.00000000000000 2.70635428172781 2.71111111111111 -4.762693e-03
004> 2.00000000000000 2.70579403496692 2.70635428172781 -5.583900e-04
005> 2.00000000000000 2.70572835635464 2.70579403496692 -6.542584e-05
006> 2.00000000000000 2.70572066096012 2.70572835635464 -7.665296e-06
007> 2.00000000000000 2.70571975936834 2.70572066096012 -8.980588e-07
008> 2.00000000000000 2.70571965373872 2.70571975936834 -1.052156e-07
009> 2.00000000000000 2.70571964136326 2.70571965373872 -1.232695e-08
010> 2.00000000000000 2.70571963991337 2.70571964136326 -1.444207e-09
011> 2.00000000000000 2.70571963974350 2.70571963991337 -1.692015e-10
Ans(011): 2.70571963974350 -1.692015e-10
001> 4.00000000000000 4.40000000000000 5.00000000000000 -1.776000e+00
002> 4.40000000000000 4.53703703703704 5.00000000000000 -5.035500e-01
003> 4.53703703703704 4.57288284549471 5.00000000000000 -1.241658e-01
004> 4.57288284549471 4.58154252927663 5.00000000000000 -2.954927e-02
005> 4.58154252927663 4.58359328153976 5.00000000000000 -6.972593e-03
006> 4.58359328153976 4.58407662563530 5.00000000000000 -1.641980e-03
007> 4.58407662563530 4.58419041748746 5.00000000000000 -3.864876e-04
008> 4.58419041748746 4.58421719997163 5.00000000000000 -9.096091e-05
009> 4.58421719997163 4.58422350320611 5.00000000000000 -2.140733e-05
010> 4.58422350320611 4.58422498664512 5.00000000000000 -5.038110e-06
011> 4.58422498664512 4.58422533576488 5.00000000000000 -1.185693e-06
012> 4.58422533576488 4.58422541792835 5.00000000000000 -2.790464e-07
013> 4.58422541792835 4.58422543726509 5.00000000000000 -6.567208e-08
014> 4.58422543726509 4.58422544181588 5.00000000000000 -1.545558e-08
015> 4.58422544181588 4.58422544288689 5.00000000000000 -3.637389e-09
016> 4.58422544288689 4.58422544313894 5.00000000000000 -8.560406e-10
Ans(016): 4.58422544313894 -8.560406e-10
考察
★自分の与えた入力データの正当性 (設定した関数に対して何故そのような入力データで良いのかを論ぜよ)
求める方程式を手計算で一階微分すると、極値が ２＋√３)、２ー√３となっていて、そのときの方程式の値が、(およそ３．７３２と０.２６８とすると)ー４.３９と１６.３９になる。よってこのあいだにひとつ解があり,その外に２つ解があることが予想され,おおよそx=10のとき方程式の値は正(y＝４４６),x＝ー１０のとき方程式の値は負(ー１６１４)となる。
求める方程式は三次方程式なので、以上の条件より、ー１０から１０の値に３つの解があると推測される。ステップ数は、極値の間を飛ばさない程度の値(xを１刻み)にしたいので、２０とした。
★自分の与えた収束条件の正当性と、解に期待できる精度
収束条件を以下の様にした。
fabs(c-a) < fabs(c)*e_pow||fabs(c-b) < fabs(c)*e_pow
e_powはここではe-10とした。
このとき、cのスケール(例えばC＝１００００００００００)かけるe_pow、に対してfabs(c-a)がおおきければ(例えばfabs(c-a)=１００００００００)、さらにはさみ撃つので、(Cのスケールのe-10分の１)程度の幅でXの精度を求められると期待でき、収束条件として適当ではないかと思う。
または(｜｜)で区切っているのは、凸のグラフと凹のグラフに対応するためである(aかbは固定される)。
★自分の得た解が実際に備えている精度
求めた解の値は、
Ans(017): -1.28994508292641 2.597140e-10
Ans(011): 2.70571963974350 -1.692015e-10
Ans(016): 4.58422544313894 -8.560406e-10
である。
上の考察により、およそ小数点以下９桁までが正しい精度だと考えられる。なぜなら、(Cのスケールのe-10分の１)程度の幅が予想されるので、小数点第１０桁目はこの幅の影響があると思うからである。
Ans(017): -1.289945082
Ans(011): 2.705719639
Ans(016): 4.584225443
までが実際に備えている精度である。
実際上の結果でC値の変化を見てみると、小数点第９桁目までは変化が終了している様に見える。
★解の持つ精度を入力から制御する方法 (可能なら各自のプログラムに組み込め)
組み込んだつもりです。
★二分法、ニュートン法と比較した収束の様子に関する考察
二分法は、一回の再帰呼出しで２進数表示で１桁づつ解の精度があがる。(二分法という手法だからこそ)
実際実行例を見るとそうなっているはず。
ニュートン法は、解にちかい初期値から始めた場合、一回の再帰呼出しで２次の収束になると言われている。
収束の様子をグラフにしようと思ったが、諸事情によりあきらめた。
最後にいちおう、ニュートン法の収束の様子を観察するためのプログラムとその結果を添付する。
★二分法、ニュートン法との使い分けに関する考察 (ガイドラインを提示せよ)
１．微分可能である。
２.導関数が求められる。
３.ある程度解の範囲がわかっている。(どこの解を示すかわからないため)
この３つが満たせるなら、収束が非常にはやいニュートン法を使うべきである。
これを１つでも満たせないなら、二分法を使うべきである。
★おまけ
NR法のメインのループに誤差のループを加えた。(変えた部分のみ)
for(ii=4;ii<15;ii++){
EPSILON=pow(10,-1*ii);
printf("ii=%d¥n",ii);
for( xidx=0; xidx<=xsteps; xidx++ ) {
x = (xfrom * (xsteps - xidx) + xto * xidx) / xsteps;
if ((y = func(x)) == 0) {
printf( "Ans(%03d): %20.14lf %16.5le %12.5le %12.5le¥n",
0, x, func(x*(1-EPSILON)), y, func(x*(1+EPSILON)) );
} else {
printf( "%03d> %20.14lf¥n", 0, x );
iterate( x, 0 );
}
}
}
これで、変化が見れるはず！！
x: from to steps> -5 5 10
ii=4
Ans(006): -1.28994508297100 3.02749e-03 -7.86716e-10 -3.02782e-03
Ans(006): -1.28994508293748 3.02749e-03 -1.03029e-13 -3.02782e-03
Ans(005): -1.28994508296585 3.02749e-03 -6.65985e-10 -3.02782e-03
Ans(004): -1.28994508326927 3.02749e-03 -7.78754e-09 -3.02783e-03
Ans(004): -1.28994508293751 3.02749e-03 -6.14619e-13 -3.02782e-03
Ans(007): -1.28994508322147 3.02749e-03 -6.66568e-09 -3.02783e-03
Ans(005): 2.70571963551646 2.03107e-03 3.15585e-08 -2.03069e-03
Ans(004): 2.70571963972096 2.03104e-03 3.55271e-15 -2.03073e-03
Ans(004): 2.70571963972096 2.03104e-03 3.55271e-15 -2.03073e-03
Ans(005): 4.58422558613660 -5.05533e-03 1.57708e-06 5.06175e-03
Ans(004): 4.58422544334095 -5.05691e-03 1.37307e-09 5.06017e-03
ii=5
Ans(006): -1.28994508297100 3.02763e-04 -7.86716e-10 -3.02768e-04
Ans(006): -1.28994508293748 3.02764e-04 -1.03029e-13 -3.02767e-04
Ans(005): -1.28994508296585 3.02763e-04 -6.65985e-10 -3.02768e-04
Ans(005): -1.28994508293748 3.02764e-04 0.00000e+00 -3.02767e-04
Ans(004): -1.28994508293751 3.02764e-04 -6.14619e-13 -3.02767e-04
Ans(008): -1.28994508293748 3.02764e-04 0.00000e+00 -3.02767e-04
Ans(006): 2.70571963972096 2.03090e-04 1.77636e-15 -2.03087e-04
Ans(004): 2.70571963972096 2.03090e-04 3.55271e-15 -2.03087e-04
Ans(004): 2.70571963972096 2.03090e-04 3.55271e-15 -2.03087e-04
Ans(006): 4.58422544321654 -5.05838e-04 1.61648e-13 5.05870e-04
Ans(004): 4.58422544334095 -5.05836e-04 1.37307e-09 5.05872e-04
ii=6
Ans(007): -1.28994508293748 3.02766e-05 0.00000e+00 -3.02766e-05
Ans(006): -1.28994508293748 3.02766e-05 -1.03029e-13 -3.02766e-05
Ans(006): -1.28994508293748 3.02766e-05 0.00000e+00 -3.02766e-05
Ans(005): -1.28994508293748 3.02766e-05 0.00000e+00 -3.02766e-05
Ans(004): -1.28994508293751 3.02766e-05 -6.14619e-13 -3.02766e-05
Ans(008): -1.28994508293748 3.02766e-05 0.00000e+00 -3.02766e-05
Ans(006): 2.70571963972096 2.03088e-05 1.77636e-15 -2.03088e-05
Ans(004): 2.70571963972096 2.03088e-05 3.55271e-15 -2.03088e-05
Ans(004): 2.70571963972096 2.03088e-05 3.55271e-15 -2.03088e-05
Ans(006): 4.58422544321654 -5.05852e-05 1.61648e-13 5.05855e-05
Ans(005): 4.58422544321652 -5.05852e-05 -3.55271e-15 5.05855e-05
ii=7
Ans(007): -1.28994508293748 3.02766e-06 0.00000e+00 -3.02766e-06
Ans(006): -1.28994508293748 3.02766e-06 -1.03029e-13 -3.02766e-06
Ans(006): -1.28994508293748 3.02766e-06 0.00000e+00 -3.02766e-06
Ans(005): -1.28994508293748 3.02766e-06 0.00000e+00 -3.02766e-06
Ans(005): -1.28994508293748 3.02766e-06 0.00000e+00 -3.02766e-06
Ans(008): -1.28994508293748 3.02766e-06 0.00000e+00 -3.02766e-06
Ans(006): 2.70571963972096 2.03088e-06 1.77636e-15 -2.03088e-06
Ans(004): 2.70571963972096 2.03088e-06 3.55271e-15 -2.03088e-06
Ans(004): 2.70571963972096 2.03088e-06 3.55271e-15 -2.03088e-06
Ans(006): 4.58422544321654 -5.05854e-06 1.61648e-13 5.05854e-06
Ans(005): 4.58422544321652 -5.05854e-06 -3.55271e-15 5.05854e-06
ii=8
Ans(007): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(007): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(006): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(005): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(005): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(008): -1.28994508293748 3.02766e-07 0.00000e+00 -3.02766e-07
Ans(006): 2.70571963972096 2.03088e-07 1.77636e-15 -2.03088e-07
Ans(005): 2.70571963972096 2.03088e-07 1.77636e-15 -2.03088e-07
Ans(005): 2.70571963972096 2.03088e-07 1.77636e-15 -2.03088e-07
Ans(007): 4.58422544321652 -5.05854e-07 3.55271e-15 5.05854e-07
Ans(005): 4.58422544321652 -5.05854e-07 -3.55271e-15 5.05854e-07
ii=9
Ans(007): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(007): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(006): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(005): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(005): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(008): -1.28994508293748 3.02766e-08 0.00000e+00 -3.02766e-08
Ans(007): 2.70571963972096 2.03088e-08 0.00000e+00 -2.03088e-08
Ans(005): 2.70571963972096 2.03088e-08 1.77636e-15 -2.03088e-08
Ans(005): 2.70571963972096 2.03088e-08 1.77636e-15 -2.03088e-08
Ans(007): 4.58422544321652 -5.05854e-08 3.55271e-15 5.05854e-08
Ans(005): 4.58422544321652 -5.05854e-08 -3.55271e-15 5.05854e-08
ii=10
Ans(007): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(007): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(006): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(005): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(005): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(008): -1.28994508293748 3.02766e-09 0.00000e+00 -3.02766e-09
Ans(007): 2.70571963972096 2.03088e-09 0.00000e+00 -2.03088e-09
Ans(005): 2.70571963972096 2.03088e-09 1.77636e-15 -2.03088e-09
Ans(005): 2.70571963972096 2.03088e-09 1.77636e-15 -2.03088e-09
Ans(007): 4.58422544321652 -5.05853e-09 3.55271e-15 5.05854e-09
Ans(005): 4.58422544321652 -5.05855e-09 -3.55271e-15 5.05853e-09
ii=11
Ans(007): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(007): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(006): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(005): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(005): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(008): -1.28994508293748 3.02766e-10 0.00000e+00 -3.02766e-10
Ans(007): 2.70571963972096 2.03086e-10 0.00000e+00 -2.03091e-10
Ans(005): 2.70571963972096 2.03089e-10 1.77636e-15 -2.03087e-10
Ans(005): 2.70571963972096 2.03089e-10 1.77636e-15 -2.03087e-10
Ans(007): 4.58422544321652 -5.05850e-10 3.55271e-15 5.05858e-10
Ans(006): 4.58422544321652 -5.05864e-10 -3.55271e-15 5.05850e-10
ii=12
Ans(007): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(007): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(006): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(005): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(005): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(008): -1.28994508293748 3.02762e-11 0.00000e+00 -3.02798e-11
Ans(007): 2.70571963972096 2.03055e-11 0.00000e+00 -2.03144e-11
Ans(005): 2.70571963972096 2.03073e-11 1.77636e-15 -2.03109e-11
Ans(005): 2.70571963972096 2.03073e-11 1.77636e-15 -2.03109e-11
Ans(007): 4.58422544321652 -5.05764e-11 3.55271e-15 5.05924e-11
Ans(006): 4.58422544321652 -5.05906e-11 -3.55271e-15 5.05835e-11
ii=13
Ans(007): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(007): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(006): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(005): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(005): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(008): -1.28994508293748 3.02691e-12 0.00000e+00 -3.02336e-12
Ans(007): 2.70571963972096 2.02505e-12 0.00000e+00 -2.03215e-12
Ans(005): 2.70571963972096 2.03215e-12 1.77636e-15 -2.03215e-12
Ans(005): 2.70571963972096 2.03215e-12 1.77636e-15 -2.03215e-12
Ans(007): 4.58422544321652 -5.05196e-12 3.55271e-15 5.06084e-12
Ans(006): 4.58422544321652 -5.06262e-12 -3.55271e-15 5.05196e-12
ii=14
Ans(007): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(007): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(006): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(005): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(005): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(008): -1.28994508293748 3.03757e-13 0.00000e+00 -3.01981e-13
Ans(007): 2.70571963972096 2.02505e-13 0.00000e+00 -2.06057e-13
Ans(005): 2.70571963972096 2.04281e-13 1.77636e-15 -2.02505e-13
Ans(005): 2.70571963972096 2.04281e-13 1.77636e-15 -2.02505e-13
Ans(007): 4.58422544321652 -5.08038e-13 3.55271e-15 5.15143e-13
Ans(006): 4.58422544321652 -5.15143e-13 -3.55271e-15 5.06262e-13
たしかに変化は見られるのだけれど、ニュートン法の収束が早すぎるのか、ほとんど変化が見られない。もっと桁数の多い計算をできるようにして試してみたい。
