[Fortranのはじめ方] Part26: 連立方程式の解法プログラム作成

Fortranを使って、科学技術計算の基本となる連立一次方程式を解くプログラムを作成してみましょう。

はじめに：連立方程式とは？ 🤔

連立一次方程式は、複数の未知数を含む複数の一次方程式の組です。例えば、以下のような形をしています。

a₁₁x₁ + a₁₂x₂ + … + a_1nx_n = b₁
a₂₁x₁ + a₂₂x₂ + … + a_2nx_n = b₂
…
a_n1x₁ + a_n2x₂ + … + a_nnx_n = b_n

ここで、a_ij は係数、b_i は定数項、x_j が未知数です。科学技術計算の多くの問題（構造解析、回路網解析、流体計算など）は、最終的にこのような連立一次方程式を解くことに帰着します。Fortranは数値計算が得意なので、このような問題を解くのに適しています 💪。

このプロジェクトでは、代表的な解法の一つである「ガウスの消去法」を用いて、連立一次方程式を解くFortranプログラムを作成します。

ガウスの消去法アルゴリズム 💡

ガウスの消去法は、連立一次方程式を行列の形 (Ax = b) で表現し、行基本変形を用いて解を求める直接法の一つです。大きく分けて二つのステップからなります。

前進消去 (Forward Elimination):
係数行列 A を上三角行列（対角成分より左下がすべて0の行列）に変形します。この操作は、右辺ベクトル b にも同時に適用します。
後退代入 (Back Substitution):
変形された上三角行列を用いて、x_n から順番に x_n-1, …, x₁ と未知数の値を計算していきます。

注意点： ピボット選択（対角成分が0または非常に小さい場合に、行を入れ替える操作）を行うことで、計算の安定性を高めることができますが、今回のサンプルコードでは簡単のため省略します。

Fortranプログラム例 💻

以下に、3元の連立一次方程式をガウスの消去法で解く簡単なFortranプログラムの例を示します。

解きたい方程式系:

2x₁ + 1x₂ + 1x₃ = 5
4x₁ + 1x₂ + 0x₃ = 6
-2x₁ + 2x₂ + 1x₃ = 1

この場合、係数行列 A と右辺ベクトル b は以下のようになります。

A = [[2, 1, 1], [4, 1, 0], [-2, 2, 1]], b = [5, 6, 1]

期待される解は x = [1, 2, 1] です。


program gaussian_elimination
  implicit none
  integer, parameter :: n = 3      ! 方程式の次元数
  real(8), dimension(n, n) :: a    ! 係数行列 A
  real(8), dimension(n) :: b       ! 右辺ベクトル b
  real(8), dimension(n) :: x       ! 解ベクトル x
  real(8) :: factor, sum
  integer :: i, j, k

  ! --- 係数行列 A と右辺ベクトル b の設定 ---
  a(1,:) = [2.0d0, 1.0d0, 1.0d0]
  a(2,:) = [4.0d0, 1.0d0, 0.0d0]
  a(3,:) = [-2.0d0, 2.0d0, 1.0d0]
  b = [5.0d0, 6.0d0, 1.0d0]

  print *, "元の係数行列 A:"
  do i = 1, n
    print '(3(F8.3, X))', a(i,:)
  end do
  print *, "元の右辺ベクトル b:"
  print '(3(F8.3, X))', b
  print *

  ! === 前進消去 ===
  do k = 1, n - 1
    ! k列目の k+1 行以降を 0 にする
    do i = k + 1, n
      ! 対角成分が0の場合はエラー処理が必要だがここでは省略
      factor = a(i, k) / a(k, k)
      ! i行目から k行目の factor 倍を引く
      a(i, k:n) = a(i, k:n) - factor * a(k, k:n)
      b(i) = b(i) - factor * b(k)
    end do
  end do

  print *, "前進消去後の係数行列 A (上三角行列):"
  do i = 1, n
    print '(3(F8.3, X))', a(i,:)
  end do
  print *, "前進消去後の右辺ベクトル b:"
  print '(3(F8.3, X))', b
  print *

  ! === 後退代入 ===
  ! まず x(n) を計算
  x(n) = b(n) / a(n, n)

  ! x(n-1) から x(1) まで計算
  do i = n - 1, 1, -1
    sum = 0.0d0
    do j = i + 1, n
      sum = sum + a(i, j) * x(j)
    end do
    x(i) = (b(i) - sum) / a(i, i)
  end do

  ! --- 結果の表示 ---
  print *, "解ベクトル x:"
  print '(3(F8.3, X))', x

end program gaussian_elimination

コンパイルと実行 🚀

上記のコードを `gauss.f90` という名前で保存し、gfortran を使ってコンパイル・実行してみましょう。


# コンパイル (gfortranの場合)
gfortran gauss.f90 -o gauss -fcheck=all -Wall

# 実行
./gauss

実行すると、以下のような出力が得られるはずです（数値はフォーマットにより多少異なります）。


 元の係数行列 A:
   2.000    1.000    1.000
   4.000    1.000    0.000
  -2.000    2.000    1.000
 元の右辺ベクトル b:
   5.000    6.000    1.000

 前進消去後の係数行列 A (上三角行列):
   2.000    1.000    1.000
   0.000   -1.000   -2.000
   0.000    0.000   -4.000
 前進消去後の右辺ベクトル b:
   5.000   -4.000   -4.000

 解ベクトル x:
   1.000    2.000    1.000

期待通り、x = [1.000, 2.000, 1.000] という解が得られました 🎉。

発展：さらに学ぶために 📚

LU分解: ガウスの消去法の前進消去の過程は、行列 A を下三角行列 L と上三角行列 U の積に分解する LU分解に対応します (A = LU)。LU分解を行っておくと、同じ係数行列 A で右辺ベクトル b だけが異なる複数の連立方程式を効率的に解くことができます。
反復法: 大規模な疎行列（要素の多くが0である行列）の場合、ガウスの消去法のような直接法よりも、ヤコビ法やガウス＝ザイデル法などの反復法が効率的な場合があります。
数値計算ライブラリ: 実用的な計算では、自分で解法を実装するよりも、高度に最適化され、安定性も考慮された数値計算ライブラリを利用するのが一般的です。Fortranでよく使われる代表的なライブラリとして LAPACK (Linear Algebra PACKage) があります。LAPACKには、連立一次方程式の解法を含む、様々な線形代数計算のサブルーチンが含まれています。

これまでのステップで学んだ Fortran の知識（配列、ループ、サブルーチン、モジュールなど）を活用すれば、より洗練された、あるいは別の解法を用いたプログラムを作成することも可能です。ぜひ挑戦してみてください！

まとめ ✨

このプロジェクトでは、Fortranを用いて連立一次方程式を解く基本的なプログラムを、ガウスの消去法を例に作成しました。

連立一次方程式が科学技術計算で重要であること
ガウスの消去法のアルゴリズム（前進消去・後退代入）
Fortranでの具体的な実装方法
数値計算ライブラリ（LAPACK）の存在

などを学びました。これは、より複雑なシミュレーションやデータ解析を行う上での基礎となります。次のステップに進む前に、ぜひこのプログラムを改良したり、別の解法を試したりしてみてください！