미분방정식을 연립방정식으로, FDM #2

이번 글은

2022.07.22 - [수학의 재미/아름다운 이론] - 미분방정식을 연립방정식으로, FDM #1

미분방정식을 연립방정식으로, FDM #1

에서 이어집니다. 위 글에서 우리는 함수의 1계 미분과 2계 미분을 차분 형태(difference form)로 바꾸는 걸 배웠습니다. 그럼 실제로 FDM으로 어떻게 미분방정식을 해결하는지 알아보도록 하죠.

 

 

일반적인 형태의 1계 미분방정식을 풀어봅시다.

 

적분인자를 이용한 풀이

 

미분방정식

$$ f'(x) + p f(x) = q ~~,~~ f(0)=r, $$

($p, q, r$은 상수)를 어떻게 풀까요? 양변에 적분 인자

$$ \exp(px)$$

를 곱해주면,

$$ e^{px} f'(x) + pe^{px} f(x) = q e^{px} $$ 

이고

$$ \left(e^{px} f(x) \right)' = q e^{px}$$

가 됩니다. 따라서 적분식으로 쓰면

$$ \int_0^x \left(e^{pt} f(t) \right)' = \int_0^x q e^{pt} dt $$

이고

$$ e^{px} f(x) - f(0)  = \frac{q}{p} \left( e^{px}-1  \right) $$

따라서 $f(0)=r$ 을 이용하여 정리하면,

$$ f(x) = \frac qp + \left(r-\frac qp\right) e^{-px} $$

 

 

FDM을 이용한 풀이

미분방정식

$$ f'(x) + p f(x) = q ~~,~~ f(0)=r \tag{1}$$

을  $f'(x)$에 대한 forward difference와 backward differece 각각을 사용하여 풀어보겠습니다.

 

풀이를 위해 $x$의 최댓값을 $x_{\max}$라 합시다. 그리고 $[0, x_{\max}]$을 $N$등분을 하여

$$ 0=x_0 < x_1 <x_2 <\cdots < x_{N-1} < x_N = x_{\max}$$

로 나눕니다. 총 $N+1$개의 $x$ 수열이 생기겠죠.

 

$x_{i+1}$과 $x_i$사이의 간격을 $h$라 합시다. 그리고 

$$f(x_i)= u_i$$

라 정의합니다(FDM을 풀 땐 $u$라는 기호를 많이 쓰더군요) 

 

Forward Difference

$x_i$에서 $f'$을 forward difference로 구하면

$$ f'(x_i ) = \frac{f(x_i +h) -f(x_i)}h = \frac{f(x_{i+1})-f(x_i)}h = \frac{u_{i+1}-u_{i}}h $$

입니다.

따라서 식(1) 은

$$ \frac{u_{i+1}-u_{i}}h + p u_i = q , x_0 =r\tag{2}$$

이라는 점화식으로 바뀌게 됩니다. 보기 좋게 정리하면

$$ u_{i+1} +(ph-1)u_i = qh ~,~ x_0 = r\tag{3}$$

입니다.

 

Backward Difference

이번에는 $x_{i+1}$에서의 $f'$의 backward difference를 사용하여 구해볼까요?

 

$$ f'(x_{i+1}) = \frac{f(x_{i+1}-h)-f(x_{i+1})}{-h} = \frac{u_{i+1}-u_{i}}h $$

이고 식(1)은

$$ \frac{u_{i+1}-u_{i}}h + p u_{i+1}= q , x_0 =r$$

이 됩니다. 식(2)과 비교하면 $p$를 계수로 갖는 항이 $u_i$에서 $u_{i+1}$로 바뀌었죠. 더 정리하면

$$ u_{i+1} = \frac1{1+ph} u_n + \frac{qh}{1+ph} \tag{4}$$

입니다.

 

식(3), (4) 모두 점화식입니다. 이 점화식을 풀면 $u_0, u_1,\cdots, u_N$을 모두 얻을 수 있습니다.

이렇게 연속적인 $x$를 잘게 나누어 수열 $\{x_i\}$로 만든 뒤 각각의 점에서 함숫값 $u_i$를 찾아내는 것이 바로 FDM 전법입니다.

 

 

Python으로 풀어보기

import numpy as np
import matplotlib.pyplot as plt

