物理を題材にしたFortran 入門-1     竹下徹

これから始めるFortran の練習コースはパリティーという雑誌に中村純さんがお書きになった「物理屋のためのFORTRAN入門」という文章のなかの教材として取り上げる物理現象 を一部お借りしています。

物理を学ぶ人のためのFortran入門webです。中村さんのおっしゃっているように動機付けなしに言語は覚えられない。少なくとも物理の学生に は物理現象がその動機付けとなるであろう。従って扱う対象を物理現象(できれば大学生が興味を持ちそうな、あるいはおもしろいと感じるような)に限定して います。

ではFortran プログラムの勉強です.さて初めは皆さんがプログラムについて何も知らないと仮定して始めます.ただしeditorは知っていて,プログラムを作るあるい は変更することはできるとして始めます.最低限のルールから始めましょう.つぎの例を見て下さい.

---------------------------------------------プログラム例1

program meher

implicit none

integer me,res

real her

me = 2

her=2.3

res=me*her

write(*,*) " result of you and her is = ",res

stop

end

--------------------------------------------プログラム例1の説明

program meher  :プログラム名がmeher であることをまず言う必要があります。

implicit none  :  おまじないの一つですが、ここで使う変数は全て自分で整数型か実数型か定義するぞという意味です。

integer me,res   :右側二つの変数 その名を me とher (両者の間だはコンマ, でつなぐ) が整数(integer)で有ることを宣言します。

real her      :同様に 変数名 her(一個だけなのでコンマなし)が実数型の変数であることを宣言する。

me = 2       :ここで整数型であることを決めた変数 me (左辺) へ 2という数値をいれます。入れるときは必ず、右辺の「物」を左辺に入れます。これは実は計算器のなかでmeと名付けられた場所(メモリー上 のアドレスのこと)に2という整数値を入れておきます。以後のプログラム中で変数meが使われると、この中に入っている数値が使われます。また例えば違う 変数名 integer x me = x とした場合、x  という名前のついた箱の中身の数値を取り出してこれと同じ数値を me という筺にいれる(従って x には元の数値が残っています)。

her=2.3    :実数型の変数の箱で名前がherという箱のなかに 2.3 という実数値を入れます。

res=me*her   : me という変数名の箱から取り出した変数の値とherという変数名の箱から取り出した変数の値のかけ算(*) を行い、その結果を整数型の箱 resのなかにいれます。ここで妙なことが起こるのです。そのためにはつぎの整数と実数のメモリー上での表し方を知らねば なりません。

write(*,*) " result of you and her is = ",res  :write 文で 変数res を印刷します。 write(*,*) の左側の*が画面への出力をしめし、右側の*がその出力の仕方がデフォルトに従うことをしめします。また" result of you and her is = " はそのまま印刷されます。

stop      :実行はこれでおしまいであることを示します。

end       :プログラムはこれでおしまいであることを示します。

------------------------------------------プログラム例1の説明 おしまい------

ただし、全ての行の先頭には6個以上の空白つまりスペースがはいっています。このファイルをmeher.fとして保存したとしましょう.最後の .fは必ずつけて下さい。Fortran-コンパイラーにこれがFortranのプログラムであることを教えています。

Fortran-コンパイラーを通して実行ファイルmeherを以下のようにして作ります.Fortranのプログラムは人間が皆さんが理解できる 形(Fortran言語)で書かれていますが、計算機は理解できません。そこで、計算機にわかるようにプログラムを変換してやるにのがFortran-コ ンパイラーの役目です。

>はUNIXが画面上に入力待ち(prompt)であることを示し、みなさんの入力の必要はありません。

>g77 meher.f -o meher

Fortran-コンパイラーはg77という実行形式のプログラムです(通常のLinuxにはこういう生前で入って居ますし、pathが通じている はずです、もしcommand not found と言われてしまったら、locate g77とやってどこにあるか調べてください)。g77 meher.f でg77とmeher.fの間にはスペースが一つ入っています。さもないと、g77meher.fというつながったプログラムを実行しようとしている事に なります。-oの前と後も同じです。-oはoutput オプションの意味で、出力ファイル名をmeherにしなさいと言っています。結果はmeherを実行してつぎの様に得られます.

>./meher

result of you and her is = 4

なぜ4になったのでしょうか?

ちなみに >g77 meher.f とやると、実行可能な出力はa.outとなります。 >./a.outとして実行してみてください。同じ結果を得るはずです。

