Fortranで数値計算｜科学技術計算の基本

Fortranで数値計算｜科学技術計算の基本

Fortranは、数値計算や科学技術計算を得意とするプログラミング言語です。
基本文法を覚えた後は、実際に数式や測定値をプログラムで処理してみると、
Fortranがどのような用途に向いているのか理解しやすくなります。

この記事では、Fortranを初めて科学技術計算に使う人向けに、
実数の精度、数学関数、配列を使った統計計算、行列計算、
数値積分、方程式の数値解法、簡単な物理計算までを順番に紹介します。

サンプルコード内には「!」を使ったコメント（注釈）を多めに入れています。
最初はコードをそのまま実行し、
数値を変更しながら結果がどう変わるか確認してみてください。

数値計算とは

数値計算とは、数式や物理法則などをコンピュータで数値として計算することです。

例えば、次のような処理が数値計算に含まれます。

• 平方根や三角関数を計算する
• 複数の測定値から平均値や標準偏差を求める
• 行列を計算する
• 積分を近似計算する
• 方程式の解を繰り返し計算で求める
• 物体の運動や温度変化を計算する

Fortranはもともと科学技術計算を強く意識して設計された言語であり、
このような処理を比較的自然に書くことができます。

実数の精度を意識する

Fortranで小数を扱う場合はrealを使いますが、
科学技術計算では「何桁程度の精度で計算するか」が重要になることがあります。

短い学習用プログラムでは通常のrealでも十分ですが、
より精度を確保したい場合はreal64を利用できます。

program precision_example
    use iso_fortran_env, only: real64
    ! iso_fortran_envはFortran標準の組み込みモジュールです。
    ! real64を利用すると、一般に倍精度相当の実数を扱えます。

    implicit none

    real(real64) :: x

    ! real64の実数定数であることを _real64 で明示します。
    x = 1.0_real64 / 3.0_real64

    print *, x

end program precision_example

use iso_fortran_env, only: real64は、
Fortran標準で用意されているreal64という種類を利用するための指定です。

科学技術計算の記事では、計算誤差の影響を減らすため、
以降の例でreal64を使う場合があります。

整数と実数の違いに注意する

数値計算で特に注意したいのが整数同士の割り算です。

integer :: a, b
real :: answer

a = 1
b = 2

answer = a / b

この場合、a / bは整数同士の割り算として処理されます。
小数を含む結果が必要な場合は、実数として計算します。

real :: a, b, answer

a = 1.0
b = 2.0

answer = a / b

科学技術計算では小数を扱うことが多いため、
変数の型を確認する習慣を付けておくとよいでしょう。

Fortranの数学関数

Fortranには、科学技術計算でよく使用する数学関数が標準で用意されています。

平方根を計算する

program sqrt_example
    implicit none

    real :: x
    real :: answer

    x = 25.0

    ! sqrt(x)で平方根を求めます。
    answer = sqrt(x)

    print *, "平方根 =", answer

end program sqrt_example

三角関数を計算する

sin、cos、tanへ渡す角度は、
通常ラジアンで指定します。
度数法の角度をそのまま入力しないよう注意してください。

例えば30度の正弦を求める場合、
30度をラジアンへ変換してからsinへ渡します。

program trig_example
    implicit none

    real :: degree
    real :: radian
    real :: answer

    real, parameter :: pi = 3.14159265

    degree = 30.0

    ! 度数法をラジアンへ変換します。
    radian = degree * pi / 180.0

    ! sinにはラジアンを渡します。
    answer = sin(radian)

    print *, "sin(30 degree) =", answer

end program trig_example

指数関数と対数を計算する

program exp_log_example
    implicit none

    real :: x

    x = 2.0

    ! eのx乗
    print *, "exp(x) =", exp(x)

    ! 自然対数
    print *, "log(x) =", log(x)

    ! 常用対数
    print *, "log10(x) =", log10(x)

end program exp_log_example

配列を使った数値計算

科学技術計算では、
1個の値だけではなく、
複数の測定値や計算値をまとめて扱うことが多くあります。

Fortranでは、配列全体に対して計算したり、
配列用の組み込み関数を利用したりできます。

測定値の合計と平均を求める

program average_example
    implicit none

    real :: data(5)
    real :: total
    real :: average

    ! 5個の測定値を配列へ保存します。
    data = [10.2, 10.5, 9.9, 10.4, 10.0]

    ! sum(data)で全要素の合計を求めます。
    total = sum(data)

    ! size(data)は配列の要素数です。
    ! real()を使い、整数を実数へ変換して割ります。
    average = total / real(size(data))

    print *, "合計 =", total
    print *, "平均 =", average

