next up previous
Next: 4 レポート課題 Up: 計算天文学 II Previous: 2 今日の講義

3 プログラム

今回は Fortran と C++ と両方出してみよう。初期条件は のところ にだけ値がある 関数的なものということにする。

なお、面倒なので D=1 ということのする。時間の単位が と思えば同 じことである。

pt

c
c     program parabolic1
c
c     Copyright 2000 J. Makino

      subroutine initparam(u, nx, dx, dt, niter, gmode)
      real*8 u(*), dx, dt, tmax
      integer nx, niter, gmode, i
      write(*,*)'Enter nx, tmax, dt:'
      read(5,*) nx, tmax, dt
      dx = 1d0/nx
      niter = (tmax+dt/2)/dt
      write(*,*)'Enter graphics mode 0: no G'
      write(*,*)'                    1: animation'
      write(*,*)'                    2: no-erase)'
      read(5,*) gmode
      write(6,*) 'nx = ', nx, ' dt = ', dt, ' niter = ', niter
      write(6,*) 'gmode = ', gmode
      do i = 1, nx+1
         u(i) = 0
      enddo
      u((nx+2)/2) = nx
      end

      subroutine push_system(u, unew, nx, dx, dt)
      real*8 u(*), unew(*), dx, dt, lambda
      integer nx, i
      lambda = dt/(dx*dx)
      do i = 2, nx
         unew(i) = u(i) + lambda*(u(i-1)-2*u(i)+u(i+1))
      enddo
      do i = 2, nx
         u(i) = unew(i)
      enddo
      end

      subroutine initgraph
      integer pgopen
      if(pgopen('?') .LE. 0) stop
      call pgask(.false.)
      end
      
      subroutine show_graphics(u, nx, dx, gmode,
     $     xmin, xmax, ymin,ymax, first)
      real*8 u(*), dx
      real*8 xmin, xmax, ymin, ymax
      integer nx, gmode, i, first
      if (gmode .eq. 0) return
      call pgbbuf
      if (first .eq. 1) then
         call pgpage
         call pgenv(real(xmin), real(xmax), real(ymin), real(ymax),0,-2)
      endif
      if (gmode .eq. 1) then
         call pgeras
      endif
      call pgbox('bcnst', 0.0, 0, 'bcnst', 0.0, 0)
      call pglab('X', 'Y', ' ')
      call pgmove(0.0, real(u(1)))
      do i = 2, nx+1
         call pgdraw(real((i-1)*dx), real(u(i)))
      enddo
      call pgebuf
      end
      
      program parabolic
      integer nmax
      parameter (nmax=10001)
      real*8 u(nmax), unew(nmax), dx, dt
      integer nx, niter, gmode, iter
      call initparam(u, nx, dx, dt, niter, gmode)
      if (gmode .ne. 0) call initgraph
      call show_graphics(u, nx, dx, gmode,
     $        0d0, 1d0, 0d0, 10d0, 1)
      do iter = 1, niter
         call push_system(u, unew, nx, dx, dt)
         call show_graphics(u, nx, dx, gmode,
     $        0d0, 1d0, 0d0, 10d0, 0)
      enddo
      end

これを適当なファイルに落して、コンパイル・実行する。コンパイルには

        f77  -o parabolic1  parabolic1.f -L/home/makino/utils/pgplot5.2 \
                -lpgplot -L/usr/X11R6/lib -lX11
とすると、 parabolic1 という名前の実行ファイルができる。

ここでは、 PGPLOT というわりと広く使われているグラフィックライブラリを 使う。これは画面、 PSファイル等様々な出力形態を同じプログラムで切替え られるので割合に便利である。

