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

 

 

이 글은

 

2023.02.06 - [수학의 재미/아름다운 이론] - 2차원 Heat Equation의 풀이 #1 : Operator Splitting Method(OSM)

2차원 Heat Equation의 풀이 #1 : Operator Splitting Method(OSM)

에서 이어집니다.

$$\frac{\partial u}{\partial t}(t,x) = \frac{\partial^2 u}{\partial x^2}(t,x)  $$

형태의 2차원 Heat Equation을 해결하기 위해서 OSM(Operator Splitting Method)라는 기법을 사용한다고 하였습니다.

OSM 을 잠시 복습해 볼까요?

 

OSM

시점과 $t$와 2차원좌표 $(x,y)$를 각각

 

$$ 0=t_0<t_1<\cdots <t_{N-1}<t_N =T, $$

$$ x_{\min} = x_0< x_1 <\cdots < x_{I-1} <x_I =x_{\max}, $$
$$ y_{\min} = y_0< y_1 <\cdots < y_{J-1} <y_J =y_{\max}, $$

라 합니다(각각의 분할 간격은 $k, h_x,h_y$입니다.)

또 2차원 편미분 방정식이 만족하는 해 $u$에 대해 시점 레벨 $n$ 즉 $t=t_n$에서  $(x_i, y_j)$ 에서의 $u$의 값을
$$ u_{ij}^n : = u(t_n, x_i,y_j)$$
라 표현하고 시점 레벨 $n$과 $n+1$ 사이에 가상의 시점 레벨 $n+\frac12$를 두어, 아래처럼 변수 별로 FDM을 진행합니다.

 

 

 

$n$부터 $n+1$ time level까지의 시간을 반으로 나누어 $x$ 변수 먼저 해결하고, 그다음 $y$변수를 해결합니다. 

다음과 같이요.

① $n$과 $n+\textstyle{\frac12}$ 사이에서
$$\frac{u_{ij}^{n+\frac12} -u_{ij}^n}{k} = \frac{u_{i+1,j}^{n+\frac12}-2u_{ij}^{n+\frac12}+u_{i-1,j}^{n+\frac12}}{h_x^2}$$
를 푼다.

② $n+\textstyle{\frac12}$ 과 $n+1$ 사이에서
$$\frac{u_{ij}^{n+1} -u_{ij}^{n+\frac12}}{k} = \frac{u_{i,j+1}^{n+1}-2u_{ij}^{n+1}+u_{i,j-1}^{n+1}}{h_y^2}$$
을 풀어 완성한다.

 

이를 정리하면, 다음과 같습니다.

$\alpha_x= k/h_x^2~,~\alpha_y=k/h_y^2$이라 하면,

① $ -\alpha_x u_{i-1,j}^{n+\frac12} +(1+2\alpha_x) u_{i,j}^{n+\frac12} -\alpha_x u_{i+1,j}^{n+\frac12} = u_{ij}^n$
   
즉, 고정된 $j$에 대해 위 식을 풀어 $u_{ 0:I, j}^n $에서 $u_{0:I, j}^{n+\frac12}$를 얻는다. 이를 모든  $0\leq j\leq J$에 대해 반복한다.

② $ -\alpha_y u_{i,j-1}^{n+1} +(1+2\alpha_y) u_{i,j}^{n+1} -\alpha_y u_{i,j+1}^{n+1} = u_{ij}^{n+\frac12}$

즉 고정된 $i$에 대해 위 식을 풀어 $u_{i, 0:J}^{n+\frac12}$에서 $u_{i,0:J}^{n+1}$을 얻는다. 이를 모든 $0\leq i\leq I$에 대해 반복한다.

 

2차원 Heat Equation 예제

 

다음의 열방정식이 있다고 해 봅시다. 자연수 $m, n$에 대해,

$$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} ~~,~~ 0\leq t\leq 1,  0\leq x \leq 1 , 0\leq y \leq 1$$ ○ 경계조건 : $ u(t, x, 0) =0, u(t,x,1) =0, u(t,0,y) =0, u(t,1,y)=0$ ○ 초기조건:  $u(0,x,y) = \sin(m\pi x) \sin(n\pi y) $

Exact Solution : $$ u(t,x,y) = \sin(m\pi x) \sin(n\pi y) e^{-(m^2+n^2)\pi t} $$

위 방정식의 정확한 해(Exact Solution)를 $t=0$에서 $t=1$까지 그려보면 다음과 같습니다($m, n=1$일 때)

 

( 위 그림을 그린 code는 아래 더 보기 글을 참고해 보시기 바랍니다.)

 

def Draw_2DHeatEquation():
    t_min = 0
    t_max = 1
    NSize = 50

    x_min = 0
    x_max = 1
    ISize = 100

    y_min = 0
    y_max = 1
    JSize = 100

    m = 1
    n = 1
    pi = np.pi

    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)

    xval, yval = np.meshgrid(xSeq, ySeq)
    closed_solution = lambda t, x, y: np.sin(m * pi * x) * np.sin(n * pi * y) * np.exp(-(m ** 2 + n ** 2) * pi ** 2 * t)
    zarray = np.zeros((NSize + 1, ISize + 1, JSize + 1))



    for i in range(NSize + 1):
        zarray[i, :, :] = closed_solution(tSeq[i], xval, yval)
        # zarray[:, :, i] = f(x, y, 1.5 + np.sin(i * 2 * np.pi / frn))

    def update_plot(frame_number, zarray, plot):
        plot[0].remove()
        plot[0] = ax.plot_surface(xval, yval, zarray[frame_number, :, :], cmap="hot")

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    plot = [ax.plot_surface(xval, yval, zarray[0, :, :], color='0.75')]
    ax.set_zlim(0, 1.1)
    ani = animation.FuncAnimation(fig, update_plot, NSize + 1, fargs=(zarray, plot))
    ani.save('test.gif', fps=5)
    plt.show()

 

 

 

