C# で Gauss-Seidel 法 によるn元連立1次方程式の単純な solver 実装

GaussSeidel.cs

using System;

namespace MathUtil
{
    public class Solver
    {
        static void Main(string[] args)
        {
            // 方程式
            // 3x + 2y + z = 10 
            // x + 4y + z = 12
            // 2x + 2y + 5z = 21

            // 解 
            // (x, y, z) = (1, 2, 3)

            var matrix = new double[,] {
                { 3, 2, 1 },
                { 1, 4, 1 },
                { 2, 2, 5 }
            };
            var b = new double[] { 10, 12, 21 };
            var eps = 1e-15;
            var result = Solver.GaussSeidel(matrix, b, 10000, eps);


            Console.WriteLine("Error:" + result.Error);
            Console.WriteLine("Count:" + result.Iterator);

            var xyz = new[] { "x", "y", "z" };
            for (int i = 0; i < result.Solution.Length; i++)
            {
                Console.WriteLine(xyz[i] + " = " + result.Solution[i]);
            }
        }

        /// <summary>
        /// <para>
        ///  ガウス＝ザイデル法（Gauss-Seidel method）
        /// </para>
        /// <para>
        /// 解が収束するのは
        ///    ・対角有利(diagonal dominant, 対角要素の絶対値>その行の他の要素の絶対値の和)
        ///    ・係数行列が対称(symmetric)かつ正定(positive definite)
        ///    ・Σ_j |a_ij/a_ii| &lt; 1 (i = 1～n, j != i) 
        /// </para>
        /// </summary>
        /// <param name="squareMatrix">n×nの係数行列</param>
        /// <param name="constantVector">n×1の定数項(b)の行列(定数項ベクトル)</param>
        /// <param name="maxIterator">最大反復数</param>
        /// <param name="eps">許容誤差</param>
        /// <param name="AbsoluteError">[option]収束判定で絶対誤差と相対誤差のどちらを用いるか示す。
        /// 真のときは絶対誤差、偽のときは相対誤差を用いる。</param>
        /// <returns>解の行列</returns>
        public static IterativeResult GaussSeidel(double[,] squareMatrix, double[] constantVector, int maxIterator, double eps, bool AbsoluteError = true)
        {
            if (squareMatrix.GetLength(0) != squareMatrix.GetLength(1))
            {
                throw new ArgumentException("引き数の係数行列が正方行列でありません。", "A");
            }
            if (squareMatrix.GetLength(0) != constantVector.Length)
            {
                throw new ArgumentException("引き数の定数項行列が係数行列の大きさと一致しません。");
            }

            // 行列の大きさ
            int n = squareMatrix.GetLength(0);
            // 解。初期値はすべて0
            double[] solution = new double[n];
            // 誤差
            double e = 0.0;
            // 現在の反復回数
            int k;

            double tmp;

            for (k = 0; k < maxIterator; ++k)
            {
                // 現在の値を代入して次の解候補を計算
                e = 0.0;
                for (int i = 0; i < n; ++i)
                {
                    tmp = solution[i];
                    solution[i] = constantVector[i];
                    for (int j = 0; j < n; ++j)
                    {
                        solution[i] -= (j != i ? squareMatrix[i, j] * solution[j] : 0.0);
                    }
                    solution[i] /= squareMatrix[i, i];

                    if (AbsoluteError)
                    {
                        // 絶対誤差
                        e += Math.Abs(tmp - solution[i]);
                    }
                    else
                    {
                        // 相対誤差
                        e += Math.Abs((tmp - solution[i]) / tmp);
                    }
                }
                // 収束判定
                if (e <= eps)
                {
                    break;
                }
            }

            return new IterativeResult(solution, k, e);
        }

        public struct IterativeResult
        {
            public IterativeResult(double[] solution, int iterator, double error)
            {
                this.Solution = solution;
                this.Iterator = iterator;
                this.Error = error;

            }

            /// <summary>
            /// 解
            /// </summary>
            public double[] Solution { get; set; }

            /// <summary>
            /// 反復回数
            /// </summary>
            public int Iterator { get; set; }

            /// <summary>
            /// 誤差
            /// </summary>
            public double Error { get; set; }
        }
    }
}