end program average_example

最大値と最小値を求める

program min_max_example
    implicit none

    real :: data(5)

    data = [10.2, 10.5, 9.9, 10.4, 10.0]

    ! maxvalは配列の最大値を求めます。
    print *, "最大値 =", maxval(data)

    ! minvalは配列の最小値を求めます。
    print *, "最小値 =", minval(data)

end program min_max_example

標準偏差を計算する

標準偏差は、データが平均値からどの程度ばらついているかを表す指標です。
ここでは学習用として、
母標準偏差に相当する形で計算します。

program standard_deviation
    implicit none

    real :: data(5)
    real :: average
    real :: variance
    real :: standard_deviation_value

    data = [10.2, 10.5, 9.9, 10.4, 10.0]

    ! 平均値を求めます。
    average = sum(data) / real(size(data))

    ! 配列全体について、
    ! 各値と平均値との差を2乗します。
    ! sum()で合計し、データ数で割ります。
    variance = sum((data - average) ** 2) / real(size(data))

    ! 分散の平方根が標準偏差です。
    standard_deviation_value = sqrt(variance)

    print *, "平均 =", average
    print *, "分散 =", variance
    print *, "標準偏差 =", standard_deviation_value

end program standard_deviation

data – averageのように、
配列と1つの値を直接計算できる点もFortranの特徴です。
配列の各要素からaverageが引かれます。

ベクトルの計算

ベクトルも1次元配列として扱えます。

program vector_example
    implicit none

    real :: a(3)
    real :: b(3)
    real :: c(3)

    a = [1.0, 2.0, 3.0]
    b = [4.0, 5.0, 6.0]

    ! 同じ位置の要素同士を足します。
    c = a + b

    print *, "a + b =", c

end program vector_example

内積を求める

Fortranにはdot_productという組み込み関数があります。
2つのベクトルの内積を求めることができます。

program dot_product_example
    implicit none

    real :: a(3)
    real :: b(3)
    real :: result

    a = [1.0, 2.0, 3.0]
    b = [4.0, 5.0, 6.0]

    ! aとbの内積を求めます。
    result = dot_product(a, b)

    print *, "内積 =", result

end program dot_product_example

行列を扱う

Fortranでは2次元配列を行列として扱うことができます。
例えば2行2列の行列は次のように宣言します。

real :: a(2, 2)

Fortranの配列は列優先で格納されますが、
初学者の段階では、
まず「a(行, 列)の形で値を指定できる」と理解しておけば十分です。

行列積を計算する

Fortranにはmatmulがあり、
行列積を直接計算できます。

program matrix_example
    implicit none

    real :: a(2, 2)
    real :: b(2, 2)
    real :: c(2, 2)

    ! reshapeを使って2×2配列を作ります。
    ! Fortranでは配列要素が列方向を先に埋める点に注意します。
    a = reshape([1.0, 3.0, &
                 2.0, 4.0], shape(a))

    b = reshape([5.0, 7.0, &
                 6.0, 8.0], shape(b))

    ! 行列積 A×B を計算します。
    c = matmul(a, b)

    print *, "C ="
    print *, c

end program matrix_example

行列の表示はそのままだと読みづらいことがあります。
本格的な行列計算では、
DO文を使って行ごとに整形して表示する方法もあります。

転置行列を求める

transposeを使うと、
行と列を入れ替えた転置行列を求められます。

real :: a(2, 3)
real :: b(3, 2)

b = transpose(a)

数値積分とは

積分を解析的に解けない場合や、
測定データから面積を求めたい場合には、
数値的に近似して積分値を求めることがあります。

最も基本的な方法の一つが台形公式です。
区間を小さく分割し、
各区間を台形として面積を足し合わせます。

台形公式で積分する

ここでは、0から1までのf(x)=x²を数値積分します。
正確な積分値は1/3ですが、
数値計算でどの程度近づくか確認できます。

program trapezoidal_rule
    use iso_fortran_env, only: real64
    implicit none

    integer :: i
    integer, parameter :: n = 1000

    real(real64) :: a
    real(real64) :: b
    real(real64) :: h
    real(real64) :: x
    real(real64) :: integral

    a = 0.0_real64
    b = 1.0_real64

    ! 区間幅 h を求めます。
    h = (b - a) / real(n, real64)

    ! 最初と最後の点は1/2の重みで加えます。
    integral = 0.5_real64 * (f(a) + f(b))

    ! 内側の点を順番に加えます。
    do i = 1, n - 1

        x = a + real(i, real64) * h
        integral = integral + f(x)

    end do

    ! 最後に区間幅を掛けます。
    integral = integral * h

    print *, "数値積分 =", integral
    print *, "1/3 =", 1.0_real64 / 3.0_real64

