Fortran実践サンプルコード10選｜二次方程式・統計・数値積分・グラフ描画

Fortran実践サンプルコード10選｜二次方程式・統計・数値積分・グラフ描画

この記事では、Fortranの基本文法を一通り学んだ人向けに、
実際に動かしながら復習できる実践サンプルコードを10個紹介します。

単純な計算から始めて、
条件分岐、繰り返し、配列、関数、ファイル出力、
数値積分、反復計算、外部ライブラリによるグラフ描画まで、
少しずつ難しくなる順番で構成しています。

各コードには「!」による注釈を多めに入れています。
最初はそのまま実行し、その後で数値や式を変更して試してみてください。

サンプル1｜三角形と円の面積を求める

最初はREAD、変数、四則演算の復習です。
三角形の底辺・高さと、円の半径を入力して面積を求めます。

program area_example
    implicit none

    real :: base
    real :: height
    real :: radius

    real :: triangle_area
    real :: circle_area

    real, parameter :: pi = 3.14159265

    print *, "三角形の底辺を入力してください。"
    read *, base

    print *, "三角形の高さを入力してください。"
    read *, height

    print *, "円の半径を入力してください。"
    read *, radius

    ! 三角形の面積 = 底辺 × 高さ ÷ 2
    triangle_area = base * height / 2.0

    ! 円の面積 = πr^2
    circle_area = pi * radius ** 2

    print *, "三角形の面積 =", triangle_area
    print *, "円の面積 =", circle_area

end program area_example

サンプル2｜摂氏・華氏・ケルビンを変換する

温度変換を使って、実数計算と複数の計算結果を扱います。

program temperature_conversion
    implicit none

    real :: celsius
    real :: fahrenheit
    real :: kelvin

    print *, "摂氏温度を入力してください。"
    read *, celsius

    ! 摂氏 → 華氏
    fahrenheit = celsius * 9.0 / 5.0 + 32.0

    ! 摂氏 → 絶対温度
    kelvin = celsius + 273.15

    print *, "華氏 =", fahrenheit
    print *, "ケルビン =", kelvin

end program temperature_conversion

サンプル3｜密度を計算し、入力値も確認する

科学計算でよく使う
密度 = 質量 ÷ 体積
を計算します。

体積が0以下の場合には計算しないように、
IF文で入力チェックを行います。

program density_example
    implicit none

    real :: mass
    real :: volume
    real :: density

    print *, "質量を入力してください。"
    read *, mass

    print *, "体積を入力してください。"
    read *, volume

    if (volume > 0.0) then

        density = mass / volume

        print *, "密度 =", density

    else

        print *, "体積は0より大きい値を入力してください。"

    end if

end program density_example

サンプル4｜二次方程式を解く

二次方程式
ax² + bx + c = 0
を解きます。

判別式
D = b² – 4ac
を使って、
2つの実数解、重解、複素数解を判定します。

program quadratic_equation
    implicit none

    real :: a
    real :: b
    real :: c

    real :: d

    real :: x1
    real :: x2

    complex :: z1
    complex :: z2

    print *, "aを入力してください。"
    read *, a

    print *, "bを入力してください。"
    read *, b

    print *, "cを入力してください。"
    read *, c

    ! a=0では二次方程式ではありません。
    if (a == 0.0) then

        print *, "aは0以外を入力してください。"

    else

        ! 判別式
        d = b ** 2 - 4.0 * a * c

        if (d > 0.0) then

            ! 異なる2つの実数解
            x1 = (-b + sqrt(d)) / (2.0 * a)
            x2 = (-b - sqrt(d)) / (2.0 * a)

            print *, "解1 =", x1
            print *, "解2 =", x2

        else if (d == 0.0) then

            ! 重解
            x1 = -b / (2.0 * a)

            print *, "重解 =", x1

        else

            ! d<0ではsqrt(d)を実数として計算できないので、
            ! complex型を使って複素数解を求めます。

            z1 = cmplx(-b, sqrt(-d)) / (2.0 * a)
            z2 = cmplx(-b, -sqrt(-d)) / (2.0 * a)

            print *, "複素数解1 =", z1
            print *, "複素数解2 =", z2

        end if

    end if

end program quadratic_equation

このサンプルでは、
Fortranがcomplex型を標準で持っていることも確認できます。

サンプル5｜測定値の平均・最大・最小・標準偏差

