質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

2回答

3380閲覧

微分方程式を解く(デバッグ)

bigdipper_1062

総合スコア5

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2019/09/26 04:59

前提・実現したいこと

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

ポワソン方程式を解くやり方(ヤコビ法: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/ツールのバージョンなど)

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

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答2

0

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

投稿2019/09/30 02:25

bigdipper_1062

総合スコア5

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

退会済みユーザー

退会済みユーザー

2019/09/30 03:25

お役に立てて光栄です!
guest

0

ベストアンサー

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

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が出ました(一番下のフォントがないのはどうでも良いです)。

C++

1T2[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つめのエラーは共に出なかったので、質問にあるエラーは、この行のオーバーフローによるものかもしれません。
かなり適当な原因究明ですがお役に立てれば幸いです。

投稿2019/09/26 14:01

編集2019/09/26 14:06
退会済みユーザー

退会済みユーザー

総合スコア0

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

bigdipper_1062

2019/09/27 22:00

ご回答ありがとうございます。 オーバーフローについては初めて知りました。 この部分のオーバーフローをうまく処理できるようにする方法はありますでしょうか? (bigfloat パッケージなどを試しているところです。)
退会済みユーザー

退会済みユーザー

2019/09/28 01:29 編集

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問