jacobi.py

# -*- coding: utf-8 -*-
import numpy as np
import pylab as plt
import copy
 
def jacobi_or_seidel(x, a, b, use_jacobi = True):
  cpx = copy.copy(x) if use_jacobi else x
  for i in range(x.size):
    cpx[i] = b[i]
    
    for j in range(x.size):
      if i == j: continue
      cpx[i] -=  a[i][j] * x[j] 
 
    cpx[i] /= a[i][i]
  return cpx
 
if __name__ == '__main__':
  #連立一次方程式ax = bををガウス-ヤコビ法とガウス-ザイデル法で解く
  a = np.array([[2, 1], [1, 2]])
  b = np.array([4, 5])
  x = np.array([1.5, 2.5]) #xの初期値
  history = []
  for i in range(b.size):
    history.append([])
 
  for i in range(20):
    print i+1, x
    for j in range(b.size):
      history[j].append(x[j])
    x = jacobi_or_seidel(x, a, b,use_jacobi = True)
 
  x = np.linspace(0, 20, 20) #“linspace(開始値，終了値，分割数)”で，線形数列を作成。
  plt.plot(x, history[0], "-", color="k", marker="*", markersize=6)
  plt.plot(x, history[1], "--", color="k", marker="*", markersize=6)
  plt.show()