contains

    function f(x) result(y)
        use iso_fortran_env, only: real64
        implicit none

        real(real64), intent(in) :: x
        real(real64) :: y

        ! 積分する関数 f(x)=x^2
        y = x ** 2

    end function f

end program trapezoidal_rule

nを大きくすると区間を細かく分割するため、
一般には近似精度が高くなります。
ただし、計算量や丸め誤差もあるため、
単純に無限に大きくすればよいというわけではありません。

方程式を数値的に解く

方程式の解を公式で求められない場合、
繰り返し計算によって近似解を求める方法があります。

代表的な方法の一つがニュートン法です。
ここでは、
x² – 2 = 0の正の解、
つまり√2を求めます。

ニュートン法で√2を求める

program newton_method
    use iso_fortran_env, only: real64
    implicit none

    integer :: i
    integer, parameter :: max_iteration = 20

    real(real64) :: x
    real(real64) :: next_x
    real(real64), parameter :: tolerance = 1.0e-12_real64

    ! 最初の予想値です。
    x = 1.0_real64

    do i = 1, max_iteration

        ! ニュートン法
        ! x_new = x - f(x) / f'(x)
        ! f(x)=x^2-2、f'(x)=2xです。
        next_x = x - (x ** 2 - 2.0_real64) / &
                       (2.0_real64 * x)

        ! 前回との差が十分小さくなったら終了します。
        if (abs(next_x - x) < tolerance) then
            x = next_x
            exit
        end if

        x = next_x

    end do

    print *, "sqrt(2)の近似値 =", x
    print *, "sqrt(2) =", sqrt(2.0_real64)

end program newton_method

toleranceは、
どの程度値が変化しなくなったら
「解に十分近づいた」と判断するかを決める値です。

繰り返し計算では終了条件が重要

数値計算では、答えに近づくまで繰り返す処理がよく使われます。
その場合、終了条件を必ず考える必要があります。

例えば次の2つを用意すると安全です。

• 誤差が一定値より小さくなったら終了する
• 最大繰り返し回数に達したら終了する

最大回数を設定しないと、
条件によってはプログラムがいつまでも終了しない可能性があります。

簡単な物理計算｜自由落下

Fortranは物理計算とも相性のよい言語です。
ここでは空気抵抗を無視し、
静止状態から物体を落としたときの落下距離を計算します。

落下距離は、
時間t、重力加速度gを使って
1/2 × g × t²で求められます。

program free_fall
    implicit none

    real :: time
    real :: distance

    real, parameter :: gravity = 9.80665

    print *, "落下時間[s]を入力してください。"
    read *, time

    ! 静止状態からの自由落下距離を計算します。
    distance = 0.5 * gravity * time ** 2

    print *, "落下距離[m] =", distance

end program free_fall

このプログラムは非常に単純ですが、
時間を少しずつ変化させて計算すれば、
物体の運動を時系列で調べることもできます。

時間ごとの位置を計算する

program free_fall_table
    implicit none

    integer :: i

    real :: time
    real :: distance

    real, parameter :: gravity = 9.80665
    real, parameter :: dt = 0.5

    ! 0秒から5秒まで0.5秒間隔で計算します。
    do i = 0, 10

        time = real(i) * dt

        distance = 0.5 * gravity * time ** 2

        print *, time, distance

    end do

end program free_fall_table

このように、
一定時間ごとに物理量を計算して出力する方法は、
さまざまなシミュレーションの基本になります。

指数減衰を計算する

放射性物質の減衰や一部の物理・化学現象では、
指数関数的な変化が現れます。

program exponential_decay
    implicit none

    integer :: i

    real :: time
    real :: amount

    real, parameter :: initial_amount = 100.0
    real, parameter :: decay_constant = 0.2
    real, parameter :: dt = 1.0

    do i = 0, 10

        time = real(i) * dt

        ! N = N0 * exp(-λt)
        amount = initial_amount * &
                 exp(-decay_constant * time)

        print *, time, amount

    end do

end program exponential_decay

数値計算結果をファイルへ保存する

科学技術計算では、
大量の結果を画面へ表示するより、
ファイルへ保存した方が扱いやすい場合があります。

次の例では、
自由落下の時間と距離をfree_fall.txtへ保存します。

