2차원 Heat Equation의 풀이 #3 : OSM 예제풀이2

 

 

이 글은

 

2023.02.09 - [수학의 재미/아름다운 이론] - 2차원 Heat Equation의 풀이 #2 : OSM 예제풀이

2차원 Heat Equation의 풀이 #2 : OSM 예제풀이

 

에 이어, 재미있는 예제를 하나 소개할까 합니다.

 

 

 

2차원 Heat Equation 예제

 

다음의 열방정식을 풀어봅시다.

$$\frac{\partial u}{\partial t} = \beta \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) ~~,~~ 0\leq t\leq 10 ~,~  0\leq x \leq 1 ~,~ 0\leq y \leq 1$$ ○ 경계조건 : $ u(t, x, 0) =30 , u(t,x,1) =100, u(t,0,y) =0,  u(t,1,y)=0$ ○ 초기조건:  $u(0,x,y) = 0, 0<x,y<1$

이 방정식의 정확한 해는 잘 모르지만, 앞서 보았던 OSM 방법으로 해결할 수 있습니다.

위 방정식의 경계조건을 그려보면 다음과 같습니다.

 

 

 

 

따라서 위 방정식은,

 

초기온도가 모두 0인 정사각형 영역에서,
- 왼쪽, 오른쪽에는 항상 0도로 유지되는 금속 막대기를 설치하고, 
- 위쪽에는 항상 100도로 유지되는 금속 막대기를 설치
- 아래쪽에는 항상 30도로 유지되는 금속 막대기를 설치
했을 때, 

주어진 영역안에서 온도가 어떻게 퍼질까? 를 구해보는 것입니다.

 

OSM FDM 을 쓰기 위한 이산화

 

지난 글의 예제 와는 달리, $\beta$라는 상수와, 경계 조건에 0이 아닌 값이 포함되어 있어,  OSM을 위한 Tridiagonal 행렬을 만들 때 주의할 점이 있습니다. 

 

$x$축 방향으로의 풀이

고정된 $j,~ 0\leq j \leq J$에 대해,

$$ \frac{u_{ij}^{n+1/2}-u_{ij}^{n}}{k} = \beta \frac{u_{i-1,j}^{n+1/2}-2u_{ij}^{n+1/2}+u_{i+1,j}^{n+1/2}}{h_x^2} $$  

 

따라서

$$ \frac{u_{ij}^{n+1/2}-u_{ij}^{n}}{k} = \alpha_x \beta ( u_{i-1,j}^{n+1/2}-2u_{ij}^{n+1/2}+u_{i+1,j}^{n+1/2} ) $$

즉,
$$ -\alpha_x \beta u_{i-1,j}^{n+1/2} +(1+2\alpha_x\beta) u_{ij}^{n+1/2} -\alpha_x \beta u_{i+1,j}^{n+1/2} =  u_{ij}^{n}$$

 

입니다. 따라서 $i=1$ 일 때는,

 

$$ -\alpha_x \beta u_{0,j}^{n+1/2} + (1+2\alpha_x\beta) u_{1,j}^{n+1/2} - \alpha_x \beta u_{2,j}^{n+1/2} = u_{1,j}^{n}$$

이고 이고 $u_{0,j}^{n+1/2} = u_{\rm left}$ 이라 하면,
$$ (1+2\alpha_x\beta) u_{1,j}^{n+1/2} - \alpha_x \beta u_{2,j}^{n+1/2} = u_{1,j}^{n} +\alpha_x \beta \cdot u_{\rm left} $$

입니다.

 

마찬가지로,  $u_{I,j}^{n+1/2} = u_{\rm right}$ 이라 하면,
$$ (1+2\alpha_x\beta) u_{I-2,j}^{n+1/2} - \alpha_x \beta u_{I-1,j}^{n+1/2} = u_{I-1,j}^{n} +\alpha_x \beta \cdot u_{\rm right} $$
입니다. 

 

위 예에서는 $u_{\rm left} = u_{\rm right} =0 $인 상황이죠.

 

$y$축 방향으로의 풀이

고정된 $i,~ 0\leq i \leq I$에 대해,

$$ \frac{u_{ij}^{n+1}-u_{ij}^{n+1/2}}{k} = \beta \frac{u_{i,j-1}^{n+1}-2u_{ij}^{n+1}+u_{i,j+1}^{n+1}}{h_y^2} $$  
로 이산화됩니다. 따라서

