pythonによる数値解析に触れてみよう！

どうも！ぷさいです。
（Twitterは（https://twitter.com/psy_intp_A））

みなさんはpythonに触れたことがありますか？
pythonとは、プログラミング言語の１つでコードの可読性が高く初心者におすすめの言語となっています。
（可読性が高い→コードの意味が読んでわかりやすい！ってこと）
pythonはいろいろな用途に使えるのですが、微分方程式を数値的に解くことにも使えます。
そこで、今回の記事では、pythonを用いた
「一次元非定常熱伝導方程式の数値解析」について紹介します。

今回の記事を読むと、
✅pythonを用いた微分方程式の数値解析の一例
✅実際の数値解析の問題例
✅実際のコード例
がわかります！

ただ、ちょっと専門的な内容かなと思うので、理系でない人や理系でもまだよくわからない人、プログラミング全然わからん人は「ふーーーん」とでも思ってみてください（笑）

一次元非定常熱伝導方程式とは

一次元非定常熱伝導方程式とは、次のような式です。

Tが温度、tが時間、xが空間方向、αが温度伝導率です。
非定常とは、時間で変化するという意味です。つまり∂/∂tの項が0でないということです。
微分方程式は、定常、つまり∂/∂tの項が0であるものが解きやすいです。式の形や処理が楽だからです。
まあ、微分方程式を解く際は定常か非定常かで解きやすさが大違い！とイメージしてくれればいいです。

なぜこれを取り上げたかというと３年後期の実験でこの数値解析を行ったからですね。
というわけで、俺の大学の同じ学科の後輩がもしこの記事見つけてたらラッキーですね（うらやましい）。

数値解析の問題設定

数値解析の課題は以下の通りです。

インターネットのサイトでは、一次元非定常熱伝導の温度分布のプログラミングが紹介されているものが多いですが、どれも温度分布のグラフを求めているもので、温度を数値解析で求めている例しか載ってないんですよね。あと、参考書も。

が、この問題設定はすこしひねって、
『ある特定の位置の温度がある温度になったときの時間を求めよ』
なんですよ。
つまり、元の微分方程式自体は温度を求める式だけど、その温度がある特定の値になったときの時間を求めろってことです。

わかりやすいイメージでいうと「z=f(x,y)」という関数があって、zとyが与えられてるからxを求めよという感じです。
で、ネットのサイトに載っているやり方はzを求めるやり方です。
（伝わってくれ）

なんかめんどくさそうですよね。（はいめんどくさいです）

では、コードに移ります。

数値解析のコードと解析結果

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

# パラメータ設定
L = 0.001
ρ = 7.28e3
C = 461
k = 48.0
T0 = 30
T1 = 20
T = 24
X = 0.0002
tolerance = 0.01
N = 200  # 空間分割数（変更可）

# 差分法の係数計算
α = k / (ρ * C)
Δx = L / N  
Δt = (0.3 * (Δx ** 2)) / α  # 安定性条件 (Δt <= 0.5 * Δx^2 / α) を満たすように Δt を設定した

r_x = (α * Δt) / (Δx ** 2)


# 初期条件
def u0(x):
    return T1 #薄い板は初めの状態が板全体が20 ℃である


x, u, new_u = [0.0] * (N + 2), [0.0] * (N + 2), [0.0] * (N + 2)

# 空間の離散化および各格子点に対応する温度の配列の初期値設定
for j in range(1, N + 1):
    x[j] = (j - 0.5) * Δx

    u[j] = u0(x[j])

x[0] = 0

x[N + 1] = 0.001

Z = 0  # (Zは時間をΔtだけ進めた回数）
while True:
    # 境界条件
    u[0] = 2*T0 - u[1] # ディリクレ条件（棒の先端（x=0)がT0=30で固定）
    u[N + 1] = u[N] # ノイマン境界条件（x=Lの箇所が断熱）
    # 各格子点の温度を更新
    for j in range(1, N + 1):
        new_u[j] = r_x * u[j - 1] + (1 - 2 * r_x) * u[j] + r_x * u[j + 1]

    for j in range(1, N + 1):
        u[j] = new_u[j]

    Z = Z + 1

    # 収束判定
    if abs(new_u[int(X / Δx)] - T) < tolerance:#toleranceの値設定が問題っぽい, tolerance=0.001にすると実行が終わらない
        break

# 時間を計算
time = Z * Δt
print("時間:", time, "秒")
# 各格子点の温度分布の確認（値による確認）
print("各格子点での温度[℃]:",u)
# x=0.0002 の温度確認（24 ℃)になっているか
print("x=0.0002 での温度[℃]:",u[int(X/Δx)])
# int(X/Δx）の箇所がx=0.0002 の箇所に対応しているかの確認
print("x=0.0002は何番目の格子点か=", int(X/Δx))

# 温度分布が適切か（境界条件を満たしているか）の視覚的確認
plt.plot(x, u)
plt.title("１次元熱伝導方程式（x=0.0002 mで24 ℃のとき）" , fontname="MS Gothic")
plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.show()

pythonの基本の文法については触れません。気になったら検索してみてください（例えば「python numpy」で検索とか）。
以下は、ちょっとかじってないとわからないと思うので、実行結果まで読み飛ばしても大丈夫です。

While True文と収束判定が今回のコードのカギです！
While True文は自分で設定した条件を満たさない限り、While文内の計算をループし続けるといったものです。
収束判定により、x=0.0002mの位置で24℃になったときにwhile文の繰り返しによる計算が終わるようになっています。
そして、コード上のZがWhile文の「ループを繰り返した回数」であり、１回の計算にΔt秒かかるとしているので、24℃になったときの時間tはt=Z×Δtで求められるようになっています。
（もっといい方法があるかもしれませんが、プログラミング「ド初心者」のワイにはわからん）

