行列の掛け算でエラーがでる numpy

Question

###問題点
whileの中での行列の計算がうまくいきません。
やりたい計算で悩んでいる部分は、漸化式がx_(n+1)=(b-Rx)/Dです。
くわしくは、ヤコビ法を見てください。
下記のようなエラーが出てきます。

yakobi.py:37: RuntimeWarning: invalid value encountered in true_divide   
  xn1=(b-np.dot(R,x))/D

エラーがでてきてもプログラム自体は動いているようなので結果を出力すると

xn1= [[ 1.   nan  inf  inf]
 [ inf -0.4 -inf  inf]
 [ inf -inf  0.   inf]
 [ inf -inf  nan  0.5]]

というものがでてきました。どうしてこのようなエラーがでてくるのか教えてください。
###コード

import numpy as np
import math
print("input n: ",end="")
n=int(input())

np.random.seed(1) #乱数の種を設定
A=np.random.randint(0,2,(n,n),int) #n行n列の行列
print("A=",A)
#対角要素にnを加える
for i in range(n):
    for j in range(n):
        if i==j:
            A[i][j]+=n
print("A=",A)

b=np.random.randint(0,2*n,n) #0以上2n未満のランダムな整数値をとる
print("b=",b)

#ヤコビ法
D=np.zeros((n,n),int)
x=np.ones((n,n))
print("x=",x)
#対角行列の生成
for i in range(n):
    for j in range(n):
        if i==j:
            D[i][j]=A[i][j]
print("D=",D)
R=A-D
print("R=",R)
i=0
while True:
    #print("i=",i)
    #xn1=np.zeros((n,n))
    xn1=(b-np.dot(R,x))/D
    #print("xn1=",xn1)
    Ax=np.dot(A,x)
    xnorm=np.linalg.norm(Ax-b)
    i+=i
    if i==100 or xnorm<=pow(10,-6):
        break
print("x=",xn1)

Answer 1

D に値が0の要素があるため、xn1=(b-np.dot(R,x))/D で0除算が発生して、エラーになっています。

print("D=",D) の結果を見るとわかりますが、D の対角成分以外は0になっていますね。

修正方法

以下をDの逆行列の行列積に変更してください

- xn1=(b-np.dot(R,x))/D
+ x = np.dot(np.linalg.inv(D), (b - np.dot(R, x)))

あと以下もおかしいですね

- i += i
+ i += 1

追記

• 行列ベクトル積、行列積は numpy.dot() を使うほかに @ 演算子でも計算可能です。

np.dot(A, x) → A @ x
np.dot(A, B) → A @ B
@ 演算子を使ったほうがプログラムの見通しがよくなるかと思います。

• 1 / D だと対角成分以外の除算を行った場合に0割が発生するので、np.linalg.inv(D) で逆行列を計算するか、D_inv = np.reciprocal(D, out=np.zeros(D.shape, float), where=D != 0) で0でない要素だけ逆数をとります。
• pow(10,-6) は、指数表記 1e-6 でも表せます

R = A - D  # 対角成分以外
D_inv = np.linalg.inv(D)

x = np.zeros_like(b)  # 初期値
for i in range(100):
    x = D_inv @ (b - R @ x)

    if np.linalg.norm(A @ x - b) <= 1e-6:
        break

Comments

[bsdfan]

書き間違いだと思いますが、pow(10, -6) は、10e-6 ではなく 1e-6 かと。

[tiitoi]

ご指摘ありがとうございます。修正しました。

[langhtorn]

ご丁寧に見やすい書き方まで教えてくださりありがとうございます。おかげで解決しました。