$$ \frac{u_{ij}^{n+1}-u_{ij}^{n+1/2}}{k} = \alpha_y\beta ( u_{i,j-1}^{n+1}-2u_{ij}^{n+1}+u_{i,j+1}^{n+1} ) $$

즉, $$ -\alpha_y \beta u_{i,j-1}^{n+1} +(1+2\alpha_y\beta) u_{ij}^{n+1} -\alpha_y \beta u_{i,j+1}^{n+1} =  u_{ij}^{n+1/2}$$

입니다. 따라서 $j=1$ 일 때는,

$$ -\alpha_y \beta u_{i,0}^{n+1} + (1+2\alpha_y\beta) u_{i,1}^{n+1} - \alpha_y \beta u_{i,2}^{n+1} = u_{i,1}^{n+1/2}$$

이고 이고 $u_{i,0}^{n+1} = u_{\rm bottom}$ 이라 하면,
$$ (1+2\alpha_y\beta) u_{i,1}^{n+1} - \alpha_y \beta u_{i,2}^{n+1} = u_{i,1}^{n+1/2} +\alpha_y \beta \cdot u_{\rm bottom} $$
입니다. 마찬가지로,  $u_{i,J}^{n+1} = u_{\rm top}$ 이라 하면,
$$ (1+2\alpha_y\beta) u_{i,J-2}^{n+1} - \alpha_y \beta u_{i,J-1}^{n+1} = u_{i,J-1}^{n+1/2} +\alpha_y \beta \cdot u_{\rm top} $$
입니다. 

위 예에서는 $u_{\rm bottom} = 30~,~ u_{\rm top} =100 $인 상황입니다.

 

[중요]
즉 요지는, 함축적 방법(implicit method)에서 tridiagonal matrix를 풀 때, 아래 소개할 code에서 known 벡터에 적절한 경계조건이 반영이 되어야 한다는 뜻입니다.

 

 

 

 

Python Code: OSM으로 예시 풀기

 

이제 python code를 통하여 위 예시를 풀어보겠습니다.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation  # 애니메이션 효과를 주기위한 animation import


def thomas(a, b, c, d):
    """ A is the tridiagnonal coefficient matrix and d is the RHS matrix"""
    """
    a is lower diagonal a2,a3,..,a_N, meaning
    b is diagonal b1,b2,b3,..,b_N meaning
    c is upper diagonal c1,c2,c3,.. c_{N-1} meaning
    """
    N = len(a)
    cp = np.zeros(N, dtype='float64')  # store tranformed c or c'
    dp = np.zeros(N, dtype='float64')  # store transformed d or d'
    X = np.zeros(N, dtype='float64')  # store unknown coefficients

    # Perform Forward Sweep
    # Equation 1 indexed as 0 in python
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    # Equation 2, ..., N (indexed 1 - N-1 in Python)
    for i in np.arange(1, (N), 1):
        dnum = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / dnum
        dp[i] = (d[i] - a[i] * dp[i - 1]) / dnum

    # Perform Back Substitution
    X[(N - 1)] = dp[N - 1]  # Obtain last xn

    for i in np.arange((N - 2), -1, -1):  # use x[i+1] to obtain x[i]
        X[i] = (dp[i]) - (cp[i]) * (X[i + 1])

    return (X)


