今回は 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) 文;
という形式で、やることは
というだけで、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 で大きい変数をとるのはあまり良い 方法ではないという意見もあるが、正しいプログラムを速く書くには安全な方 法である。