program save_free_fall
    implicit none

    integer :: i
    integer :: unit_number

    real :: time
    real :: distance

    real, parameter :: gravity = 9.80665
    real, parameter :: dt = 0.1

    ! 出力ファイルを作成します。
    open(newunit=unit_number, file="free_fall.txt", &
         status="replace", action="write")

    ! 見出しを書き込みます。
    write(unit_number, *) "time distance"

    do i = 0, 100

        time = real(i) * dt

        distance = 0.5 * gravity * time ** 2

        ! 時間と距離を1行ずつ保存します。
        write(unit_number, *) time, distance

    end do

    close(unit_number)

    print *, "free_fall.txtへ保存しました。"

end program save_free_fall

計算誤差とは

コンピュータは、すべての実数を完全に正確な値として保存できるわけではありません。
そのため、計算結果には小さな誤差が含まれることがあります。

例えば、0.1のような十進小数でも、
コンピュータ内部の2進数では有限桁で正確に表せない場合があります。

科学技術計算では、
次のような誤差を意識する必要があります。

• 浮動小数点による丸め誤差
• 数値積分などの近似誤差
• 繰り返し計算による誤差の蓄積
• 入力データそのものが持つ測定誤差

初学者の段階では、
「コンピュータの計算結果は必ず数学的に完全な値になるわけではない」
ということを覚えておけばよいでしょう。

実数同士を==で直接比較するときの注意

浮動小数点数には丸め誤差があるため、
計算結果が理論上同じ値でも、
内部ではわずかに異なる場合があります。

そのため、次のように完全一致だけで判定するよりも、
差の絶対値が十分小さいかを確認する方法が使われます。

real :: a
real :: b
real :: tolerance

a = 0.1 + 0.2
b = 0.3

tolerance = 1.0e-6

! aとbの差が十分小さいかを確認します。
if (abs(a - b) < tolerance) then
    print *, "ほぼ等しい"
end if

オーバーフローや不正な計算にも注意する

非常に大きな値を計算したり、
数学的に定義できない値を計算したりすると、
正常な数値が得られない場合があります。

例えば、次のようなケースです。

• 0で割る
• 実数として負の値の平方根を求める
• 非常に大きな指数関数を計算する
• log(0)や負の実数の対数を求める

入力値を確認するIF文を入れるなど、
計算前に条件を確認すると安全です。

数値計算で関数を使うメリット

数値積分や方程式の解法では、
計算対象となる数式を関数として分けておくと便利です。

function f(x) result(y)
    implicit none

    real, intent(in) :: x
    real :: y

    y = x ** 2 + 2.0 * x + 1.0

end function f

後から計算する式を変更したい場合も、
関数の中身を変更するだけで済みます。

Fortranの関数やサブルーチンは、
単にコードを整理するためだけでなく、
数値計算プログラムを再利用しやすくするためにも重要です。

科学技術計算でよく使う組み込み関数

ここまででできること

この記事の内容を理解すると、
Fortranだけでも次のような基礎的な科学技術計算ができます。

• 平方根・三角関数・指数・対数の計算
• 測定値の平均・最大・最小・標準偏差の計算
• ベクトルの加算と内積
• 行列積や転置行列の計算
• 台形公式による数値積分
• ニュートン法による方程式の近似解
• 自由落下などの簡単な物理計算
• 指数減衰の時系列計算
• 計算結果のファイル保存

さらに発展すると何ができる？

Fortranの数値計算をさらに発展させると、
次のようなテーマへ進むことができます。

• 連立一次方程式
• 固有値・固有ベクトル
• 微分方程式
• Runge-Kutta法
• モンテカルロ法
• 偏微分方程式
• 熱伝導計算
• 振動・波動シミュレーション
• 流体力学
• 有限要素法

より高度な行列計算では、
BLASやLAPACKなどの外部ライブラリを利用することもあります。
外部ライブラリについては、別の記事で導入方法と使い方を扱うことができます。

まとめ

Fortranは、
単純な四則演算だけでなく、
配列、数学関数、行列演算など、
科学技術計算に役立つ機能を標準で備えています。

特に、
配列全体をまとめて計算できることや、
dot_product、matmul、sumなどの組み込み関数を利用できることは、
数値計算を行う上で便利です。

また、台形公式やニュートン法のような数値計算法も、
DO文、IF文、関数などの基本文法を組み合わせることで実装できます。

科学技術計算では、
計算式そのものだけでなく、
実数の精度、近似誤差、終了条件、異常値などにも注意する必要があります。

基本的な数値計算に慣れたら、
次はFortran 77・Fortran 90/95・Modern Fortranのコードを比較し、
昔のFortranと現在のFortranで何が変わったのか確認してみると、
Fortranという言語の発展も理解しやすくなります。