1

ピボットせずにガウスアルゴリズムを実装しました。

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

def gauss_solve(A,b):
    """
    args: coefficient matrix A of dim(nxn) and vector b of dim(n) 
    of a system of linear equations with n unknowns.
    note: no zeroes on the main diagonal of A allowed!

    returns: vector x of dim(n) which solves the SLE
    """
    while np.ndim(A) != 2 or A.shape[0] != A.shape[1]:
        A = input(["The matrix you entered is not square, specify new input matrix A: "])
#    print "A ok."
    while np.ndim(b) != 1 or A.shape[1] != b.shape[0]:
        b = input(["The dimension of the constant vector b is incorrect, please specify new input vector b"])
#    print "b ok."
    if np.linalg.det(A) == 0:
        return "This linear system doesn't have a single unique solution."
#    print "System does have solution: "
    n = len(b)
    for i in xrange(n): # create triangular matrix
        if A[i,i] == 0:
            return "This implementation doesn't allow A to have zero entries on the main diagonal."
        A[i] = A[i]/float(A[i,i])
        b[i] = b[i]/float(A[i,i])
        for l in xrange(i+1,n):
            A[l] -= A[i]*A[l,i]
            b[l] -= b[i]*A[l,i]
    r = np.zeros(n) # result
    for i in xrange(n):
        r[-(i+1)] = b[-(i+1)] - np.dot(r,A[-(i+1)])
    return r

def test_gauss():
    m = 10
    e = 0.1
    A = sp.rand(m,m)
#    A,b = np.array([[e,1.],[1.,1.]]),np.array([1.,e])
    b = sp.rand(m)
    print gauss_solve(A,b)
    print "Build-in function says: \n", np.linalg.solve(A,b)

test_gauss()

Atest-function は、およびのランダムなエントリを生成できますb。すべてが完全にうまく機能すると思いますが、ここにマトリックスがあり、予期しない結果が生じます。

A = [[e 1] [1 1]]
b = [1 e]

e != 1解析解は

x = [-1 e+1]

しかし、いくつかの値を試してみましeたが、分析解が得られません。組み込み機能でさえsolve(A,b)失敗します。xたとえばの最初のエントリは常に0です (ただし、 と-1は完全に独立している必要がありますe)。なぜこれが起こるのか誰か説明できますか？

4

1 に答える 1

1

の新しい値を使用して更新しているためA、 との並列更新bは正しくありません。次の行を置き換える必要があります。bA

A[i] = A[i]/float(A[i,i])
b[i] = b[i]/float(A[i,i])

次のようなもので：

divisor = A[i,i]
A[i] = A[i]/float(divisor)
b[i] = b[i]/float(divisor)

同様に、次の行：

A[l] -= A[i]*A[l,i]
b[l] -= b[i]*A[l,i]

と

multiplier = A[l,i]
A[l] -= A[i]*multiplier
b[l] -= b[i]*multiplier

元のコードでは、 の行は何もしbません (浮動小数点の精度の問題は無視します): コードの最初のセクションは で除算b[i]し1.0、2 番目のセクションはから0.0時間を減算します。b[i]b[l]

于 2012-12-01T18:00:58.527 に答える