第2-3回　定常拡散方程式[python]

$$\frac{d ^2 u(x)}{d x^2}=0$$

今回解くべき方程式は上記の式であり、式変形を行うことで下記の通り変形できました。

$$\bf Au=b$$

$$\bf A= \left[\begin{array}{cccccc}1&0&&&& \\ 1&-2&1&&& \\ &\ddots& \ddots & \ddots && \\&& \ddots & \ddots & \ddots &\\ &&&1&-2&1 \\ &&&&0&1 \end{array}\right] ,\bf b= \left[\begin{array}{c} u_{0,bd} \\ 0 \\ \vdots \\ \vdots \\0 \\ u_{n-1,bd} \end{array}\right] $$

今回は以下のように条件を設定します．
x方向の長さ : \(L=1.0\)
x方向に取る点の数 : \(101\)
境界条件 : \(u(0)=0, u(L)=100\)

問題設定は次のようになります。

長さ1.0 mの物体があり、両端でそれぞれ0 Kと100 Kの物体と接しています。この条件で十分時間が経ち，定常状態に達しました。このときの温度分布はどのようになっているでしょうか．

まず解析解を境界条件から求めることができます．

$$u(x)=Ax+B$$

$$u(0)=0 ⇔ B=0$$

$$u(L)=100 ⇔ A=\frac{100}{L}$$

したがって次のようになります。

$$u(x)=\frac{100}{L}x$$

次に数値解です。結果は次のようになりました。

解析解と一致しました。今回はガウスの消去法をメインで扱いましたが、今後は他の解法も少しずつ書いていきたいと思います。

by hide

import numpy
from scipy.sparse import dia_matrix, identity
import matplotlib.pyplot as pyplot
import matplotlib.animation as animation

# -----------------------------
# -- initial value           --
# -----------------------------
x_L = 1.0                      # x_L : x方向の長さ
x_point = 101                  # x_point : ｘ方向の空間分割数
dx = x_L/(x_point-1)

# -----------------------------
# -- create matrix           --
# -----------------------------

"""
Au=bという形になるようにAとbを作成．
またプロットのため，位置xを作成.
"""

""" x """
x = []
for i in range(x_point):
    x.append(dx*i)

""" A """
data = numpy.array(
    [numpy.ones(x_point), -2.0*numpy.ones(x_point), numpy.ones(x_point)])

data[0][x_point-2] = 0
data[1][0] = 1
data[1][x_point-1] = 1
data[2][1] = 0

offsets = numpy.array([-1, 0, 1])
A = dia_matrix((data, offsets), shape=(x_point, x_point))
A = A.todense()

print("-- A --", A, "", sep='\n')

""" b """
b = numpy.array([0.0]*x_point)
b[x_point-1] = 100.0

print("-- b --", b, "", sep='\n')


# -----------------------------
# -- function                --
# -----------------------------

def gaussian_elimination(A_list, b_list):
    A_array = numpy.array(A_list)
    b_array = numpy.array(b_list)

    """ 前進消去 """
    for i in range(1, A_array.shape[0]):
        for j in range(i, A_array.shape[0]):
            hi = A_array[j, i-1] / A_array[i-1, i-1]
            A_array[j] = A_array[j] - A_array[i-1] * hi
            b_array[j] = b_array[j] - b_array[i-1] * hi

    """ 後退代入 """
    for i in range(A_array.shape[0]-2, -1, -1):
        for j in range(i, -1, -1):
            hi = A_array[j, i+1] / A_array[i+1, i+1]
            A_array[j] = A_array[j] - A_array[i+1] * hi
            b_array[j] = b_array[j] - b_array[i+1] * hi

    """ 解の算出 """
    true_x = []
    for i in range(A_array.shape[0]):
        true_x.append(b_array[i] / A_array[i, i])
    print(true_x)

    return true_x

# -----------------------------
# -- main                    --
# -----------------------------

u = gaussian_elimination(A, b)

pyplot.xlabel("length L [m]")
pyplot.ylabel("temperature [K]")

pyplot.plot(x, u, 'b-')

pyplot.show()