複数の測定値を配列へ入れ、
平均、最大値、最小値、標準偏差を求めます。

program statistics_example
    implicit none

    real :: data(5)

    real :: average
    real :: variance
    real :: standard_deviation

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

    ! 平均
    average = sum(data) / real(size(data))

    ! 分散
    variance = sum((data - average) ** 2) / real(size(data))

    ! 標準偏差
    standard_deviation = sqrt(variance)

    print *, "平均 =", average
    print *, "最大値 =", maxval(data)
    print *, "最小値 =", minval(data)
    print *, "標準偏差 =", standard_deviation

end program statistics_example

sum、size、maxval、minvalのような
Fortran標準の配列関数をまとめて確認できます。

サンプル6｜素数を判定する

DO文、IF文、mod、exitを組み合わせた例です。

program prime_check
    implicit none

    integer :: number
    integer :: i

    logical :: is_prime

    print *, "2以上の整数を入力してください。"
    read *, number

    is_prime = .true.

    if (number < 2) then

        is_prime = .false.

    else

        ! 2からsqrt(number)まで調べれば十分です。
        do i = 2, int(sqrt(real(number)))

            if (mod(number, i) == 0) then

                ! 割り切れたので素数ではありません。
                is_prime = .false.

                ! これ以上調べる必要がないのでDOを終了します。
                exit

            end if

        end do

    end if

    if (is_prime) then
        print *, number, "は素数です。"
    else
        print *, number, "は素数ではありません。"
    end if

end program prime_check

サンプル7｜台形公式で数値積分する

ここから少し数値計算らしい内容になります。
0から1までの
f(x)=x²
を台形公式で積分します。

program trapezoidal_integration
    implicit none

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

    real :: a
    real :: b
    real :: h
    real :: x
    real :: integral

    a = 0.0
    b = 1.0

    ! 1区間の幅
    h = (b - a) / real(n)

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

    do i = 1, n - 1

        x = a + real(i) * h

        integral = integral + f(x)

    end do

    integral = integral * h

    print *, "数値積分 =", integral
    print *, "正確な値  =", 1.0 / 3.0

contains

    function f(x) result(y)
        implicit none

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

        ! 積分する関数
        y = x ** 2

    end function f

end program trapezoidal_integration

関数f(x)の中身を変更すれば、
別の関数の数値積分にも利用できます。

サンプル8｜ニュートン法で平方根2を求める

ニュートン法を使って、
x² - 2 = 0の解を求めます。

program newton_method
    implicit none

    integer :: i

    real :: x
    real :: next_x

    real, parameter :: tolerance = 1.0e-6
    integer, parameter :: max_iteration = 100

    ! 最初の予想値
    x = 1.0

    do i = 1, max_iteration

        ! f(x)=x^2-2
        ! f'(x)=2x
        !
        ! ニュートン法：
        ! x_new = x - f(x)/f'(x)

        next_x = x - (x ** 2 - 2.0) / (2.0 * 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)

end program newton_method

この例では、
反復計算・誤差判定・最大繰り返し回数
という数値計算の基本を学べます。

サンプル9｜自由落下を計算してCSVへ保存する

時間ごとの自由落下距離と速度を計算し、
free_fall.csvへ保存します。

このサンプルでは、
DO文、物理式、ファイル操作をまとめて使います。

program free_fall_csv
    implicit none

    integer :: i
    integer :: unit_number

    real :: time
    real :: distance
    real :: velocity

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

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

    ! 1行目に項目名を書きます。
    write(unit_number, '(A)') "time,distance,velocity"

    ! 0秒から10秒まで0.1秒刻みで計算します。
    do i = 0, 100

        time = real(i) * dt

        ! 静止状態からの自由落下
        distance = 0.5 * gravity * time ** 2

        ! 速度 v = gt
        velocity = gravity * time

        ! CSV形式で保存します。
        write(unit_number, '(F8.3,",",F12.5,",",F12.5)') &
            time, distance, velocity

    end do

    close(unit_number)

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

end program free_fall_csv

CSV形式にしておけば、
表計算ソフトなどでも開きやすくなります。

サンプル10｜PLplotを使ってy=x²のグラフを描く

最後は外部ライブラリを利用する例です。
標準Fortranにはグラフ描画機能はありませんが、
PLplotのような科学技術向け描画ライブラリを利用すると、
Fortranプログラムからグラフを描画できます。

