今週のプログラムその5

解説

前回使った二分法のプログラムは、同じ関数を三回も書かなければならないという、なかなかダサいものだった。 そこで、今回はそれを一回ですませるために、Cでは"関数"、FORTRANでは"関数副プログラム"というものを導入しよう。

sin,cos,tanなど、よく使われる数学関数は、あらかじめ用意されている。 たとえば、Cなら

#include<stdio.h>
#include<math.h> 

main()  
{ double x;
  x = 0.1;
  printf("%f",sin(x));
}

      program main
      real*8 x
      x = 0.1
      write(*,*)sin(x)
      end

自分が作った関数に名前をつけて、sinやcosと同じように使えれば便利だろう。 たとえば、func(x)という関数を作ったとする。 これを使うには、Cでは

#include<stdio.h>

double func(double);

main()  
{ double x;
  x = 0.1;
  printf("%f",func(x));
}

      program main
      real*8 x
      real*8 func
      x = 0.1
      write(*,*)func(x)
      end

さて、あとは、関数本体を作らなくてはならない。これは、基本的にはmainプログラムとは独立した一個のプログラムだが、mainの後ろにつけ加えておけばよい。 別のファイルにする方法については説明しないので、興味のある人は調べてください。

前回の二分法に使った関数をfuncという名前にすることを考えよう。まずCでは

double func(double x)
{
  return pow(x,3) + 5 * x*x - 3*x - 5;
}

double func(double x)
{
  double f;
  f = pow(x,3) + 5 * x*x - 3*x - 5;
  return f;
}

      function func(x)
      real*8 func                        
      real*8 x 
      func= x**3 + 3*x**2 + 5*x + 1
      end 

      function func(x)
      real*8 func                        
      real*8 x,v
      f= x**3 + 3*x**2 + 5*x + 1
      func = f
      end 

なお、関数(関数副プログラム)中に現われる変数名は"局所的"である。つまり、 関数の中でxという名前の変数を使っていても、それはmainプログラム中に現われる同じ名前の変数xとは関係ない。 逆に、関数を定義しているところで変数名がxになっているからといって、 mainプログラム中でも同じ名前にしなくてはならないわけではない。 要するに、関数の中でどういう名前の変数が使われているかを気にせずに、sinやcos などと同じように使ってよい、ということ。

なお、関数(関数副プログラム)はmainと同様にプログラムなので、forやifといった制御もできるし、また、変数として配列を使うこともできる。 そのあたりについては次回考えることにしよう。

ちなみに、program12のニュートン法のプログラムは、収束の判定条件を手抜きしているため、解が0に近い場合にうまくいかない。 0に近い場合もうまく扱えるようにプログラムを工夫してみるべし

• C版
• FORTRAN版

まずはC版

program

/* とりあえず、関数の使い方を見るだけのために、
解説に使ったプログラムの完全版を*/
#include<stdio.h>                       いつものおまじない
#include<math.h>              数学関数を使うおまじない

double func(double);         倍精度関数funcを使うというプロトタイプ宣言

main()                      mainプログラムの始まり
{ 
  double x;
  x = 0.1;
  printf("x=%f  f=%f\n",x,func(x));     funcを使ってみる
  x = 0.2;
  printf("x=%f  f=%f\n",x,func(x));     念のため、変数の値を変えて、もう一度funcを使う
}                           mainプログラムはここまで

double func(double x)                    関数funcの始まり
{
  return pow(x,3) + 5 * x*x - 3*x - 5;     値を計算して、mainに戻る
}                                       関数funcの終わり

program 12

/* 関数の利用例として 
   ニュートン法による非線形方程式の解を求める */
#include<stdio.h>           いつものおまじない
#include<math.h>              数学関数を使うおまじない
#define NMAX 30               プログラム中のNMAXをすべて30に

