singular matrix

Question

一次元拡散方程式のシミュレーションのプログラミングを作成しています。一通りシミュレーションの計算までプログラムを組んで実行しようとしたところ、「singular matrix」というエラーが発生してしまいました。行列に問題があるときにでるエラーとのことですが、どこに問題があるかさっぱりわかりません。どなたかご指摘いただけないでしょうか？

import numpy as np
import matplotlib.pyplot as plt
import math
import matplotlib.animation as animation

T = 2
nx = 41
nt = 201
dx = math.pi/(nx-1)
dt = 0.01
du = 0.3
alpha = 0.1
beta = 1.0
gamma = 10.0
X = np.linspace(0,math.pi,nx)

#初期条件
u1 = []
u2 = []
u3 = []
for i in range(nx):
  x = dx*i
  u1.append(math.cos(x)+math.cos(7*x))
  u2.append(math.cos(x)+math.cos(7*x))
  u3.append(math.cos(x)+math.cos(7*x))

fig = plt.figure()
ims = []
for t in range(0,nt):
  for i in range(1,nx-1):
    u1[i] = u1[i]+alpha*dt/(dx*dx)*(u1[i+1]-2*u1[i]+u1[i-1])
    u2[i] = u2[i]+beta*dt/(dx*dx)*(u2[i+1]-2*u2[i]+u2[i-1])
    u3[i] = u3[i]+gamma*dt/(dx*dx)*(u3[i+1]-2*u3[i]+u3[i-1])
    im1 = plt.plot(X,u1,label="D=0.1")
    im2 = plt.plot(X,u2,label="D=1.0")
    im3 = plt.plot(X,u3,label="D=10.0")
    ims.append(im1+im2+im3)


plt.grid()
plt.xlabel("x")
plt.ylabel("u")
anim = animation.ArtistAnimation(fig, ims)
plt.show()

Comments

[jbpb0]

エラーが出るのは、コードのどこを実行した時ですか？

[lucy-0333]

それが、プログラムを実行すると同じようなエラーが無数に出てくる上に、line1782,line97,line158,line2176…などと、どこからわからないlineばかりが指定されて、どこでエラーが起きてるかすらわからない状況です…

[jbpb0]

fig =...以降のグラフ作る部分をコメント化して、下記のように変えてください
(インデントは元のを維持したまま)

#fig = plt.figure()
#ims = []
for t in range(0,nt):
for i in range(1,nx-1):
u1[i] = u1[i]+alpha*dt/(dx*dx)*(u1[i+1]-2*u1[i]+u1[i-1])
u2[i] = u2[i]+beta*dt/(dx*dx)*(u2[i+1]-2*u2[i]+u2[i-1])
u3[i] = u3[i]+gamma*dt/(dx*dx)*(u3[i+1]-2*u3[i]+u3[i-1])
#im1 = plt.plot(X,u1,label="D=0.1")
#im2 = plt.plot(X,u2,label="D=1.0")
#im3 = plt.plot(X,u3,label="D=10.0")
#ims.append(im1+im2+im3)
#plt.grid()
#plt.xlabel("x")
#plt.ylabel("u")
#anim = animation.ArtistAnimation(fig, ims)
#plt.show()

上記の下に、下記を追加してください

print(u1)
print(u2)
print(u3)

そして実行して、上記追加したprint()の結果を見てください

[lucy-0333]

確認したところ、u3について[2.0,nan,nan,nan…,]という形で欠損値？のようなものが生まれていました。

[lucy-0333]

わかりました。安定条件を考えずにやっていたのが原因のようです。何とか解決させることができました！ありがとうございました！

Answer 1

fig =...以降のグラフを作る部分を全部コメント化して、下記のように変えてください

#fig = plt.figure()
#ims = []
for t in range(0,nt):
  for i in range(1,nx-1):
    u1[i] = u1[i]+alpha*dt/(dx*dx)*(u1[i+1]-2*u1[i]+u1[i-1])
    u2[i] = u2[i]+beta*dt/(dx*dx)*(u2[i+1]-2*u2[i]+u2[i-1])
    u3[i] = u3[i]+gamma*dt/(dx*dx)*(u3[i+1]-2*u3[i]+u3[i-1])
    #im1 = plt.plot(X,u1,label="D=0.1")
    #im2 = plt.plot(X,u2,label="D=1.0")
    #im3 = plt.plot(X,u3,label="D=10.0")
    #ims.append(im1+im2+im3)

#plt.grid()
#plt.xlabel("x")
#plt.ylabel("u")
#anim = animation.ArtistAnimation(fig, ims)
#plt.show()

上記の下に、下記を追加してください

print(u1)
print(u2)
print(u3)

そして実行して、上記追加したprint()の結果を見てください
そうすれば、何が悪いか分かると思います