def HeatEqaution_2D_ImplicitFDM_AnotherExam():
    # u_t(t,x,y) = beta *(u_xx(t,x,y) + u_{yy}(t,x,y)) 0<x,y<1, t>0
    # u(0,x,y) = 0
    # u(t,0,y) = u(t,1,y) =0
    # u(t,x,0)= 30,  u(t,x,1) = 100
    

    beta = 0.01
    t_min, t_max, NSize = 0, 10, 100  # 시점 [0,10]을 100등분
    x_min, x_max, ISize = 0, 1, 100  # x=[0,1]을 100등분
    y_min, y_max, JSize = 0, 1, 100  # y=[0,1]을 100등분

    tSeq = np.linspace(t_min, t_max, NSize + 1)
    xSeq = np.linspace(x_min, x_max, ISize + 1)
    ySeq = np.linspace(y_min, y_max, JSize + 1)

    hx = xSeq[1] - xSeq[0]  # time level 분할 간격
    hy = ySeq[1] - ySeq[0]  # x 분할간격
    k = tSeq[1] - tSeq[0]  # y 분할간격

    alphaX, alphaY = k / hx ** 2, k / hy ** 2  # 연립방정식 풀이에 나오는 α 계산

    xval, yval = np.meshgrid(xSeq, ySeq)

    # u 그리드 판 설계
    u = np.empty((NSize + 1, ISize + 1, JSize + 1))

    # -----initial condition--------------------
    u[0, :, :] = 0

    # -----boundary condition--------------------
    u_left= 0        # x=0 인 부분의 경계조건
    u_right = 0      # x=1 인 부분의 경계조건
    u_top =100       # y=1 인 부분의 경계조건
    u_bottom = 30    # y=0 인 부분의 경계조건

    u[:, 0, :] = u_left  # u(t,x_min,y)=0
    u[:, -1, :] = u_right  # u(t,x_max, y) =0
    u[:, :, 0] = u_bottom  # u(t,x,y_min)=0
    u[:, :, -1] = u_top  # u(t,x,y_max)=0


    # x direction Thomas algorithm matrix
    diagX = np.ones(ISize - 1) * (1 + 2 * alphaX * beta)  # 식 ①로 만든 행렬의 계수 (Tridiagonal matrix)
    underX = np.ones(ISize - 1) * (-alphaX * beta)
    overX = np.ones(ISize - 1) * (-alphaX * beta)

    # y direction Thomas algorithm matrix
    diagY = np.ones(JSize - 1) * (1 + 2 * alphaY * beta)  # 식 ②로 만든 행렬의 계수 (Tridiagonal matrix)
    underY = np.ones(JSize - 1) * (-alphaY * beta)
    overY = np.ones(JSize - 1) * (-alphaY * beta)

    for tval in range(NSize):  # 각 시점에 대해 (t=0,1,..,NSize-1 까지)
        u_temp = u[tval, :, :].copy()  # time level이 고정된 2차원 격자판 설정

        # x축 방향으로 먼저 해결. y값을 고정시키고 x축으로 implicit method 적용
        for j in range(1, JSize):  # 고정된 j에 대해 (j=1,2,..,JSize-1)
            knownX = u[tval, 1:ISize, j].copy()  # Ax =b , A: tridiagonal matrix, b=known
            knownX[-1] += alphaX * beta * u_right   # [중요] 참고
            knownX[0] += alphaX * beta * u_left     # [중요] 참고 u[1:ISize]의 첫번째, 마지막 원소에 경계조건 반영
            u_temp[1:ISize, j] = thomas(underX, diagX, overX, knownX)  # Thomas algorithm
            u_temp[0, j] = u_left  # boundary 값 설정
            u_temp[-1, j] = u_right  # boundary 값 설정

        # 그다음, y축 방향으로 해결, x 값을 고정 후 y축으로 implicit method 적용
        for i in range(1, ISize):
            knownY = u[tval, i, 1:JSize].copy()
            knownY[-1] += alphaY * beta * u_top     # [중요] 참고
            knownY[0] += alphaY * beta * u_bottom   # [중요] 참고
            u_temp[i, 1:JSize] = thomas(underY, diagY, overY, knownY)
            u_temp[i, 0] = u_bottom  # boundray 값 설정
            u_temp[i, -1] = u_top     # boundary 값 설정

        u[tval + 1, :, :] = u_temp.copy()  # 그 다음 time level에 적용함

    def plotheatmap(u_k, k):
        # Clear the current plot figure
        plt.clf()

        plt.title('Temperature at t = {:.2f}'.format(tSeq[k]))
        plt.xlabel("x")
        plt.ylabel("y")

        # This is to plot u_k (u at time-step k)
        plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
        plt.colorbar()

    def animate(k):
        plotheatmap(u[k, :, :].T, k)

    anim = animation.FuncAnimation(plt.figure(), animate, interval=500, frames=NSize, repeat=False)
    anim.save("heat_equation_solution.gif", fps=10)

    plt.show()

if __name__ == '__main__':
    HeatEqaution_2D_ImplicitFDM_AnotherExam()

 

 

애니메이션 부분은  이 블로그를 참고했습니다. 

실행시켜 볼까요?

 

 

 

 

시간이 흐를수록 위의 뜨거운 열이 아래로 전달되는 과정을 볼 수 있습니다.

 

그럴싸하네요.