double func(double);          倍精度関数funcを使うというプロトタイプ宣言
double dfunc(double);         倍精度関数dfuncを使うというプロトタイプ宣言
main()
{
  double x,xtemp;
  int irepeat;

/* 解探索の初期値 */
  scanf("%lf",&x);

/* Newton法を最大nmax回くりかえす */
  for(irepeat=1; irepeat<=NMAX; ++irepeat){

/* 導関数の値が小さすぎるときは、精度が出ないので計算終了 */
    if(fabs(dfunc(x)) < 1e-10){
       printf("too small derivative at x=%17.12f\n",x);
       return;
    }

/* 導関数の値が大丈夫なら解の候補を計算 */
     else {
        xtemp = x - func(x)/dfunc(x);
     }

/*  xとxtempの差が与えた精度以内なら終了 */
     if(fabs((xtemp-x)/x) <= 1e-13){
       printf("converged: x= %17.12f\n",xtemp);
       return;
     }
     else {
       x=xtemp;
       printf("%d step x=%17.12f\n",irepeat,xtemp);
     }
   }

/*  nmax回のくりかえしで解に収束しなかったときは終了 */
   printf("not converged after %d steps\n",NMAX);
   printf("current value of x is %17.12f\n",xtemp);
}                                  ここまでがmainプログラム

/*  解を得たい関数 */
double func(double x)               関数funcの定義
{
 return x+cos(x);
}                                   ここまで
 
/*  解を得たい関数の導関数 */
double dfunc(double x)              関数副dfuncの定義
{
  return 1-sin(x);
}                                  ここまで

FORTRAN版

program

C とりあえず、関数の使い方を見るだけのために、
C 解説に使ったプログラムの完全版を
      program main                     mainプログラムの始まり
      real*8 x
      real*8 func                      funcが倍精度関数であることを宣言   
      x = 0.1
      write(*,*)"x=",x,"  f=",func(x)               funcを使ってみる
      x = 0.2
      write(*,*)"x=",x,"  f=",func(x)        念のため、変数の値を変えて、もう一度funcを使う
      end                         mainプログラムはここまで

      function func(x)                   関数funcの始まり
      real*8 func                        funcが倍精度であることを宣言          
      real*8 x                       funcが受け取る変数xが倍精度であることを宣言
      func=  x**3 + 5 * x*x - 3*x - 5  ここで関数値を計算
      end                               funcはここまで

program 12

C 関数副プログラムの利用例として
C ニュートン法による非線形方程式の解
      program main
      parameter(nmax=30)        主プログラム中のnmaxの値を30に(ここでの定義は主プログラム中でのみ有効。関数副プログラム中にnmaxという変数があっても影響をうけない)
      real*8 x,xtemp
      real*8 func,dfunc        倍精度関数を使う宣言

C 解探索の初期値
      read(*,*)x

C Newton法を最大nmax回くりかえす
      do irepeat = 1,nmax

C 導関数の値が小さすぎるときは、精度が出ないので計算終了
        if(abs(dfunc(x)).lt.1e-10)then
          write(*,*)"too small derivative at x= ",x
          stop

C 導関数の値が大丈夫なら解の候補を計算
        else
          xtemp = x - func(x)/dfunc(x)
        endif

C  xとxtempの差が与えた精度以内なら終了
        if(abs((xtemp-x)/x).le. 1e-13)then
          write(*,*)"converged: x= ",xtemp
          stop
        else
          x=xtemp
          write(*,*)irepeat," step  x=",xtemp
        endif
      enddo

C  nmax回のくりかえしで解に収束しなかったときは終了
      write(*,*)"not converged after ",nmax," steps"
      write(*,*)"current value of x is ",xtemp
      end                           mainプログラムはここまで

C  解を得たい関数
      function func(x)                関数副プログラムfuncの定義
      real*8 func                     funcが倍精度関数であることを宣言
      real*8 x                        xが倍精度変数であることを宣言
      func= x+cos(x)                funcの値を計算
      end                           関数funcはここまで

C  解を得たい関数の導関数
      function dfunc(x)                関数副プログラムdfuncの定義
      real*8 dfunc
      real*8 x
      dfunc= 1-sin(x)
      end                            関数dfuncはここまで