def first_ode_test():
    # Solve:
    # f'(x) + p f(x) = q, f(0)= r 
    
    p, q, r = 3, 2, 1
    x_min, x_max = 0, 10
    Nsteps = 30
    x = np.linspace(x_min, x_max, Nsteps + 1)
    sol = q / p + (r - q / p) * np.exp(-p * x)

    h = x[1] - x[0]
    u = np.zeros(Nsteps + 1)
    v = np.zeros(Nsteps + 1)
    u[0] = r
    v[0] = r
    for i in range(Nsteps):
        u[i + 1] = q * h - (p * h - 1) * u[i]
        v[i + 1] = 1 / (1 + p * h) * (v[i] + q * h)
    plt.plot(x, sol, label='exact solution')
    plt.plot(x, u, label='forward FDM')
    plt.plot(x, v, label='backward FDM')
    plt.legend()
    plt.show()

    print('maximum difference between exact sol and forward  FDM: {:.3f} '.format(np.max(abs(sol - u))))
    print('maximum difference between exact sol and backward FDM: {:.3f} '.format(np.max(abs(sol - v))))

 

코드를 간략히 살펴봅시다.

 

    p, q, r = 3, 2, 1		# 미분방정식 세팅
    x_min, x_max = 0, 10	# x의 초기값은 0이고 최댓값을 10으로 설정
    Nsteps = 30				# x_min과 x_max 사이를 균등격자로 나누는 갯수
    x = np.linspace(x_min, x_max, Nsteps + 1)	#수열 {x_i} 생성

 

---

    sol = q / p + (r - q / p) * np.exp(-p * x)	#exact solution, 적분인자 풀이법 참고

 

---

    h = x[1] - x[0]		# 수열사이의 간격 -> FDM difference 를 결정하는 수			
    u = np.zeros(Nsteps + 1)	# forward difference로 구한 함숫값
    v = np.zeros(Nsteps + 1)	# backward difference로 구한 함숫값
    u[0] = r	# 초기값 세팅
    v[0] = r	# 초기값 세팅
    for i in range(Nsteps):
        u[i + 1] = q * h - (p * h - 1) * u[i]		# forward difference의 점화식
        v[i + 1] = 1 / (1 + p * h) * (v[i] + q * h)	# backward difference의 점화식

 

 

결과를 보겠습니다.

 

maximum difference between exact sol and forward  FDM: 0.123 
maximum difference between exact sol and backward FDM: 0.044 

Process finished with exit code 0

○ forward 방법은 아예 비슷하지도 않습니다. 반면 backward 방법은 비슷은 하지만 이격은 약간 있죠(최대 거리 0.044)

○ 아무래도 Nstep을 너무 조금 줘서 $h$가 크기 때문일 것입니다. 

 

그러면 $h$가 더 커지면 결과가 더 이상할까요? Nstep에 10을 줘보겠습니다.

 

 

maximum difference between exact sol and forward  FDM: 341.333 
maximum difference between exact sol and backward FDM: 0.067 

Process finished with exit code 0

○ forward 방법은 심하게 요동치네요. 실제 해와 backward 방법이 0으로 보일 정도로 요동을 쳐버립니다.

○ 반면에 backward 방법은 그림 속엔 보이지 않지만, maximum 차이가 0.067로서 어느 정도 만족스럽습니다.

 

 

 

그렇다면 Nstep을 확 쪼개 보겠습니다. 10000 정도 줘보겠습니다.

 

maximum difference between exact sol and forward  FDM: 0.000 
maximum difference between exact sol and backward FDM: 0.000 

Process finished with exit code 0

○ 이제 속이 시원할 정도로 똑같이 나오는군요.

○ Step을 아주 잘게 쪼개니 forward나 backward나 해를 잘 설명하는 듯합니다.

 

 

 

결과 분석

이 밖에도 Nstep을 바꾸어 가면서 했더니 잘게 쪼개면 forward 방법이나 backward 방법 간 차이가 없습니다. 

그런데 아무래도 잘게 쪼개면 계산량이 많아지겠죠. 그래서 조금 성기게 잘랐더니

backward 방법은 그래도 해를 잘 따라가는데 forward는 들쑥날쑥 변동이 심합니다.

미리 언급을 하자면 각각의 명칭이 아래와 같이 다릅니다.

 

차분방식	정식 이름
forward difference	explicit method (명시적 방법)
backward difference	implicit method (함축적 방법)

 

두 방법 중 어떤 방법이 계산이 편할까요? 또는 어느 방법이 안정적인 해를 얻을 수 있을까요?

다음 글에서 계속 이어나가겠습니다.