Jacobi method
from pprint import pprint
from numpy import array, zeros, diag, diagflat, dot
def jacobi(A,b,N=25,x=None):
"""Solves the equation Ax=b via the Jacobi iterative method."""
# Create an initial guess if needed
if x is None:
x = zeros(len(A[0]))
# Create a vector of the diagonal elements of A
# and subtract them from A
D = diag(A)
R = A - diagflat(D)
# Iterate for N times
for i in range(N):
x = (b - dot(R,x)) / D
return x
A = array([[2.0,1.0],[5.0,7.0]])
b = array([11.0,13.0])
guess = array([1.0,1.0])
sol = jacobi(A,b,N=25,x=guess)
print "A:"
pprint(A)
print "b:"
pprint(b)
print "x:"
pprint(sol)
输出结果
A:
array([[ 2., 1.],
[ 5., 7.]])
b:
array([ 11., 13.])
x:
array([ 7.11110202, -3.22220342])