それはFortranの計算では整数型の変数(ここではmeとres)と実数型の変数(ここではher)があり両者はハッキリ区別されていることに 注意してください.上のプログラム例1では、 integer me,res  としてmeとresという2つの変数を整数型だぞとしています。同様に real her では herという名の変数を実数型だぞとしています。 なぜこのような型を定義しなければならないかは、この章の最期に述べます。そこでたとえば,実数型x実 数型=実数型になりますし,実数型x整数型=実数型,整数型x整数型=整数型という演算をします.このプログラムでは,res=me*herという整数型 x実数型の結果は実数型となるところをresという整数型変数に代入して,2x2.3=4.6 をまるめて4となっています.普通に我々が期待する結果と異なる理由を理解できましたか?

では正しい結果を得るように,皆さんはこのプログラムを改造してください.

ここで少しFortranの規則についてふれます.プログラムは program から end までです.つまりこの間が本当の実行文が書いて有る領域であることがわかります。また全ての実行文は1行に1つだけとしています。整数型変数であること は,integer 変数名; ここでは (integer me,res)で、meとres という整数型の変数を定義しています。 同様に,real 変数名; (real her )で変数herが実数型であることが分かります.Fortranでは登場するすべての変数をこのように宣言する必要はありませんが、こう書くべきと思いま す.write(*,*) は画面に以下のことを表示させる関数です。" result of you and her is = ",res で、double-quataion"の間の形式で画面に出力します。ですから画面に result of you and her is = と表示されます。変数はresを表示しています。

さてつぎの例を見て下さい.

---------------------------------------------プログラム例2

program sum

implicit none

integer k,sum,kmax

kmax=10

sum = 0

do k=0,kmax

sum=sum+k*k

write(*,*) "  k and sum ",k,sum

enddo

stop

end

---------------------------------------------

ここでは整数型変数sumがkの2乗の和になっている事が読みとれますか?.ちなみにfortranではkの2乗を k*kとかくこともできますし、 k**2 とかくこともできます・ ただし,初期値 sumに0が代入されています.これを抜かして走らせてみて下さ い。ひょっとするととんでもない結果を得て、「計算機がくるってる!」と叫ぶかもしれませんが、それは皆さんのミスです。計算機は自動的に変数の定義をす ればゼロをいれてくれるとはかぎりません。ここでループ(繰り替えし)はdoからenddoの間で実行されます.まずk=0の場合の計算が実行されます。 sumはゼロのままです。次に k=k+1 が実行され k=1 でsumが1になります.このとき,kの最大値は10に等しくなるまで計算が実行されま す.つまり次の式を計算しているわけです。do文の kはかならず整数型の変数にしてくだ さい。さもないと終わらないかもしれません。こうして整数の2乗和を計算できます.皆さんは初期値が0でない場合や,終了が10でない場合,あるいはもっ と進んで1/k 等の実数の和をつくプログラムを書くことができます.

ですがちょっと注意をしておきます。1/k は実数ですが、kを整数で切る、つまり program中で integer k としている限り、 1/k も整数扱いを計算機はします。それは上で学んだとおりです。ですから、和を取るとき、つまり、 sum=sum+1/k では、答えはずううううとゼロです。やってみてください。 これを逃げるには、   1.0/k と書いてみてください。 1.0 /k は計算機が実数と思ってくれるンですよ。 

 では k! を計算する プログラムを作れと言われたらどうしましょう。ヒントは今までは和を取ってきたが、今度は積をとるのだ!です。

----------------------------------------

program factorial

implicit none

  integer ifact,nmax,l,k
ifact=1
write(*,*) 'input max loop '
read(*,*) nmax
do k=1,nmax
ifact=ifact*k
write(*,"(i5,i15)") k,ifact
  enndo
stop
end

--------------------------------------------

コンパイルして実行してみましょう: > g77 factorial.f -o factrial > factorial

 最大値に15をいれてみると、あれれれ、なんか変です

input max loop
15
2 2
3 6
4 24
5 120
6 720
7 5040
8 40320
9 362880
10 3628800
11 39916800
12 479001600
13 1932053504
14 1278945280
15 2004310016

