今週のプログラムその4

今週はif文の簡単な使い方とその応用として非線形方程式を二分法で解くプログラムを作ろう

ポイント

実数(浮動小数点)演算の精度は計算機の語長(数字ひとつが何ビットで表されるか)で決まる。一語32ビットの機械の場合、実数は普通10進数で7桁程度の精度がある。しかし、計算によってはこれでは精度が不足する場合があるので、数字ひとつを2語で表す(32ビットの機械なら64ビットで)ように指定する倍精度実数を使うことができる。たとえば、数aと配列b(要素数10)を倍精度数として宣言するには
double a,b[10]; (Cの場合） real*8 a,b(10) (FORTRANの場合)
などとする。昔と違って今の計算機は倍精度実数の計算だからといって特に計算時間がかかるわけではないので、科学計算の場合は常に倍精度を指定すべきだろう。
数字を0.000000000123などと書いていては大変なので、1.23e-10という書きかたができる。つまり10^-10をe-10と表記するわけだ。
if文は条件によってプログラムの実行範囲を変えるための文である。
Cでは
if(条件1){ 実行部分1 } else if(条件2){ 実行部分2 } else{ 実行部分3 }
という形式で書く。
この場合、条件1が満たされれば実行部分1が実行される。条件1が満たされず条件2が満たされれば実行部分2が実行される。条件1も2も満たされなければ実行部分3が実行される。
実行部分は一文だけのプログラムでも長いプログラムでもかまわない。また
else if(条件){ 実行部分 }
はひとつとは限らず、異なる条件をいくつ並べてもよい。その場合、条件判断は上から順に行なわれ、最初に満たされた条件に対応する実行部分が実行される。else以下はどの条件も満たされなかった場合に実行される文である。
なお、if(条件){ 実行部分 }だけが必須であり、else ifとelseに関する部分はなくても構わない。
条件を記述する関係演算子には、値が等しい(==)、等しくない(!=)、左辺が右辺より大きい(>)、小さい(<)、左辺が右辺以上(>=)、以下(<=)がある。複数の条件の論理和(||)・論理積(&&)等も使える。条件の詳細は参考書にあたること。

同じプログラムをFORTRANで書くと
if(条件1)then 実行部分1 else if(条件2)then 実行部分2 else 実行部分3 endif
となる。Cとほとんど同じ形式だが、thenを書く必要があることと最後を endifで終える点に注意。Cと同様
else if(条件)then 実行部分
はいくつあってもいいし、ひとつもなくてもいい。elseの部分もなくてもよい。ただし、どんな場合でもendifは必須。
条件を記述する関係演算子には、値が等しい(.eq.)、等しくない(.ne.)、左辺が右辺より大きい(.gt.)、小さい(.lt.)、左辺が右辺以上(.ge.)、以下(.le.)がある。複数の条件の論理和(.or.)・論理積(.and.)等も使える。詳細は参考書にあたること。

program9は少々構造が複雑である。構造だけを取り出すと、Cなら


一番外側がfor
その中に、
  if(条件1){実行1}
  else if(条件2){実行2}   
  else {実行3}
という構造があり、さらに実行2の中に
  if(さらに条件1){さらに実行1}
  else{さらに実行2}
という構造を含んでいる。

というforと二段階のifの入れ子構造になっていることを確認してほしい。

FORTRANで表すと


一番外側がdo
その中に、
  if(条件1)then 実行1
  else if(条件2)then 実行2   
  else 実行3
  endif
という構造があり、さらに実行2の中に
  if(さらに条件1)then さらに実行1
  else さらに実行2
  endif
という構造を含んでいる。

となる。

Cのプログラム10中ではべき乗を計算する関数powを使っている。FORTRANではx**yと書けば"xのy乗"を表すのに対し、Cではpowを用いてpow(x,y)としなくてはならない。ただし、xとyは倍精度実数である。したがって、厳密には、2乗をpow(x,2)ではなく(これでは2が整数とみなされる)、 pow(x,2.)やpow(x,2e0)などと書く必要がある(ただし、どうやら、pow(x,2)でもちゃんと計算してくれる機械もあるようだ)。もっとも、 Cの場合、xの二乗くらいならx*xと書くほうが簡単だし、見やすいと思う。他にどういう関数があるかは解説書を参照のこと。

なお、Cではこういった数学関数を使うために#include<math.h>というおまじないが最初に必要。また、コンパイル時には
#cc -o Prog prog.c -lm
と数学関数ライブラリを使うための宣言-lmが必要になる (OSや環境の設定によっては、なくてもいい場合もある)。
Cには大きくわけてK&RとANSIの二種類があり、古いCコンパイラでは何も指定しないとK&Rとみなされる場合がある。例題はすべてANSI準拠なので(といっても、前回までのプログラムはどちらでも同じ)、コンパイラによってはANSI Cであることを宣言する必要がある。例えば
cc -Aa -o prog prog.c -lm
など(ここで -Aa がANSI Cを宣言するオプション)。コンパイル時に警告メッセージが出ないなら気にする必要はない。
Cプログラムでは
printf("solution is between %17.12f and %17.12f\n",x1,x2)
のように、%fではなく%12.7fと書いたところがある。これは出力のしかたを指定するもので、今の例では実数ひとつを17桁分の幅の中に書き、小数点以下12桁を表示するという意味。17桁には小数点やマイナス記号も含まれる。この手の書式指定は自習してもらうつもりだったが、%fではあまりにも出力される桁数が少ない場合があるので、例外的に使用した。

FORTRANでも同じことはできるが、こちらはwrite(*,*)でも充分な桁数を出力する場合が多いようなので、あえて解説しない。

ところで、program10では、解を求めたい関数がプログラム中の三か所に書かれている。これはエレガントではないし、別の関数を解きたい場合に三か所書き直さなくてはならないので、間違う可能性も高い。

関数を一回だけ書いてすます方法として、Cならdefine文を使う手があるだろう

#define FUNC pow(x,3)+5*x*x-3*x-5
と書けば、プログラム中のFUNCはすべて目的の関数に置き換えられるはずだ。ただし、変数名がすべてxになるので、x1やx2といった変数をいちいちxに代入してから関数を使うことになる。余裕がある人は試してみてください(実は、defineのもっと高度な使い方をすれば、本当の関数のように定義することもできる。興味のある人は調べて下さい)。

もっとも、この手はエレガントというよりトリッキーなものであって、あまり薦めない。正しくは、"関数"を定義するべきであり、それについては次回以降考えることにする。

FORTRANではこのようなdefineの使い方に相当するものが用意されていない。文関数というものを使えば、もっときれいに関数が定義できるので、興味のある人は調べて下さい。ただし、FORTRANでも"関数副プログラム"を用いるのが一番よい手段であり、それについては次回以降考えることにする。

なお、プログラム自体は

まずはC版

program 9


/* ifの使い方を知るためだけの、殆ど無意味なプログラム
  (何をやってるか、考えてみてください) */
#include<stdio.h>
main()
{
 int i;       整数変数iを宣言
 double x;    倍精度変数xを宣言
 for(i=1; i<=10; ++i){   以下を10回繰り返し
  scanf("%lf",&x);       倍精度変数xをキーボードから読み込む。倍精度数に読み込むときは"%f"ではなく"%lf"とする(lはlongの意味)
  if( x > 10 ){                        x > 10のとき
     printf("%f > 10\n",x); }    
  else if( x > 5 ){                    x <= 10でかつx > 5のとき
     printf("10 >= %f > 5\n",x); }  
  else if( x > 0 ){                     x <= 5(x <= 10は当然成立している)でかつx > 0のとき
     printf("5 >= %f > 0\n",x); }   
  else {                                 上のどの条件も満たさない(つまりxは負)とき
     printf("%f <= 0\n",x); 
     return;        プログラム終了
    }               else以下の実行はここまで(if--else if--else if--else という構造もここで終わり)
 }              forループの終了
}              プログラムの終わり

program 10


/* if文の利用例として2分法で非線形方程式の解を求める。
 下の例題では
 x^3 + 5x^2 - 3x -5 = 0
 の解を求めている(^はべき乗を表す)。
 */
#include<stdio.h>
#include<math.h>               数学関数を使うおまじない
#define NMAX 50                プログラム中のNMAXをすべて50に置き換え
main()                     ここから主プログラム
{
  double x1,x2,x3;       倍精度変数x1,x2,x3の宣言
  double f1,f2,f3;          倍精度変数f1,f2,f3の宣言(上のx1,x2,x3,f1,f2,f3をひとつの宣言文で宣言してもよいし、いくつかに分けてもよい。ここではx1,x2,x3とf1,f2,f3のふたつの組にわけて宣言した。意味は同じなので、プログラムを見やすくするためだけ)
  int irepeat;              整数変数irepeatの宣言

/* 解探索の範囲を指定 */
  scanf("%lf%lf",&x1,&x2);    ふたつの倍精度数をキーボードから読み込む。倍精度数に読み込むときは"%f"ではなく"%lf"とする(lはlongの意味)

/* 2分法を最大NMAX回くりかえす */
  for(irepeat=1; irepeat<=NMAX; ++irepeat){
    f1 = pow(x1,3.) + 5*x1*x1 - 3*x1 -5;     x1での関数値(pow(x,y)はxのy乗を計算する関数)
    f2 = pow(x2,3.) + 5*x2*x2 - 3*x2 -5;     x2での関数値

/*  x1とx2の差が与えた精度以内なら終了 */
    if( fabs((x1-x2)/x1) <= 1e-10){       x1とx2の相対差が1e-10以下なら収束とみなす 
       printf("convergence\n");       convergenceという文字を出力
       printf("x=%17.12f f=%17.12f\n",x1,f1);    x1の値とそこでの関数値f1を出力
       printf("x=%17.12f f=%17.12f\n",x2,f2);    x2の値とそこでの関数値f2を出力

       return;                        プログラムの実行を終了(mainプログラム中のreturn文はプログラムを終了する)
    }

/*  収束してないと判断されたとき。
  f1とf2の符号が逆なら範囲を二分する */
    else if(f1*f2 < 0){         f1とf2の積が負なら
       x3 = (x1+x2)/2;              x1とx2の中点の値をx3に
       f3 = pow(x3,3.) + 5*x3*x3 - 3*x3 -5;     x2での関数値

       if(f1*f3 < 0){           f1とf3の積が負なら
         x2=x3;                     x1はそのままで、x3を新しいx2とする
       }
       else{                           f1とf3の積が負でなければ  
         x1=x3;                     x2はそのままで、x3を新しいx1にする
       }                          elseの終わり(if--else構造がここで終わる)
       printf("%d step x1= %17.12f f1= %17.12f  x2= %17.12f  f2= %17.12f\n",irepeat,x1,f1,x2,f2);         ステップ数、x1,f1,x2,f2の値を出力
    }

/*  f1とf2の積が正なら、この範囲に解がないと判断して終了 */
    else{
      printf("no solution is found between %17.12f and %17.12f\n",x1,x2);
      return;                       プログラムの実行を終了
    }                      elseの終わり(if--else if--else構造がここで終わる)
  }                                 ここまでがforループの範囲

/*  NMAX回のくりかえしで解に収束しなかったときは終了 */
  printf("no convergence after %d steps\n",NMAX);

}   プログラムの終わり

FORTRAN版

program 9


C ifの使い方を知るためだけの、殆ど無意味なプログラム 
C  (何をやってるか、考えてみてください) */
      program main
      real*8 x              倍精度変数xを宣言

      do i=1,10               これ以下を10回繰り返し
        read(*,*)x            倍精度変数xをキーボードから読み込む。
        if( x.gt.10)then       x > 10のとき
            write(*,*)x," > 10"
        else if(x.gt.5)then    x <= 10でかつx > 5のとき
            write(*,*)"10 >= ",x," > 5"
        else if(x.gt.0)then    x <= 5(x <= 10は当然成立している)でかつx > 0のとき
            write(*,*)"5 >= ",x," > 0"
        else                   上のどの条件も満たさない(つまりxは負)とき
            write(*,*)x," <= 0 "
            stop               プログラム終了(elseの場合、この二行が実行される)
        endif                  if--else if--else if--else という構造の終わり
      enddo                    forループの終了
      end                      プログラム終了

program 10


C if文の利用例として2分法で非線形方程式の解を求める。
C 下の例題では
C x^3 + 5x^2 - 3x -5 = 0
C の解を求めている(^はべき乗を表す)。

      program main
      parameter(nmax=50)       プログラム中のnmaxの値をすべて50に
      real*8 x1,x2,x3       倍精度変数x1,x2,x3の宣言
      real*8 f1,f2,f3       倍精度変数f1,f2,f3の宣言(上のx1,x2,x3,f1,f2,f3をひとつの宣言文で宣言してもよいし、いくつかに分けてもよい。ここではx1,x2,x3とf1,f2,f3のふたつの組にわけて宣言した。意味は同じなので、プログラムを見やすくするためだけ)

C 解探索の範囲を指定
      read(*,*)x1,x2        ふたつの倍精度数をキーボードから読み込む。

C 2分法を最大nmax回くりかえす
      do irepeat = 1,nmax
        f1 = x1**3 + 5*x1**2 - 3*x1 -5    x1での関数値
        f2 = x2**3 + 5*x2**2 - 3*x2 -5    x2での関数値

C  x1とx2の差が与えた精度以内なら終了
        if(abs((x1-x2)/x1).le. 1e-10)then   x1とx2の相対差が1e-10以下なら収束とみなす 
          write(*,*)"convergence"      convergenceという文字を出力
          write(*,*)"x1 = ",x1," f1 = ",f1   x1の値とそこでの関数値f1を出力
          write(*,*)"x2 = ",x2," f2 = ",f2   x2の値とそこでの関数値f2を出力
          stop                             プログラムの実行を終了

C  収束してないと判断されたとき。
C  f1とf2の符号が逆なら範囲を二分する
        else if(f1*f2.lt.0)then            f1とf2の積が負なら
            x3 = (x1+x2)/2                 x1とx2の中点の値をx3に
            f3 = x3**3 + 5*x3**2 - 3*x3 -5    x3での関数値
            if(f1*f3.lt.0)then            f1とf3の積が負なら
              x2=x3                       x1はそのままで、x3を新しいx2にする
            else                          f1とf3の積が負でなければ

              x1=x3                       x2はそのままで、x3を新しいx1にする

            endif                    if--else構造の終わり
       write(*,*)irepeat," step  x1=",x1,"  f1=",f1,"  x2=",x2,"  f2=",f2    ステップ数、x1,f1,x2,f2の値を出力

C  f1とf2の積が正なら、この範囲に解がないと判断して終了
        else
          write(*,*)"no solution is found between ",x1," and ",x2
          stop                            プログラムの実行を終了
        endif                           if--else if--else構造の終わり
      enddo                               ここまでがdoループの範囲

C  nmax回のくりかえしで解に収束しなかったときは終了
      write(*,*)"not converged after ",nmax," steps"
      end                                 プログラムの終わり