中身の解説は、見ての通りという気もする。次に C++ のプログラムをつける

    1   //     program parabolic1
    2   //
    3   //     Copyright 2000 J. Makino
    4   #include "cpgplot.h"
    5   #include  <stdlib.h>
    6   #include  <math.h>
    7   #include  <stdiostream.h>
    8	
    9   void initparam(double  u[],
   10                int & nx,
   11                double &dx,
   12                double & dt,
   13                int & niter,
   14                int &gmode)
   15   {
   16       cerr << "Enter nx, tmax, dt:";
   17       double tmax;
   18       cin >> nx >> tmax >> dt;
   19       dx = 1.0/nx;
   20       niter = (tmax+dt/2)/dt;
   21       cerr << "Enter graphics mode 0: no G\n"
   22            << "                    1: animation\n"
   23            << "                    2: no-erase)\n";
   24       cin >> gmode;
   25       cout << "nx = " << nx<< " dt = " << dt << " niter = "<< niter<< "\n";
   26       cout << "gmode = " <<  gmode << "\n";
   27       int i;
   28       for (i=0;i<=nx; i++) u[i] = 0;
   29       u[(nx+1)/2] = nx;
   30   }
   31        
   32   void push_system(double u[],
   33                    double unew[],
   34                    int nx,
   35                    double dx,
   36                    double dt)
   37   {
   38         double lambda = dt/(dx*dx);
   39         int i;
   40         for (i=1; i<nx; i++){
   41             unew[i] = u[i] + lambda*(u[i-1]-2*u[i]+u[i+1]);
   42         }
   43         for(i=1; i<nx; i++){
   44            u[i] = unew[i];
   45         }
   46   }
   47        
   48   void initgraph()
   49   {
   50       if(cpgopen("?") !=1 ) exit(-1);
   51       cpgask(0);
   52   }
   53        
   54   void show_graphics(double u[],
   55                      int nx,
   56                      double dx,
   57                      int gmode,
   58                      double xmin,
   59                      double xmax,
   60                      double ymin,
   61                      double ymax,
   62                      int first)
   63   {
   64       if (gmode == 0) return;
   65       cpgbbuf();
   66       if (first == 1){
   67            cpgpage();
   68            cpgenv(xmin, xmax, ymin, ymax,0,-2);
   69       }
   70       if (gmode == 1) {
   71            cpgeras();
   72       }
   73       cpgbox("bcnst", 0.0, 0, "bcnst", 0.0, 0);
   74       cpglab("X", "Y", " ");
   75       cpgmove(0.0, u[0]);
   76       int i;
   77       for(i=1;i<nx+1;i++){
   78           cpgdraw(i*dx, u[i]);
   79       }
   80       cpgebuf();
   81   }
   82   
   83   void main()
   84   {
   85   #define NMAX 10001
   86       static double u[NMAX], unew[NMAX];
   87       double dx, dt;
   88       int nx, niter, gmode, iter;
   89       initparam(u, nx, dx, dt, niter, gmode);
   90       if (gmode !=  0) initgraph();
   91       show_graphics(u, nx, dx, gmode, 0.0, 1.0, 0.0, 10.0, 1);
   92       for(iter = 1; iter <=niter; iter++){
   93            push_system(u, unew, nx, dx, dt);
   94            show_graphics(u, nx, dx, gmode,0.0, 1.0, 0.0, 10.0, 0);
   95       }
   96   }

コンパイルには

g++ -O2  -I/home/makino/utils/pgplot5.2  -o parabolic1   parabolic1.C \
            -L/home/makino/utils/pgplot5.2 -lcpgplot -lpgplot \
            -L/usr/X11R6/lib -lX11 -lf2c -lm

とりあえず、ここ数回で使う C++ の機能はほぼこのプログラムで網羅されて いるので、頭から見ていくことにする。

``//'' で始まる行はコメントである。 Fortran での ``C'' と同じ。

#include "cpgplot.h" ``#'' で始まる行は、「プリプロセッサコマ ンド」というものである。これにはいくつかあるが、ここで使っている #include というのは、次の行に進む前に指定した名前のファイルを読み込 めという指示である。読み込まれるファイルの中身はなんでもいいが、ここで は PGPLOT ライブラリが準備したいろいろな定義が書かれたファイル cpgplot.h を読み込む(インクルードする)。

インクルードされるファイルをインクルードファイル、あるいはヘッダ (header) ファイルという。これのありかは、 UNIX システムの場合は通常 /usr/include であるが、コンパイラに与えるオプション -I に よって捜す場所の追加を指定できる。

5-7 行目は C++ の標準的な機能を使うための指定である。

9-30行目が最初の「関数」 initparam である。ここでは Fortran のも のをほとんど逐語的に変換している。 Fortran との違いを説明していくと、

と、ここまでは形式的な違いだが、次はちょっと意味的な違いがあるので注意 したい。

呼ばれた関数の中で値を変更し、その結果を呼んだほうの関数に反映さ せたい時には、型名と引数名の間に & をつける。これをつけないと、 関数が呼ばれる時に値が呼ばれたほうが持っている変数にコピーされ、そっち が書き変わっても呼んだほうには影響しない。

C言語では上の機能は存在しなかった。そのかわりに、変数の「アドレス」を 扱う演算が可能である。例えば変数 x に入力する関数 readx(x) は、 Fortran なら

      subroutine readx(x)
      real * 8 x
      read(5,*) x
      write(6,*) 'x = ', x
      end
      .....
      call readx(x)
      .....
となるが、これが C だとこんな感じになる。
void readx(double * xptr)
{
    scanf("%lf", xptr);
printf("x = %le\n", *xptr);
}
.....
readx(&x);
.....
&x x のアドレスであ。 C では(C++ でも)任意の型の変数の 「アドレス」である型(ポインタ)が使えるようになっていて、例えば int 型の変数へのポインタは int * ptr とかいう風に宣言できる。で、そ の変数が指している場所の値というのが、 *ptr で使える。

