6. 수치적 근사법으로 해 구하기-오일러의 방법 (Numerical Approximation-Euler's Method)

 

지금껏 일계 미분방정식의 해를 구하는 여러가지 방법을 알아보았다.

선형인 경우엔 적분 인자를 곱하여 풀었고

변수분리형인 경우엔 변수를 분리한 후 적분하여 풀었고

완전 미분방정식인 경우에는 원함수가 무엇인지 적분으로 추정해내서 풀었다.

 

그리고 여기서 다음과 같은 일계 미분방정식은 특정 조건을 만족하면

항상 유일한 해가 존재한다는 정리를 소개했었다.

$$ \dfrac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 $$

 

하지만 위에 나열하였듯 풀이 방법이 존재하는 경우는 아주 특수한 경우인 것이고

대부분의 경우에는 풀이법이 알려져있지 않다.

그럼에도 해가 존재한다는게 증명이 되어있으니 수치적 근사하는 방법으로라도

해를 찾아볼 수는 있을 것이다.

이번 글에서 수치적 근사 방법 중 하나인 오일러의 방법을 알아볼 것이다.

 

 

 

---

 

 

 

다음과 같은 미분 방정식을 보자.

$$ \dfrac{dy}{dt} = 3 - 2t - 0.5y $$

 

이 미분방정식은 사실 선형 일계 미분방정식이므로 적분인자를 통해 해를 구할 수 있다.

하지만 그러한 방법을 모른다고 가정하고 근사 해를 구하려면 어떻게 해야하는지 살펴보자.

이 미분방정식의 방향장(Direction Field) 를 그려보면 다음과 같이 나타난다.

 

해가 어떻게 그려지는지 모양을 보여주기 위해

초기값 $(2, -1)$ 을 지나는 해를 초록색으로 그려넣었다.

 

보다시피 해가 방향장의 선이 표현되는 대로 쭉 이은 것 처럼 생겼다.

여기서 해 곡선을 방향장에 표현된 기울기를 갖는 선분 조각들을 이으면

실제 해랑 비슷한 모양이 나온다는 아이디어를 얻을 수 있다.

이 방법이 오일러가 1768년 쯤 만들어낸 방법이라고 한다.

 

이제 그림을 자세히 들여다 보자.

선분 조각 하나를 만들고자 한다. $(2, -1)$ 부분에서 시작해보자.

$(2, -1)$ 에서 찾고자 하는 해 $y$ 의 기울기 $y'$ 는

주어진 미분방정식에 값을 대입하여 얻을 수 있다. 이 때 $y'(2) = -0.5$ 이므로

$(2, -1)$ 을 지나고 기울기가 $-0.5$ 인 직선 $ y = -0.5(t - 2) - 1 $ 으로 선분 조각을 시작하자. 

 

 

편의상 진행 step의 크기를 1로 잡아서 $2 \le t \le 3$ 에서 끊어냈다.

다음 선분 조각은 이 선분의 끝부분인 $(3, -1.5)$ 에서 시작한다.

이 점에서의 기울기는 역시 주어진 미분방정식에 점을 대입하여 구할 수 있다.

$ y'(3) = -2.25 $ 이므로, 다음 선분 조각은 직선 $y = -2.25(t-3)-1.5$ 에서 얻어낸다.

 

아래 그림은 두 단계를 진행한 근사 해와 실제 해이다.

 

 

실제 해랑 선분 조각을 이어 만든 근사 해랑 오차가 좀 나보이는 모습이다.

하지만 단계의 크기를 $1$ 보다 작게 한다면,

훨씬 실제 해와 근접한 모습을 얻어낼 수 있을 것이라 예상할 수 있다.

아래의 사진은 단계의 크기를 $0.5$ 로 잡고 만들어낸 근사 해의 모습이다.

 

 

 

---

 

 

 

이제 이 방법을 일반화해서 표현해보자.

시작 점을 $(t_0, y_0)$ 이라고 하면 첫 선분을 만드는 직선의 방정식은 다음과 같다.

$$ y = y_0 + f(t_0, y_0)(t - t_0) $$

 

 

첫 단계 진행 후 도달하는 점을 $(t_1, y_1)$ 이라고 하고

두 단계 진행 후 도달하는 점을 $(t_2, y_2)$,

$n$ 단계 진행 후 도달하는 점을 $(t_n, y_n)$ 이라고 하면 다음과 같은 점화식을 세울 수 있다.

$$ y_{n+1} = y_n + f(t_n, y_n)(t_{n+1} - t_n), \; n = 0, 1, 2, ... $$

 

간단히 표현하기 위해 $f(t_n, y_n) = f_n$ 이라고 하면 다음과 같이 표현된다.

$$ y_{n+1} = y_n + f_n \cdot (t_{n+1} - t_n) $$

혹은 더욱 간단히 하기 위해 진행 단계의 크기 $(t_{n+1} - t_n) = h$ 라고 정의하면

$$ y_{n+1} = y_n + f_n \cdot h $$ 

 

한글로 풀어 쓴다면

다음 단계의 $y$ 값 $=$ 현재 단계의 $y$ 값 + 현재 점에서 기울기 $\times$ 단계의 크기

라고 할 수 있겠다.

 

 

 

 

---

 

 

 

이 문제는 원서 11판 기준 연습문제에 15번으로 수록된 문제이다.

 

오일러 방법의 수렴성 : 
미분방정식 $y' = f(t, y), \; y(t_0) = y_0 $ 의 오일러 방법으로 수치적 근사시킨 해가
특정 조건에서 단계의 크기인 $h$ 를 $0$ 으로 보냄으로써 실제 해랑 정확히 같도록 수렴하게 된다.
다음 예에서 그 방법을 보여줄 것이다.
$$ y' = 1 - t + y, \quad y(t_0) = y_0 $$

$ a. $ 이 미분방정식의 실제 해가 $y = \phi(t) = (y_0 - t_0)e^{t - t_0} + t$ 가 됨을 보여라.

$ b. $ 오일러의 방법을 이용하여 다음을 보여라.
$$ y_k = (1 + h)y_{k-1} + h - ht_{k-1}, \quad k = 1, 2, 3, ... $$

$ c. $ $y_1 = (1 + h)(y_0 - t_0) + t_1$ 임을 참고하여 정수 $n$ 에 대해 다음이 성립함을 귀납법으로 보여라.
$$ y_n = (1 + h)^n (y_0 - t_0) + t_n $$

$d. $ $t > t_0$ 을 만족하는 고정점 $t$ 를 생각하자. 주어진 $n$ 에 대해 $h$ 를 $h = \frac{t - t_0}{n}$ 으로 선택하자.
그러면 모든 $n$ 에 대해 $ t_n = t $ 를 만족함을 알 수 있고, $n \to \infty$ 이면 $h \to 0$ 임도 알 수 있다.
이러한 $h$ 를 문제 $c$ 의 결과에 대입하고 $n \to \infty$ 하여 $y_n \to \phi(t)$ 가 됨을 보여라.
힌트 : $\displaystyle \lim_{n \to \infty} (1 + \frac{a}{n})^n = e^a$

$y' = 1 - t + y$ 를 표준형으로 정리하면 $y' - y = 1 - t$ 이므로 일계 선형 미분방정식이다.

따라서 적분인자 방법으로 해를 구할 수 있다.

즉, $\mu(t)$ 가 존재하여 주어진 미분방정식은 $ [\mu(t)y(t)]' = \mu(t) \cdot (1 - t)$ 가 성립한다.

 

 

$$ \mu = e^{\int -1 dt} = e^{-t} $$

이므로 주어진 미분방정식은 다음과 같다.

$$ [e^{-t} y]' = e^{-t}(1 - t) $$

양변을 $t$ 에 대해 적분하면

$$ e^{-t} y = \int e^{-t}(1 - t) dt = te^{-t} + C $$

즉 일반 해는 다음과 같다.

$$ y(t) = \phi(t) = t + Ce^t $$

 

 

초기값 $y(t_0) = y_0$ 을 대입하면

$$ y_0 = t_0 + Ce^{t_0} $$

정리하여 $C$ 를 구하면

$$ C = \dfrac{y_0 - t_0}{e^t_0} $$

 

 

구한 해에 대입하면 다음과 같다.

$$ y = \phi(t) = t + (y_0 - t_0)e^{(t - t_0)} $$

오일러의 방법 공식은 다음과 같다고 하였다.

$$ y_{n+1} = y_n + f_n \cdot h $$

 

$f_n = f(t_n, y_n) = 1 - t_n + y_n$ 이므로 $y_{n+1}$ 은 

$$ \begin{align} y_{n+1} &= y_n + (1 - t_n + y_n)h \\ &= (1 + h)y_n + h - ht_n, \quad (n = 0,1,2...)  \end{align} $$

 

$n + 1$ 대신 $k, \quad (k = 1,2,3)$ 을 대입하면

$$ y_k = (1+h)y_{k-1} + h - ht_{k-1} $$ 

우선 왜 $y_1 = (1+h)(y_0 - t_0) + t_1$ 로 주어졌는지 설명하자면

$b.$ 번 문제에서 $y_k = (1+h)y_{k-1} + h - ht_{k_1}$ 이라고 구했었다.

 

여기에 $ k = 1 $ 을 대입하면 다음과 같다.

$$ y_1 = (1+h)y_0 + h - ht_0 $$

$(1+h)$ 에 $t_0$ 을 묶기 위해 다음과 같이 $t_0$ 을 한 번 빼준 후 $-(1+h)t_0$ 로 묶어주고 $t_0$ 을 다시 더해주자.

$$ y_1 = (1+h)y_0 + h - (1+h)t_0 + t_0 $$

$$ \Longrightarrow y_1 = (1+h)(y_0 - t_0) + t_0 + h $$

$t_1 = t_0 + h$ 이므로

$$ y_1 = (1+h)(y_0 - t_0) + t_1 $$

 

 

바로 위에서 $n = 1$ 인 경우는 성립함을 보였다.

이제 $n$ 번째 항에서 성립한다고 가정할 때 $n+1$ 에서도 성립하는지 보이자.

 

$b$ 결과에 의해 $y_{n+1} = (1 + h)y_n + h - ht_n$ 이다.

한편 $n$ 일 때 성립한다고 가정했으므로 $y_n = (1+h)^n (y_0 - t_0) + t_n$ 이고 이를 대입하면

$$ \begin{align} y_{n+1} &= (1+h){ (1+h)^n (y_0 - t_0) + t_n } + h - ht_n \\ &= (1+h)^{n+1}(y_0 - t_0) + (1+h)t_n + h - ht_n \\ &= (1+h)^{n+1} (y_0 - t_0) + t_n + h \end{align} $$

$t_{n+1} = t_n + h$ 이므로

$$ y_{n+1} = (1+h)^{n+1}(y_0 - t_0) + t_{n+1} $$

 

이렇게 $n$ 일 때 성립한다고 가정했을 때 $n+1$ 일 때도 성립함을 증명하였다.

따라서 수학적 귀납법에 의해 명제는 참이다.

문제가 무엇을 말하는지부터 해석해보자.

주어진 initial point $t_0$ 보다 좌표축 상 오른쪽에 해당하는 고정점 $t$ 를 선택했다.

왼쪽에 해당하는 고정점의 경우에 대해서는 유사하게 증명할 수 있기 때문에 오른쪽만 선택해서

증명해 보이는 것으로 예상할 수 있다.

 

또 $h$ 를 살펴보면, $t_0$ 와 $t$ 사이의 거리를 $n$ 등분한 크기로 잡았음을 볼 수 있다.

$t_{n+1} = t_n + h$ 임을 생각해보면 $t_0$ 에서 $n$ 단계 진행하면 $t_0 + nh$ 인데

주어진 $h$ 를 대입해보면 $t_0 + t - t_0 = t$ 이므로 $t_n = t$ 가 되는 것은 자명하다. 

$n \to \infty$ 면 $h \to 0$ 인것도 자명하다.

 

문제에서 마지막에 보여라고 한 것인 $y_n(t) \to \phi(t)$ 이 것은 한글로 풀어서 해석하자면,

오일러 방법으로 근사시킨 $y_n(t)$ 이 임의로 택한 $t$ 값에서 실제 해인 $\phi(t)$ 랑

같은 값을 나타내도록 수렴한다는 말이다.

이게 된다는 것을 보이면 실제 해를 구하는 방법이 밝혀지지 않은 경우에도 오일러 방법으로 근사시킨 후

$n \to \infty$ 를 하여 실제 해의 값을 구해낼 수 있게 된다.

 

 

 

풀이

 

$$ y_n = (1 + h)^n (y_0 - t_0) + t_n $$

위 식에 $h = (t - t_0) / n$ 를 대입하면 다음과 같다.

$$ y_n = \left( 1 + \dfrac{t - t_0}{n} \right)^n (y_0 - t_0) + t $$

$h$ 를 주어진 바와 같이 잡으면 $t_n = t$ 라고 하였으므로 $t_n$ 을 $t$ 로 바꿔주었다.

 

힌트의 식에 의해 $n \to \infty$ 극한에서는 다음과 같다.

$$ \lim_{n \to \infty} y_n = e^{t - t_0} (y_0 - t_0) + t $$

우변은 문제 $a$ 에서 보였던 실제 해랑 같음을 볼 수 있다.

 

따라서 $n \to \infty$ 이면 $y_n(t) = \phi(t)$ 이다.

참고로 힌트의 식이 왜 저렇게 나왔는지가 궁금하다면

$\left( 1 + \frac{a}{n} \right)^n$ 에 자연로그를 씌운 식을 $g(n)$ 이라 한 후

$\displaystyle \lim_{n \to \infty} g(n)$ 을 로피탈 법칙을 적용하여 구해보면 쉽게 증명할 수 있다.