ここでは、
y=x²のデータを作り、
PLplotで線グラフとして描画します。

program plot_example
    use plplot
    implicit none

    integer, parameter :: n = 101
    integer :: i

    ! PLplotで使用する実数型です。
    real(kind=plflt) :: x(n)
    real(kind=plflt) :: y(n)

    ! -5から5までの点を作ります。
    do i = 1, n

        x(i) = -5.0_plflt + &
               10.0_plflt * real(i - 1, kind=plflt) / &
               real(n - 1, kind=plflt)

        ! y = x^2
        y(i) = x(i) ** 2

    end do

    ! PNG出力を選択します。
    ! 利用可能なデバイス名は
    ! PLplotのビルド環境によって異なる場合があります。
    call plsdev("pngqt")

    ! 出力ファイル名を指定します。
    call plsfnam("quadratic_graph.png")

    ! PLplotを初期化します。
    call plinit()

    ! 新しいページを開始します。
    call pladv(0)

    ! グラフ表示領域を設定します。
    call plvpor(0.15_plflt, 0.85_plflt, &
                0.15_plflt, 0.85_plflt)

    ! x軸 -5～5、y軸 0～25
    call plwind(-5.0_plflt, 5.0_plflt, &
                 0.0_plflt, 25.0_plflt)

    ! 軸を描画します。
    call plbox("bcnst", 0.0_plflt, 0, &
               "bcnstv", 0.0_plflt, 0)

    ! 軸ラベルとタイトルを付けます。
    call pllab("x", "y", "y = x^2")

    ! x,yの点を線で結びます。
    call plline(x, y)

    ! 描画処理を終了します。
    call plend()

    print *, "quadratic_graph.pngを作成しました。"

end program plot_example

PLplotを使うときの注意

この10番だけは、
PLplotを別途インストールし、
コンパイル時にPLplotをリンクする必要があります。

具体的なコンパイルコマンドは、
PLplotのインストール方法やOSによって異なります。

例えば、pkg-configが利用できる環境では、
インストールされたPLplotの設定を確認しながら
コンパイルオプションを取得できる場合があります。

PLplotはFortranバインディングを提供しており、
標準的なXYプロットだけでなく、
等高線、3D表面、メッシュ、棒グラフなども扱えます。

グラフ描画サンプルで学べること

• 配列へx座標とy座標を保存する
• DO文で関数値を計算する
• 外部ライブラリのMODULEをuseする
• 外部ライブラリのサブルーチンをcallする
• 計算結果をグラフとして可視化する

10個のサンプルで学べる内容

どの順番で試すとよい？

Fortranを始めたばかりなら、
1番から順番に実行するのがおすすめです。

1～3では基本的な入出力と計算、
4～6では条件分岐・配列・繰り返し、
7～9では数値計算・物理計算・ファイル操作、
10では外部ライブラリへ進みます。

つまり、この10本を順番に試すことで、
これまで学んだFortranの主要機能を一通り復習できます。

サンプルコードを改造してみる

コードをそのまま動かした後は、
少しずつ変更してみると理解が深まります。

• 面積計算に台形や球の体積を追加する
• 温度変換をメニュー形式にする
• 二次方程式の係数をファイルから読み込む
• 統計計算のデータ数を増やす
• 素数を1個ではなく一定範囲で一覧表示する
• 数値積分の関数をsin(x)へ変更する
• ニュートン法で別の方程式を解く
• 自由落下に初速度を追加する
• CSVへ加速度も保存する
• PLplotでsin(x)や自由落下データを描画する

まとめ

この実践サンプル集では、
Fortranの基本文法から数値計算、
ファイル操作、外部ライブラリまでを
10個のプログラムでまとめて確認しました。

特に二次方程式ではIF文とcomplex型、
統計計算では配列、
数値積分とニュートン法では反復計算、
自由落下ではファイル出力、
最後のグラフ描画では外部ライブラリを利用しています。

Fortranは画面レイアウトやGUIを作ることよりも、
数値を入力し、計算し、大量のデータを処理することを得意とする言語です。

今回の10個のコードを自分で書き換えながら試すことで、
Fortranの基本的な使い方を実際の計算へ応用できるようになります。

参考資料

• PLplot公式サイト
• Fortran-lang｜PLplot