gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

u'''
  連立方程式の解法
    テスト用の式
      2*x1 + 3*x2 +   x3 =  4
      4*x1 +   x2 - 3*x3 = -2
       -x1 + 2*x2 + 2*x3 =  2
'''

import unittest
from math import fabs


class TestEquation(unittest.TestCase):

    EPS = 1e-8

    def testEquation(self):
        def is_list_below_eps(l1, l2):
            for v1, v2 in zip(l1, l2):
                if fabs(v1 - v2) > self.EPS:
                    return False
            return True

        # 使用する方程式です
        expressions = ((2.0,  3.0,  1.0,  4.0),
                       (4.0,  1.0, -3.0, -2.0),
                       (-1.0, 2.0,  2.0,  2.0))
        # 方程式の解です
        answer = (2.0, -1.0, 3.0)

        # 各関数の結果です
        result1 = gauss_jordan(expressions)
        result2 = pivot(expressions)
        result3 = gaussian_elimination(expressions)

        # 本テストケースの本体
        self.assert_(is_list_below_eps(answer, result1))
        self.assert_(is_list_below_eps(answer, result2))
        self.assert_(is_list_below_eps(answer, result3))


def gauss_jordan(expressions):
    u'''
        ガウス=ジョルダン法
         ＼ l →
        i,k a11 * x1 + a12 * x2 + a13 * x3 = b1 …  (1)
        ↓   a21 * x1 + a22 * x2 + a23 * x3 = b2 …  (2)
            a31 * x1 + a32 * x2 + a33 * x3 = b3 …  (3)

            一つ目の式の第一項から項の選択を斜め進めます。
            a[i][i] と同じ項で他式を減算し、
            他の式からその項を排除する。
    '''
    results = map(lambda t: list(t), expressions)
    for i in xrange(len(expressions)):
        k = results[i][i]
        for l in xrange(len(results[i])):
            # a[i][i] == 1 になるように行を除算する
            results[i][l] /= k
        for j in xrange(len(results)):
            # a[i] で他の式を減算する
            if i == j:
                continue
            k = results[j][i]
            for l in xrange(len(results[j])):
                results[j][l] -= results[i][l] * k
    return [d[-1] for d in results]


# ピボット選択法
def pivot(expressions):
    u'''
        ピボット選択法
            |0| から最も遠い係数の項を選択し、
            他の式からその項を排除する。
    '''
    results = map(lambda t: list(t), expressions)
    for i in xrange(len(expressions)):
        k = m = 0
        for l, v in enumerate(results[i][:-1]):
            if abs(k) < abs(v):
                m, k = l, v
        for l in xrange(len(results[i])):
            # a[i][i] == 1 になるように行を除算する
            results[i][l] /= k
        for j in xrange(len(results)):
            # a[i] で他の式を減算する
            if i == j:
                continue
            k = results[j][m]
            for l in xrange(len(results[j])):
                results[j][l] -= results[i][l] * k

    answer = [0] * (len(results[0]) - 1)
    for i in xrange(len(results)):
        for l in xrange(len(results[i]) - 1):
            if results[i][l]:
                answer[l] = results[i][-1]
    return answer


# ガウスの消去法
def gaussian_elimination(expressions):
    """
        ガウスの消去法

        未知数 x1, x2, x2 についての式 (1), (2), (3)

        # 式
        a11 * x1 + a12 * x2 + a13 * x3 = b1 ... (1)
        a21 * x1 + a22 * x2 + a23 * x3 = b2 ... (2)
        a31 * x1 + a32 * x2 + a33 * x3 = b3 ... (3)

        # 前進消去
        (2)' = (2) - (1)*(a21/a11)
        (3)' = (3) - (1)*(a31/a11)

        (3)'' = (3)' - (2)'*(a32'/a22')

        # 後退代入
        x3 = b3''/a33''
        x2 = (b2' - a23'*x3) / a21
        x1 = (b1 - a12*x2 - a13*x3) / a11
    """
    results = [list(_) for _ in expressions]
    pivots = []

    # 前進消去
    for i in xrange(len(results) - 1):
        value = max(results[i][:-1], key=lambda x: fabs(x))
        pivot = results[i].index(value)
        pivots.append(pivot)
        for j in xrange(i + 1, len(results)):
            k = results[j][pivot] / results[i][pivot]
            for l in xrange(0, len(results[j])):
                results[j][l] -= results[i][l] * k
    value = max(results[-1][:-1], key=lambda x: fabs(x))
    pivot = results[-1].index(value)
    pivots.append(pivot)

    # 後退代入
    for i in range(len(results))[::-1]:
        for l, j in list(enumerate(pivots))[i + 1:]:
            results[i][-1] -= results[i][j] * results[l][-1]
            results[i][j] = 0
        results[i][-1] /= results[i][pivots[i]]
        results[i][i] = 1

    return [results[pivots[i]][-1] for i in xrange(len(results))]


def main():
    # a11 * x1 + a12 * x2 + a13 * x3 = b1
    # a21 * x1 + a22 * x2 + a23 * x3 = b2
    # a31 * x1 + a32 * x2 + a33 * x3 = b3
    expressions = ((2.0,  3.0,  1.0, 4.0),
                   (4.0,  1.0, -3.0, -2.0),
                   (-1.0, 2.0,  2.0, 2.0))
    print 'Expressions:'
    print expressions
    print 'Result:'
    print gauss_jordan(expressions)
    print pivot(expressions)
    print gaussian_elimination(expressions)


if __name__ == '__main__':
    unittest.main()