微分方程式を解く（デバッグ）

Question

前提・実現したいこと

現在、python3で以下のような２次元熱拡散方程式を解くプログラムを作っています。

ポワソン方程式を解くやり方（ヤコビ法：URL）を参考に解いてみようとしたのですが、ポワソン方程式では電荷の項が定数であるのに対し、添付の方程式では、T^4 の項が含まれており、おそらくここの処理でうまく行っていないようです。

https://qiita.com/sci_Haru/items/6b80c7eb8d4754fb5c2d

初歩的な質問かもしれませんが、プログラムがうまく実行されない原因について教えていただけると幸いです。

また、このやり方で画像のような方程式が解けるのかについて自信がないので、そもそもの方法が間違っているなどであればご指摘いただけると助かります。

発生している問題・エラーメッセージ

次のようなエラーメッセージが出て、メインのところでiterationが行われていないようです。
invalid value encountered in subtract

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation # アニメーション作成のためのメソッドをインポート

anim = [] 

# 
delta_L=0.01 #グリッド幅
LL = 1 # 正方形の幅
L = int(LL/delta_L)

V = 300.0 # 境界条件
convegence_criterion = 10**-5

T = np.zeros([L+1,L+1])
T[0,:] = V # 境界条件
T2 = np.empty ([L+1,L+1])


e = 1.                   #輻射率
sb = 5.67               #ステファンボルツマン定数(W/m^2*K^4)
T0 = 300                #バックグラウンド室温(K) 
k = 72                  #白金の熱伝導率(W/m*K) 
tf = 2.5*10**(-6)       #箔の厚さ(m)
Prad = 5.9*10**(-3)     #レーザーパワー強度(W)
    
#レーザーパワーの設定
rad = np.zeros([L+1,L+1])
laser_pix = 2

CC=int(L/2 - laser_pix/2)
CC2=int(L/2 + laser_pix/2)

# レーザーパワーを埋め込む
for i in range(CC, CC2+1):
    for j in range(CC,CC2+1):
        rad[i,j] = Prad/(k*tf)
        
#rad = np.empty([L+1,L+1])
#rad[:,:] = Prad/(k*tf)

#for plot
im=plt.imshow(T,cmap='hsv')
anim.append([im])

# メイン
delta = 1.0
n_iter=0
conv_check=[]
while delta > convegence_criterion:
    if n_iter % 100 ==0:  # 収束状況のモニタリング
        print("iteration No =", n_iter, "delta=",delta)
    conv_check.append([n_iter, delta])
    for i in range(L+1):
        for j in range(L+1):
            if i ==0 or i ==L or j==0 or j==L:
                T2[i,j] = T[i,j]
            else: 
                T2[i,j] = (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1])/4 + e*sb*(T[i,j]**4 - T0**4)/(k*tf) + rad[i,j]
    delta = np.max(abs(T-T2))

    T, T2 = T2, T
    n_iter+=1

    im=plt.imshow(T,cmap='hsv')
    anim.append([im])


#for plot    

pp=plt.colorbar (orientation="vertical") # カラーバーの表示 
pp.set_label("T", fontname="Ariel", fontsize=24)
plt.xlabel('X',fontname="Ariel", fontsize=24)
plt.ylabel('Y',fontname="Ariel", fontsize=24)
#ani = ArtistAnimation(fig, anim, interval=10, blit=True,repeat_delay=1000)
#ani.save("t.gif", writer='imagemagick')  

plt.show()

試したこと

型が一致しているかの確認

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Answer 1

ご回答を受けて、変数の値（境界条件等）を色々変えてみたところうまく実行できたようです。
オーバーフローは知らなかったので大変助かりました。

Comments

[退会済みユーザー]

お役に立てて光栄です！

Answer 2

該当のソースコードをこちらで実行してみたところ、

iteration No = 0 delta= 1.0
test.py:59: RuntimeWarning: overflow encountered in double_scalars
T2[i,j] = (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1])/4 + esb(T[i,j]4 - T04)/(k*tf) + rad[i,j]
test.py:60: RuntimeWarning: invalid value encountered in subtract
delta = np.max(abs(T-T2))
***\matplotlib\font_manager.py:1241: UserWarning: findfont: Font family ['Ariel'] not found. Falling back to DejaVu Sans.
(prop.get_family(), self.defaultFamily[fontext]))

とwarningが出ました(一番下のフォントがないのはどうでも良いです)。

T2[i,j] = (T[i+1,j] + T[i-1,j] + T[i,j+1] + T[i,j-1])/4 + e*sb*(T[i,j]**4 - T0**4)/(k*tf) + rad[i,j]

(質問にはございませんが一応)こちらでオーバーフローしているようです。

それと、この行をT2[i,j] = 1に変更したところ、1つめ、2つめのエラーは共に出なかったので、質問にあるエラーは、この行のオーバーフローによるものかもしれません。
かなり適当な原因究明ですがお役に立てれば幸いです。

Comments

[bigdipper_1062]

ご回答ありがとうございます。
オーバーフローについては初めて知りました。

この部分のオーバーフローをうまく処理できるようにする方法はありますでしょうか？
（bigfloat パッケージなどを試しているところです。）

[退会済みユーザー]

うまく処理する方法について色々考えてみたのですが、私自身この方程式をしっかり理解できていないこともあり、うまい方法は思いつきませんでした。
ただ、このTやT2が物体の温度であるとすると(そう読み取りましたがそれは正しいでしょうか…)その値がオーバーフローするのは到底おかしなことだと思いますので、もしかしたら式に問題があるのかもしれないと思いました(デバッグとしてそのT、T2の値を見てみたのですがとてつもなく大きくなっていました)。