タイトルそのままです．偏微分方程式をPythonで解きました．

問題

円柱座標系で，$z$軸に印加された磁場$B_z$中に静止した無限長の半径$a$，抵抗率$\eta$の円柱プラズマを考える．このプラズマに$z$軸方向の電流$I_z$を流す時，方位角方向の磁場$B_\theta$が発生する．プラズマに浸透する$B_\theta$と電流密度$J_z$の時空間発展を求める．
（金属に電流をいきなり流した時の過渡応答と同じになるっぽい．いわゆる表皮効果）

解きたい微分方程式

〜〜導出略〜〜

$$
\frac{\partial B_\theta}{\partial t} = \frac{\eta}{\mu_0}\frac{\partial}{\partial r}\left(\frac{1}{r} \frac{\partial (rB_\theta)}{\partial r}\right)
$$
$\mu_0$は真空の透磁率．また電流密度$J_z$は$\nabla \times B=\mu_0 J_z$より，
$$
J_z = \frac{1}{\mu_0r}\frac{\partial}{\partial r}(rB_\theta)
$$
$r>a$では$J_z=0$なので$B_\theta(r)=\mu_0I_z/2\pi r$．また，定常状態（$\partial/\partial t=0$）で，左辺を解くと，$B_\theta(r)=(\mu_0I_z/2\pi a^2)r$になる（$r=a$で$B_\theta(r)=\mu_0I_z/2\pi a$）．$r$を$a$，$t$を$a^2\mu_0/\eta$，$B_\theta(r)$を$B_\theta(a)$，$J_\theta$を$I_z/\pi a^2$でそれぞれ規格化して整理すると，

\frac{\partial B_\theta}{\partial t} = \frac{\partial}{\partial r}\left(\frac{1}{r} \frac{\partial (rB_\theta)}{\partial r}\right)\\
= \left\{\frac{\partial^2}{\partial r^2} + \frac{1}{r} \frac{\partial}{\partial r} -\frac{1}{r^2} \right\}B_\theta

$$
J_z = \frac{1}{2}\frac{\partial}{\partial r}(rB_\theta)
$$

これを数値計算で解く．

離散化

以下，添字が煩雑なので$B_\theta$を$B$と書きます．時刻$t^j$，位置$r_i$における値を$B_i^j$，時間と距離の幅をそれぞれ$\Delta t$, $\Delta r$とする．
左辺について，
$$
(\text{左辺})=\frac{B_i^{j+1}-B_i^{j}}{\Delta t}
$$

右辺について
1階中心差分商と，
$$
\frac{\partial B}{\partial r} \approx \frac{B_{i+1}^j-B_{i+1}^j}{2\Delta r}
$$
2階中心差分商
$$
\frac{\partial B}{\partial r} \approx \frac{B_{i+1}^j-2B_{i}^j+B_{i-1}^j}{\Delta r^2}
$$
を使って書き下すと，
$$
(\text{右辺})=\frac{B_{i+1}^j-2B_{i}^j+B_{i-1}^j}{\Delta r^2}
+ \frac{1}{r_i}\frac{B_{i+1}^j
- B_{i+1}^j}{2\Delta r}-\frac{1}{r_i^2}B_i^j\
= \left(\frac{1}{\Delta r^2}-\frac{1}{2\Delta r}\right)B_{i-1}^j
- \left(-\frac{2}{\Delta r^2}+\frac{1}{r_i^2}\right)B_{i}^j
+ \left(\frac{1}{\Delta r^2}+\frac{1}{2\Delta r}\right)B_{i+1}^j\
=a_iB_{i-1}^j - b_iB_{i}^j + c_iB_{i+1}^j
$$
ここで$a_iB_{i-1}^j - b_iB_{i}^j + c_iB_{i+1}^j$を$\delta B_i^j$と定義する．

ここで，左辺と右辺の添字$j$の取り方で前進差分と後進差分がある．

前進差分

$j$を$\left(B_i^{j+1}-B_i^{j}\right)/\Delta t = \delta B_i^j$ととる方法で，時刻$j$の値から$j+1$を計算する（陽解法）．$B_i^{j}$を移項すると漸化式になるので，計算が簡単．
しかし，$\Delta r$に対して$\Delta t$を細かくとらないと数値的に不安定になる．

後進差分

$j$を$\left(B_i^{j+1}-B_i^{j}\right)/\Delta t = \delta B_i^{j+1}$ととる方法．時刻$j$の値を満たすような$j+1$での値を計算する（陰解法）．連立方程式を解く必要があり，ちょっとしんどい．

θ法

前進差分と後進差分のブレンド．
$$
\left(B_i^{j+1}-B_i^{j}\right)/\Delta t = (1-\theta)\delta B_i^j + \theta\delta B_i^{j+1}
$$
時間方向に2次精度になるらしい．参考

コード

差分方程式解く部分

import numpy as np
from tqdm import tqdm

def theta_method(n=10, m=100, theta=0.0, t_max=1.0):
    n = int(n)
    m = int(m)
    r, h = np.linspace(0.0001, 1, n + 1, retstep=True)
    t, tau = np.linspace(0, t_max, m + 1, retstep=True)

    a = 1 / h ** 2 - 1 / (2 * r * h)
    b = 2 / h ** 2 + 1 / r ** 2
    c = 1 / h ** 2 + 1 / (2 * r * h)

    B_theta = np.zeros((n + 1, m + 1))
    w = np.zeros(n + 1)
    A = np.zeros((n + 1, n + 1))
    A[0, 0] = 1
    A[n, n] = 1

    for i in range(1, n):
        A[i, i - 1] = -tau * theta * a[i]
        A[i, i] = 1 + tau * theta * b[i]
        A[i, i + 1] = -tau * theta * c[i]

    B_theta[:, 0] = np.zeros_like(r)

    for j in tqdm(range(m)):
        w[0] = 0
        w[n] = 1
        w[1:n] = (1 - theta) * tau * a[1:n] * B_theta[:-2, j] \
                 + (1 - (1 - theta) * tau * b[1:n]) * B_theta[1:n, j] \
                 + (1 - theta) * tau * c[1:n] * B_theta[2:, j]
        v = np.linalg.solve(A, w)
        for i in range(n + 1):
            B_theta[i, j + 1] = v[i]

    return r, h, t, tau, B_theta, np.apply_along_axis(lambda B: calc_J_z(B, r), 0, B_theta)