이제 예제의 열방정식을 OSM 방법으로 풀어보겠습니다.

 

 

 

 

Python Code: OSM으로 2D Heat Equation 풀기

 

def thomas(a, b, c, d):
    # 생략 (다른 글 참조)

def HeatEqaution_2D_ImplicitFDM():
    # u_t(t,x,y) = u_xx(t,x,y) + u_{yy}(t,x,y) 0<x,y<1, t>0
    # u(0,x,y) = sin(mπx)sin(nπy)
    # u(t,0,y) = u(t,1,y) u(t,x,0)=u(t,x,1) = 0
    # closed form solution : u(t,x,y)= sin(mπx)sin(nπx)exp(-(m^2+n^2)π^2 t)

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

    m, n = 1, 1
    pi = np.pi      # pi 값 설정

    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)
    # exact solution
    closed_solution = lambda t, x, y: np.sin(m * pi * x) * np.sin(n * pi * y) * np.exp(-(m ** 2 + n ** 2) * pi ** 2 * t)
    
    # u 그리드 판 설계
    u = np.empty((NSize + 1, ISize + 1, JSize + 1))

    # -----boundary condition--------------------
    u[:, 0, :] = 0  # u(t,x_min,y)=0
    u[:, -1, :] = 0  # u(t,x_max, y) =0
    u[:, :, 0] = 0  # u(t,x,y_min)=0
    u[:, :, -1] = 0  # u(t,x,y_max)=0

    # -----initial condition--------------------
    u[0, :, :] = closed_solution(t_min, xval, yval).transpose()
    
    # u[0,:,:]은 vertical 방향이 x, horizontal 방향이 y이고
    # xval, yval은 2차원 좌표축처럼 vertical 방향이 y축, horizontal이 x축
    
    # x direction Thomas algorithm matrix
    diagX = np.ones(ISize - 1) * (1 + 2 * alphaX) # 식 ①로 만든 행렬의 계수 (Tridiagonal matrix)
    underX = np.ones(ISize - 1) * (-alphaX)
    overX = np.ones(ISize - 1) * (-alphaX)

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

    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]  # Ax =b , A: tridiagonal matrix, b=known
            u_temp[1:ISize, j] = thomas(underX, diagX, overX, knownX)  #Thomas algorithm
            u_temp[0, j] = 0              # boundary 값 설정
            u_temp[-1, j] = 0             # boundary 값 설정
        
        # 그다음, y축 방향으로 해결, x 값을 고정 후 y축으로 implicit method 적용
        for i in range(1, ISize):
            knownY = u[tval, i, 1:JSize]
            u_temp[i, 1:JSize] = thomas(underY, diagY, overY, knownY)
            u_temp[i, 0] = 0              # boundray 값 설정
            u_temp[i, -1] = 0             # boundary 값 설정

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

    
    # exact solution과 FDM 값이 얼마나 차이나는지 알아볼 예정
    errorVec = []
    for tIdx in range(NSize + 1):           # time level을 변화하며
        cs_res = closed_solution(tSeq[tIdx], xval, yval)  # 고정된 time에 대한 exact solution
        fdm_res = u[tIdx, :, :].T                         # 위 고정된 time에서의 FDM값
        error = np.linalg.norm(cs_res - fdm_res) / (ISize + 1) / (JSize + 1)
        # 두 array의 차이의 Frobenius norm을 계산하여 행렬 size로 나눠줌
        # 그러면 각 (x,y)좌표에서 excat solution과 FDM사이의 평균 오차로 볼 수 있음.
        errorVec.append(error)      # 위에서 구한 error을 time별로 적립
        

    print('Mean Error : {:.4f}'.format(np.mean(errorVec)))   # 위에서 구한 error의 mean 표시

    fig = plt.figure()        # Graph를 위한 fig(캔버스) 설정
    ax = fig.add_subplot(projection='3d')  # 그 fig라는 캔버스에 그래프 공간 ax 설정, 3D graph

    def update_plot(idx):
        ax.clear()          # 그림 그리는 공간 ax 지움
        ax.plot_surface(xval, yval, u[idx, :, :].T, cmap='hot') # 3차원 그래프 plotting, 
        ax.set_zlim(0, 1.1) # z축 값 변화 고정

    ani = animation.FuncAnimation(fig, update_plot, NSize + 1)
    ani.save('test2d.gif', fps=24)  # 1초당 frame 24로 고정
    plt.show()

 

 

1D 함축적 방법에서 다루었던 스킬을 썼습니다. 위의 식 ①, ② 를 Thomas Algorithm으로 푸는 것이었죠. 

1D 함축적 방법은 

2022.08.01 - [수학의 재미/아름다운 이론] - FDM #7, Heat Equation의 풀이(3)

FDM #7, Heat Equation의 풀이(3)

을 참고해 보시기 바라고,  3중 대각행렬을 푸는 Thomas Algorithm은

 

2022.07.29 - [수학의 재미/행렬 이론] - 3중 대각행렬의 풀이

3중 대각행렬의 풀이

을 참고하시기 바랍니다(위의 코드 생략된 부분을 찾을 수 있습니다.)

 

 

 

결과를 보실까요?

 

Mean Error : 0.0003
# exact solution과 FDM 결과값의 평균 error가 아주 준수합니다. 거의 붙어있죠

 

위의 Exact Solution 형태와 비슷한 움직임

 

 

 

어떤가요? 중심부가 달궈진 프라이팬을 실온에 놔뒀을 때, 열이 식는 모습이 보이는지요.

 

다음 글에서는 다른 예제로 재미있는 열 방정식을 풀어보겠습니다.