このコードのもう一つのポイントは「境界条件」です。
今回の問題設定では、この微分方程式の数値解を出すうえで「ディリクレ境界条件」と「ノイマン境界条件」が必要となります。
「ディリクレ境界条件」は、ある場所での温度が時間によらず一定である条件です。今回だと、x=0m（端っこ）で30℃（t>0）が条件です。
「ノイマン境界条件」は、ある場所での熱流束が時間によらず一定という条件です。つまり、∂T/∂x（温度勾配）が一定という条件です。とくにこの∂T/∂xが0である場合を「断熱条件」といいます。今回の問題だと、x=0.001m（x=0mとは逆の端っこ）で∂T/∂x=0が条件です。

また、空間の離散化には「数値流束」という考え方を使っています。
コードの
「# 空間の離散化および各格子点に対応する温度の配列の初期値設定」
のx[j]のところです！
このようなコードの書き方をすることで、境界条件のコードの書き方も楽になります。
プログラミングで微分方程式を数値的に解くには、コンピュータは微分ができないので（数値しか扱えない）、「離散化」という手法を使って、数値的に解けるようにする必要があります。

空間の離散化と境界条件のコードの書き方について詳しく知りたい人はこちら
→https://sites.google.com/site/dueyama/home-jp/python_simulation/6#h.jep7xs6wjv35

また、数値解を求めるにあたって、「Von Neumannの安定性解析」という考えを使っています。具体的には次の式で表されます。

この式は空間刻みΔxと時間刻みΔtは独立に決められないことを表しています。
今回の場合では、この式を満たしながら計算を進めるようにコードを書かないとめちゃめちゃな計算結果となってしまいます。
この条件式の導出が気になる人はこちら！！
→https://takun-physics.net/9643/

というわけで、実行結果を見てみるとこんな感じです。

実行結果

＜グラフからわかること＞
✅「Von Neumannの安定性解析」を満たして計算している
→グラフが滑らかになっている
✅グラフのx=0の位置の温度が30℃になっている
→「ディリクレ境界条件」が満たされている
✅グラフのx=0.001の位置の温度勾配∂T/∂xの傾きが０になっている
→「ノイマン境界条件（断熱条件）」が満たされている（∂T/∂x=0）

＜実行結果からわかること（画像の下）＞
✅時間：0.001917…秒
→x=0.0002 の位置での温度が24 ℃になるのにかかる時間（今回求めたかったもの）
✅x=0.0002での温度[℃]：23.990…℃
→x=0.0002で温度が24℃になっていることが確認できる（誤差0.01は許している）
✅x=0.0002は何番目の格子点か
→空間の離散化により、各位置を分割しているため、注目しているx=0.0002の位置の格子点位置が正しく設定されているか？の確認
→この実行結果はx=0.01の200分割であり、x=0.0002の位置はその1/5の位置である。よって40番目は200×1/5より、x=0.0002の格子点位置は正しく反映されていると判断できる。

以上のことから今回の問題設定の結果を正しく反映していると判断できます！
それっぽいコードが書けて実行結果がとりあえず出たら満足するのではなく、その実行結果が妥当そうか判断するための結果を出力するためのコードも書いてしっかり「自分で」チェックしないといけません。。マジ大変。。

コードを見ると簡単そうだけど、苦労したんですよ…。
だって、「数値流束」の考え方も「境界条件」のプログラミング上での実装の仕方もWhile True文使えば時間を求められるのも何も知らない・大学側は何も教えてくれない状態だからね….
コードのフォーマットも何もない。。
聞けばいいじゃんとの意見もあるかもだけど、他も忙しすぎて聞いてる暇がない。。。（質問をしっかり考える余裕が時間的余裕がない）
しかも、実行して本当にコードが妥当そうかチェックするのにも時間かかるし…。
クソ大学F○ck you！！

みなさんもぜひコードを実行してみてください！
実行環境はいろいろありますが、おすすめは「Pycharm」です。
詳しくはこの方のyoutubeに載っています！

https://www.youtube.com/watch?v=FaI8wcC1PXI

感想

これは通っている大学の３年の後期の実験で課された課題なのですが、僕が通っている大学ではプログラミングを用いた計算数理の講義がまったくないくせに課されました。
この課題出すのおかしくないか（#^ω^）これ、詰むぞ。死人がでるぞ。
→以後、死者を出したくない。
との思いで本記事を書きました（笑）
あと、書いてみると、わからない人に向けてどこまで書けばいいかわからなくなってくるなと思い、専門科目の解像度が上がれば上がるほど、それを知らない人の立場の目線の解像度が落ちやすいなと思いました。。
丁寧に書きすぎても読む側の読む気失せやすいし、手間もかかりすぎるし。。。
だとしても、この課題を課した教授はあまりにひどいと思いますけどね。
実験項目なのに、指導書もなければ実験の指示はぺらっぺらのpdf１ページだけとかなめてるだろ。。。
コードのフォーマットとかコード例ぐらい教えろや！！！！！
情報学部じゃないんですけど。。。
聞ける友達もいないので、泣きながらコードを調べて考えて書きました。
プログラミングはがうまく書けるかどうかは試行錯誤の末の運みたいなところある。こんなのやってられんわ。

本当は最後の愚痴を言いたいがために書きました（笑）
おわり！