明らかに 13! で間違っているのが判りますか?計算機が間違いをするはずがないと思っていませんか? プログラムは正しく書いたし、実行時の  k が小さい値では結果も正しかったのですから、途中から間違いが起こったのは、計算機の限界を超えてしまったのです。計算機のこのプログラムでは integer ifact としたことで実はifactという変数の場所をメモリー上にとるとき有る広さを約束に応じて(約束した覚えはないかもしれませんが:しちゃっ たんです)確保しています。それはここでは4バイト、つまり32ビットです。32ビットで表すことのできる最大の整数値は2の31乗=2**31です、 32乗じゃあないのかって?それは1ビットは符号に使っているので、一つビットが減って2の31乗 2*831=2147483647おおよそ2.147x10**10、文字で10乗を **10のように書きます。10桁しかないのです。13!= 6227020800 なのですが、この値は2**31を越えています。そうするともう計算機が何を計算して返してくるのかは、皆目見当がつきません。実 際そうなりました。 ここで言いたいことは、計算機の能力の限界が有ることを知り使わないといけない。計算機は正しい計算をするというのは、正確ではな く、プログラマー(あなた)が正しい動かせ方をしたら、正しい計算をする程度のものでしょう。さてもう一つおまけで、先ほどの k**2 の和を取るプログラムで kmaxをやたら、大きくしてください。でも2**31以下にしてくださいね。それでも間違った結果が現れるのをみるはずです。

プログラムでは write(*,"(i5,i15)") k,ifact というのが使ってあります。これはformat出力といい、画面に出た数値をきれいに並べて見栄えを良くし、内容の把握をしやすくするた めの仕掛けです。右そろいになっているでしょう。 write(*,*) k,ifact では左そろいになちゃうのです。ちなみにi5とかi15は i が整数型の変数を印刷する事を示し、5や 15はそのからむ数(何桁分印刷するか)に対応しています。つまり i5 では 整数型 5桁の数値を左詰めで員サルしてくれます。

さてここでの課題は、 の計算をせよです、この計算はexp(x)であることはTalyer 展開から知っていることでしょう。

これを用いてつぎの問題を考えてみましょう.ある年のある場所の例えば虫の固体数Njとして,つぎの年の固体数Nj+1は最も簡単な仮定として,前 年の固体数に比例する.比例定数をcとしてすなわち,  Nj+1= cNj,

しかし固体数のモデルとしては明らかにこれでは荒すぎます.環境の為に制限が存在し(例えば餌がなくなる),固体数はあまり増えなくなる.すなわり 限界Nmaxに近づくと増加は抑制されるとしたほうが,現実的でしょう.そこで,このモデルを,

Nj+1=cNj(Nmax-Nj) と書き換えます.この両辺をNmaxで割って,xj=Nj/Nmaxとすると,xj+1=axj(1-xj), a=cNmax,

として,つぎのプログラムを見て下さい.

--------------------------------------------プログラム例3

program chaos

implicit none

integer n,nmax

real x,a

nmax= 100

x=0.1

write( *,*) " input initial value a : "

read (*,*) a

write( *,*) " inputs are 繰り返し、a xの初期値 ",nmax,a,x

write( *,*) " n: x \n"

do n=0, nmax

x= a*x*(1.0-x)

write(*,*) n,x

enddo

stop

end

------------------------------------------------

ここで,readは,キーボードからのプログラムへのデータ入力を行います.この場合aの値を入力させています.

この結果がaに依存して収束の様子が違ってくることを確かめて下さい.つまりいろいろのaの場合を試し、年毎の個体数の変動が収束する場合と、振動 するばあいが有ることを確かめて下さい。

aがおおよそ3より大きくなると,xは収束せず振動します.我々の数理モデル(漸化式)Nj+1=cNj(Nmax-Nj)は相当簡単なものだった のですが、結果は非常に複雑なふるまいをすることが見て取れますね。自然界で複雑な振る舞いを発見しても、それは必ずしも複雑な機構を反映するとは限りま せん、という訳でしょうか?

この回の宿題

プログラム例2でS=(1/1+1/2+1/3+1/4+...........+1/x)を計算するプログラムを作りx=23の結果を計算する。

次にプログラム例3で収束の境目となるaの値を可能な限りの精度で求める。

====================付録=====================================

この付録では、計算機の「計算」で要注意事項を学んでもらいます。

