ここでは、この先使用するルンゲ・クッタ法が記述されているサブルーチンの使い方を説明します。まず、適当なディレクトリ（例えば runge-kutta）をつくり、ここからRunge_Kutta.zipをダウンロードし、解凍してrunge.f90,
xrunge.f90 を入手して下さい。
1. コンパイル方法
このプログラムは二つのファイルから構成されています。ですから、これまでのように単体だけコンパイルするだけではダメです。次のようにコンパイルして下さい。
% gfortran -o rk4.exe runge.f90 xrunge.f90% ls% rk4.exe runge.f90 xrunge.f90 |
これは二つの Fortran90 のファイルをいっぺんにコンパイルしているわけです。このようにすると実行ファイル rk4.exe が作られます。
2. 基本的な使い方
基本的な使い方ですが、runge.f90 の方は手を加える必要はありません。手を加えてしまうと思わぬ動作をする可能性がありますので、こちらは手を触れないようにしてください。
いくつかの重要な変数、パラメータを解説します。
neq について、方程式の個数というのは、連立させる方程式がいくつあるか、ということです。このサンプルで解いているのは単振動の式
dv/dt = - (k/m) x
dx/dt = v
v: 質点の速度 x: 質点の位置
ですが、これは速度と位置の連立微分方程式になっているわけですから、式の数は2つ、となるわけです(プログラム中では k/m=1 としている)。もっと複雑な運動になれば方程式の数が増えることもありますので、その時は neq を式の数に合わせてやれば良いのです。
nmax について、これは何回計算をするかということを定めたパラメータです。これを大きくすればたくさん計算をすることになり、少なくすれば計算回数が減るわけです。これは見たい物理現象がどの程度の時間スケールを持っているかによって変えます。（もちろん、時間刻み幅との兼ね合いもありますので、それは計算毎に考えなくてはいけません）
x0 について、これは最初の時間を示しています。サンプル中では最初の時間は 0 となっていますが、これは変えなくてもいいのではないかと思います。
dx について、時間の刻み幅ですが、ルンゲ・クッタのルーチンの中で刻み幅は自動的に制御されていますので、これはこのままにしておいて下さい。
さて、プログラム中に y0(neq),yn(neq) という変数が出てくると思います。これは、y0 のほうが更新される前の値、つまり前ステップでの値、yn の方が更新された値、つまり今のステップで得られた値です。また、このサンプルプログラム中では y?(1) は速度 v, y?(2)
は位置 x を指しています。
微分方程式を記述している部分は、 function ex2(i,x,y) というところです。これは i の値によって方程式を指定しています。サンプル中では i=1 の時が dv/dt=-x
を表しており、i=2 のときが dx/dt=v を表しています。この ex2 を書き換えれば、基本的にはどのような常微分方程式でも（数値的に）解けることになります。
具体的にどのように書き換えられるかについては、演習時間に、TA が作ったいくつかのサンプルを示しますので、今回与えられたサンプルと見比べて、どのように書き換えれば良いかということを考えて下さい。