バーラル @ ウィキ

Cで方程式の数値解法


※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

はさみ撃ちによる方程式の解を求めるプログラミング

ただし、UNIXで書いたので、若干違うところがあるので注意!

#include plugin Error : 指定されたページはございません。
<stdio.h>
#include plugin Error : 指定されたページはございません。
<math.h>
#include plugin Error : 指定されたページはございません。
<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
  • e with arg 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




考察


★自分の与えた入力データの正当性 (設定した関数に対して何故そのような入力データで良いのかを論ぜよ)
求める方程式を手計算で一階微分すると、極値が 2+√3)、2ー√3となっていて、そのときの方程式の値が、(およそ3.732と0.268とすると)ー4.39と16.39になる。よってこのあいだにひとつ解があり,その外に2つ解があることが予想され,おおよそx=10のとき方程式の値は正(y=446),x=ー10のとき方程式の値は負(ー1614)となる。
求める方程式は三次方程式なので、以上の条件より、ー10から10の値に3つの解があると推測される。ステップ数は、極値の間を飛ばさない程度の値(xを1刻み)にしたいので、20とした。

★自分の与えた収束条件の正当性と、解に期待できる精度
収束条件を以下の様にした。
fabs(c-a) < fabs(c)*e_pow||fabs(c-b) < fabs(c)*e_pow
e_powはここではe-10とした。
このとき、cのスケール(例えばC=10000000000)かけるe_pow、に対してfabs(c-a)がおおきければ(例えばfabs(c-a)=100000000)、さらにはさみ撃つので、(Cのスケールのe-10分の1)程度の幅でXの精度を求められると期待でき、収束条件として適当ではないかと思う。
または(||)で区切っているのは、凸のグラフと凹のグラフに対応するためである(aかbは固定される)。


★自分の得た解が実際に備えている精度
求めた解の値は、
Ans(017): -1.28994508292641 2.597140e-10
Ans(011): 2.70571963974350 -1.692015e-10
Ans(016): 4.58422544313894 -8.560406e-10
である。
上の考察により、およそ小数点以下9桁までが正しい精度だと考えられる。なぜなら、(Cのスケールのe-10分の1)程度の幅が予想されるので、小数点第10桁目はこの幅の影響があると思うからである。
Ans(017): -1.289945082
Ans(011): 2.705719639
Ans(016): 4.584225443
までが実際に備えている精度である。

実際上の結果でC値の変化を見てみると、小数点第9桁目までは変化が終了している様に見える。


★解の持つ精度を入力から制御する方法 (可能なら各自のプログラムに組み込め)
組み込んだつもりです。

★二分法、ニュートン法と比較した収束の様子に関する考察
二分法は、一回の再帰呼出しで2進数表示で1桁づつ解の精度があがる。(二分法という手法だからこそ)
実際実行例を見るとそうなっているはず。
ニュートン法は、解にちかい初期値から始めた場合、一回の再帰呼出しで2次の収束になると言われている。

収束の様子をグラフにしようと思ったが、諸事情によりあきらめた。
最後にいちおう、ニュートン法の収束の様子を観察するためのプログラムとその結果を添付する。

★二分法、ニュートン法との使い分けに関する考察 (ガイドラインを提示せよ)
1.微分可能である。
2.導関数が求められる。
3.ある程度解の範囲がわかっている。(どこの解を示すかわからないため)
この3つが満たせるなら、収束が非常にはやいニュートン法を使うべきである。
これを1つでも満たせないなら、二分法を使うべきである。

★おまけ
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

たしかに変化は見られるのだけれど、ニュートン法の収束が早すぎるのか、ほとんど変化が見られない。もっと桁数の多い計算をできるようにして試してみたい。