その1:計算する変数や定数の型について

   計算機で言語をもちいて計算するとき、まず変数の名前をきめて登録します。このとき、同時に変数の型を決めます。通常は整数型か実数型かで す。なぜならこれらを表す内部のビットの並びが整数と実数では全く異なるからです。さらに整数、実数にも長さ(長いほど精度がよい)を決めます。通常の PC あるいはコンピュータは32bitコンピュータなので、32ビットbit で表される長さになります。32ビットで整数はどのように表されるかを説明しましょう。

計算機の内部表現の最小単位はビットです。ビットは1あるいは0です。これをON/OFFとよんだり、True/Fauls あるいは真/偽 とよんだりしますが、要はどちらかしか取らない物です。これは数学では2進数です。1ビットが表すことのできる組み合わせは2つ(0か 1)です。正確には(0(10進数)か1(10進数))です、と書くべきです。なぜなら私たちの通常の数の呼び方は10進法に基づいており、それをわざわ ざいいません。なにも言わねば10進法であることが当たり前だからです。しかし、ここではわざわざ書きます、なぜなら2進法を議論するからです。

さて10進法では0,1,2,3,4,。。。。9の10個の数値を使います。だから10進法なのですが、2進法では0,1の2種類です。これは1桁 を表すことが出来ます。2桁目の登場により、表す事のできる組み合わせの数は飛躍的に増加します。10進法では0,1,、、9ときてつぎ9より1多い数値 を10として桁を一つ増やします。(当たり前ですね)。更に1多い数値は11さらに12,13,、と数を増やして行けます。19まで来ると、2桁目を1あ げて20としますね。つまり二桁の任意の数が ij (i,j=0,1,2,...9) のとき、ij=i*10^1+j*10^0 と書くことができます。ここで^はべき乗を表します。これではわかりにくいかもしれませんので、通常の式にしましょう。と上の式を読んでください。でもこれでは一般 化しにくいのでi,jを書き換えますよ。こ れが10進数2桁の表記です。各桁の数値は整数でN1,N2と書かれています。あとはn+1桁に一般化します。

2進数の値も同様にやってみましょう。1桁は0,1の2個で表すことができます。と右辺を10進数で書くことができます。0, 1と1づつ加えてきてそのつぎは10進数の2は2進数では桁を1つ増やして対応します。つまり10(2進数)=10(2)と書きます。2進数10(2) は10進数(右辺)でと かけますし、3(10)はと 表します。桁数を増やすと、各桁の値はビットと呼ばれます。
10進数 2進数 ビット3 ビット2 ビット1 ビット0

0

0

0

0

0

0

1

1

0

0

0

1

2

10

0

0

1

0

3

11

0

0

1

1

4

100

0

1

0

0

5

101

0

1

0

1

6

110

0

1

1

0

7

111

0

1

1

1

8

1000

1

0

0

0

9

1001

1

0

0

1

10

1010

1

0

1

0

11

1011

1

0

1

1

12

1100

1

1

0

0

13

1101

1

1

0

1

14

1110

1

1

1

0

15

1111

1

1

1

1
2進数4桁(4ビット)で表すことのできる組み合わせは2の4乗=16通りあり、10進数とすると、0から15までとなります。では32ビットで整数はい くつまでと行きたいところですが、世の中には整数でも正の数と負の数がありますので、32ビットのうちで一番左側を符号ビットとして使いますので、実は数 を表すのは31ビットつまり、0から2の31乗までとマイナス2の31乗までの結局(2の32乗ー1)個ということになります。最大値である2の31乗は 10進数では2147483647という10桁の10進数です。たった10桁(10進数)しか表せないのです。2進数を10進数に変換する事は容易です、 なぜなら各 ビットの0あるいは1をnに代入して右辺を計算すればよりからです。例えば2進11001(2)は

この方法でと なり10進25であることが解ります。その反対に10進数を与えられて2進数をだすには上の式の各項の係数を求めることになりますが、実は簡単な方法があ ります。それは与えられた10進数を2で割ってあまりをそれぞれの桁で記憶しておきます。例えば、25(10)の場合最後に商が0になる前にやめてこれを左から 上へ向かって読むと11001(2)=25(10)をえることが出来ます。この方法は上の展開式の各項を計算しているだけですね。