なんだかよくわからないと思うが、ややこしいだけなので落ちついて考えれば わかる。この講義ではなるべく ポインタとかは避けたいと思うが、 C/C++ 言語の柔軟性の多くがこのポインタの 自由な扱いから来るということも確かである。人のプログラムを見る時には注 意が必要になる。

で、上のポインタで用が足りるのは確かだが、上のようなプログラムは読みや すくないのもまた確かである。また、 * & を付け忘れるといった 間違いもおきやすい。ということで、 C++ では関数宣言のほうで & を つけることで、関数に渡ってきたものを実質的にはポインタとして扱うがプログラ ムとしてはそれを意識しないで書けるようになっている。

なお、 Fortran ではそもそもこの当たりが曖昧で、上の readx

             call readx(5.0d0)
とかやった時になにが起きるのかはよくわからない。 C++ ではコンパイル時 にこのようなエラーは検出される。

16行目はダブルクォートの中の文字を画面に出す。'' cerr << '' と `` cout << '' はほとんど同じ意味だが、 UNIX等の OS でリダイレクトする時 に違いが出てくる。

C++ では Fortran と違って改行に意味はないので、文の終りはセミコロンで 示す。

17行目。 C++ では使いたい変数をプログラムのなかのどこでも好きなように 宣言できる。

18行目は入力。出力とは向きが逆。

19, 20 は計算で、この辺は Fortran と同じ。セミコロンがつくだけ。

21-23 はまた出力だが、最後に n をつけないと改行され ないということに注意。 Fortran だとなんでも改行されるが、 C++ だと「改行 しない」ということができる。

24は入力。

28 行目で配列の全要素を 0 にしている。 C++ の for は Fortran の DO loop と同様であるが、使い方がはるかに自由である。つまり、

for(式1; 式2; 式3) 文;

という形式で、やることは

  1. 式1を実行
  2. 式2を評価、成り立っていたら次に進む
  3. 文を評価(実際の計算)
  4. 式3を実行
  5. step 2 に戻る

というだけで、Fortran のように「制御変数」とかそういう特別な概念がある わけではない。例えば

for(;;);
というのも文法としては正しい(但し、これは無限ループである)。

なお、複数の文をまとめるには

{
  式1;
  式2;
}
というふうに中括弧で括る。

なお、「文」と「式」をだいぶいい加減に使っているが、この辺はわりとやや こしいのでちゃんと知りたいひとはC言語の文法書を見てみること。 C++ は初 心者向けの本というのがあまりよいものがないので、、、

配列要素は u[i]等。添字が 0 から始まるのが Fortran と違う。

次の関数 push_system は特に説明の必要はないであろう。

50行目。 if 文。 Fortran では .eq. とかそういうなんだかわからな いものを使うが、 C(++) でもやはり良くわからない ==, !=, <, >, >=, <= といったものを使う。 != が「等しくない」の他は意味は大体想像 できる通りであるが、注意して欲しいのは「等しい」が = ではなくて == で あること。

これは、「あらゆる式が値をもつ」という C 言語の特性によっている。 if 文のなかにかかれるのは整数値をもつ式ならなんでも良く、論理式でないとい けないとかそういう制約はない。さらに、 C言語では、代入式も値をもつ。 したがって、たとえば if (i=5) ... と書くと、これは i が 5 かどうかを調べるのではなく、まず i に 5 を代入して、しかるのちに その代入式の値を評価することになる。で、代入式の値はその代入された値と 同じということに決めてあるので、その値は 5 ということになる。

これは便利なときもある(例えば、 x = y = z = 0 とか書ける)が、 時々よくわからないミスを引き起こすのも確かである。

なお、サブルーチン(値を返さない関数)は、51行目のように call と かつけないでいきなり名前を書く。これは、これで「式」であって文法として は正しいということになる。

文法事項は大体これで終り。あと #define というのは、マクロ定義と いうものである。 #defineの後の語がプログラムにでてきたらその後のに置き換えられる。 これは C 言語の機能で、 C++ ではあまり使う必要はないが一応ここでは使っ てみた。

あ、その次の static というのにちょっと注意して欲しい。 Cは再帰呼 びだしが可能な言語であり、関数のローカル変数は呼ばれた時に生成される (自動変数)。

static をつけておくと、どこか固定されたアドレスにとられる。

main も関数であるが、どうせ一度しか呼ばれないのでその中の変数に は static とつけ ても意味は普通はないが、通常の UNIX 等の OS では自動変数に割り当てられて いるメモリが小さいので static をつけないで大きな配列を宣言すると実行時 にエラーになったりする。まあ、 static で大きい変数をとるのはあまり良い 方法ではないという意見もあるが、正しいプログラムを速く書くには安全な方 法である。



Jun Makino
Wed Oct 11 11:55:04 JST 2000