def calc_J_z(B_theta, r):
    return B_theta / r + np.gradient(B_theta, r)

描画

if __name__ == '__main__':
    while True:
        cmd = input("quit(q) or [enter]>>")
        if cmd == "q":
            break
        else:
            # t_list = np.array(list(map(float, input("t>>").split(","))))
            t_list = np.array([0.0001, 0.001, 0.01, 0.1, 0.2, 0.5, 1])
            n = int(input("n>>"))
            m = int(input("m>>"))
            theta = float(input("theta>>"))

            cmd = input(f"lambda={n ** 2 * t_list.max() / m:.2} continue?(Y/n)>>")
            if cmd == "n":
                continue

            r, h, t, tau, B_theta, J_z = theta_method(n=n, m=m, theta=theta, t_max=t_list.max())
            ind = np.round(t_list / tau).astype(int)
            B_theta_list = B_theta.T[ind]
            J_z_list = J_z.T[ind]
            B_theta_inf = r
            J_z_inf = calc_J_z(r, r)

            fig, axes = plt.subplots(1, 3, figsize=(12, 6), sharex=True)
            fig.suptitle(rf"$n={n}$, $m={m}$, $\lambda={tau / h ** 2:.2}$, $\theta={theta:.2}$")
            axes[0].set_title(r"$B_\theta$")
            axes[0].set_ylabel(r"$B_\theta$")
            axes[0].set_ylim(0.0, 1.1)

            axes[1].set_title(r"$B_\theta$(steady)$-B_\theta$")
            axes[1].set_ylabel(r"$B_\theta$")
            axes[1].set_ylim(0.0, 1.1)

            axes[2].set_title(r"$J_z$")
            axes[2].set_ylabel(r"$J_z$")
            axes[2].set_ylim(0., 5)

            axes[0].set_xlim(0.0, 1.1)
            axes[0].plot(r, B_theta_inf, ls="--", c="k", label="steady")
            axes[1].plot(r, B_theta_inf, ls="--", c="k", label="steady")
            axes[2].plot(r, J_z_inf, ls="--", c="k", label="steady")

            for B, J, t_ in zip(B_theta_list, J_z_list, t_list):
                axes[0].plot(r, B, label=f"t={t_:.2}")
                axes[1].plot(r, B_theta_inf - B, label=f"t={t_:.2}")
                axes[2].plot(r, J, label=f"t={t_:.2}")

            for ax in axes:
                ax.legend()
                ax.set_xlabel(r"$r$")
                ylim = ax.get_ylim()
                ax.vlines(1, ylim[0], ylim[1], ls="-", colors="k", lw=1)
                ax.set_ylim(0, ylim[1])

            fig.tight_layout()
            fig.savefig(f"B_theta_penetration_{n}_{m}_{theta}.png", bbox_inches="tight", pad_inches=0.05, dpi=200)
            plt.show()

            # animation
            fig, axes = plt.subplots(1, 3, figsize=(12, 6), sharex=True)

            # pick up 100 times from t_list in log scale
            t_list_animation = np.unique(np.logspace(0, np.log10(m), MAX_FRAMES, dtype=int))
            frame_number = len(t_list_animation)
            skip = m // frame_number


            def update(frame):
                j = t_list_animation[frame]
                axes[0].cla()
                axes[1].cla()
                axes[2].cla()

                axes[0].set_title(rf"$B_\theta$ ($t={t[j]:.2}$)")
                axes[0].set_ylabel(r"$B_\theta$")
                axes[0].set_ylim(0.0, 1.1)

                axes[1].set_title(rf"$B_\theta$(steady)$-B_\theta$ ($t={t[j]:.2}$)")
                axes[1].set_ylabel(r"$B_\theta$")
                axes[1].set_ylim(0.0, 1.1)

                axes[2].set_title(rf"$J_z$ ($t={t[j]:.2}$)")
                axes[2].set_ylabel(r"$J_z$")
                axes[2].set_ylim(0., 5)

                axes[0].plot(r, B_theta.T[j])
                axes[1].plot(r, B_theta_inf - B_theta.T[j])
                axes[2].plot(r, J_z.T[j])

                axes[0].set_xlim(0.0, 1.1)
                axes[0].plot(r, B_theta_inf, ls="--", c="k", label="steady")
                axes[1].plot(r, B_theta_inf, ls="--", c="k", label="steady")
                axes[2].plot(r, J_z_inf, ls="--", c="k", label="steady")

                for ax in axes:
                    ax.legend()
                    ax.set_xlabel(r"$r$")
                    ylim = ax.get_ylim()
                    ax.vlines(1, ylim[0], ylim[1], ls="--", colors="k")
                    ax.set_ylim(0, ylim[1])

            print("making animation...")
            anim = FuncAnimation(fig, update, frames=range(frame_number), interval=200)
            print("saving animation...")
            anim.save(f"B_theta_penetration_{n}_{m}_{theta}.gif", writer="pillow", fps=10)
            print("done.")