一方実数はさらに注意が必要です。実数とは少数点以下のの数値のことですので、まずは 小数点以下のみでできている10進の数値を考えます。上で やったようにと 拡張できます。小数点下n桁も同様ですね。さてこれを2進数にする事を考えましょう。こんどは2で割ってゆく訳に行きません。余りが整数にならないからで す。どうしてかは置いといて割る代わりに掛けます。理由は考えてくださいね。さて0.25に2をどんどん掛けてゆきます。これは簡単な例ですが、点線で囲んで上から読 んだ01が小数点以下の表記になります。つまり 0.25(10)=0.01(2)という訳です。小数点以下が全部0になったらこれ以上の計算はありませ ん。0.01(2)を上の2進数にする方法でであることが確かめられます。他の例で例えば 0.56を試してみてください。いつまでたっても最後0になりません。つまり、10進数0.56を2進数で表すことが出来ないのです。これは計算機の重大 問題です。計算が正しくないのです!(といっても決定的に違うのではなくほんのちょっぴり違うという微妙なでも注意の必要な事態となります)この方法によ り、例えば25.25(10)=11001.01(2)と書くことできます。ですが、計算機で実数を扱う上でこの方法(固定小数点法)はあまり効率がよく ありません。小数点いくらまで桁を使うのか、また上の桁は31ビット使っても10進でたったの10桁という問題をクリアしねけれなりません。そこで多くの 計算機ではつぎの方法が取られています。これを浮動小数点法とよびIEEE754(ここでは正規表現だけを使います)という統一規格が存在してつぎのよう に定めています。小数点以下の数を仮数と呼びます。これを23ビット使って表します。また小数点より上の整数部を指数部とよび8ビットを使います。さらに 1ビット先頭に符号ビットを付けます。32ビットを左から[符号ビット(1)][指数部(8)][仮数部(23)]=[s][e][a]のようにな内部表 現になっています。ただし、指数部には注意が必要です。それは指数部=0000 0000 (2) =0(10) は2の128乗に対応することです。つまり仮数部のビット並びをaとし、指数部をe、全体符号をsとすると、これは実数 (-1)^s×[1.a](2)×2^{[e](10)−127} を表わすということです。例を挙げましょう、65.25(10)は32ビット表現でどうなるでしょうか。正の数ですからs=0。65.25(10)= 1000001.01(2)と固定小数点で表されます。これを正規化して65.25(10)=1.00000101(2)x2^6と書けますから、e (10)-127=7よってe(10)=133(10)=10000101(2)これが指数部です。最後に00000101(2)はそのまま仮数部なので 結果32ビットは 0 10000101 00000101000000000000000 と表されます。コンピュータでこのような0000並びは書く桁が多すぎますので、16進がよく使われます。16進数は  0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F の16個が使われます。16進数一桁は2進数4桁の整数を表現できます。例えば1001(2)=9(16), 1010(2)=A(16),1011(2)=B(16),1011(2)=C(16),1100(2)=D(16),1101(2)=E(16),1111 (2)=F(16)となります。これを使うと

65.25(10)= 0 10000101 00000101000000000000000(2)=42828000(16)と書かれます。

でももしこれを整数として読んでしまうとどうなるでしょうか?

だからプログラム言語ではあらかじめ変数の型をあなたが(人が)決めるのです。こいつは整数だって。

それは計算機は万能ではありません。特に「計算」に関して。

次のプログラムを見て下さい。

--------------------------

program realacc

implicit none

integer n

real x,z

real*8 y

x = 0.0

y = 0.0

z = 100000.0

do n=0,10000

x = x+z

y = x+z

write(*,*) " real sum x and double-real y ",n,x/z,y/z

enddo

stop

end

------------------------------

このプログラムでは、実数で単精度(real)のxと実数で倍精度(real*8)のyが、z=10万という単精度の実数を複数回足しあげられてい ます。最終的に印刷writeするときはこれをzで割っていますから、答えは足しあげた回数になるはず(実数を実数で割った実数としています)です が、、、、、、

実際これを実行すると、5000回ぐらいで妙な事が起こります。自分で確かめてください。理由は単精度の計算と倍精度の計算で桁落ちと呼ばれる現象 が起こっているためにこんな事が見えてしまいます。こんな簡単な例でもわかるように、一言で計算すると言っても、計算機の能力や特徴をよく理解したプログ ラムを組まないと、結果は信用できないものに成りかねません。また皆さんはよくわかった例で計算が正しいことを確かめたうえで、とても手計算ではできない ことを計算器にやらせる事がはじめてできることを理解しておいて下さい。

===================付録おしまい==========================