수치 미분 방정식 풀이법을 이용하여 최종값 조건이 주어진 문제 풀기 : 슈팅 방법

지난 두 번의 포스팅을 통해 수치해석적인 방법을 통하여 미분 방정식을 푸는 방법을 소개하였습니다. 오일러 방법과 룽게-쿠타 방법은 연속적인 값을 갖는 시간을 불연속화 하여 $t_n = t_0 + n \cdot h$와 같이 표현하고, 이 때의 함수값이나 미분값 $y(t_n), y'(t_n), y''(n), ...$을 이용하여 수열 $\{y_{n}\}$이 만족해야 하는 점화식을 구한다는 공통점이 있었습니다. 점화식을 컴퓨터를 이용하여 반복 계산하게 되는데, 초기 상태 $t_0$에서의 초기값 $y(t_0)= y_0, y'(t_0) = y'_0 , ... $ 이 주어진다면, 이로 부터 $y_1, y'_1, y''_1, ...$을 구할 수 있고, 이 방법을 계속해서 진행하면 일반적인 $y_n$에 대한 값을 얻을 수 있습니다. 컴퓨터 프로그래밍만 성공적으로 한다면, 초기 조건이 주어진 임의의 미분 방정식을 풀 수 있습니다. (사실, "임의의" 미분 방정식은 아니고, 지난 포스팅에서 알아 본 것들과 같이 수치적으로 풀이가 가능한 형태의 미분 방정식이라고 해야 맞습니다)

 

초기값 문제와 최종값 문제 : 최종값 문제는 초기값 문제로 바꿔서 풀기 = 슈팅 방법 (shooting method)

 

지난 두 번의 포스팅에서 다룬 문제에서 미분 방정식의 경계 조건은 초기값에서의 함수값 혹은 함수의 미분값 이었습니다. 즉, 정의역 $t \in [t_0, t_N]$ 일때, 함수 $y(t)$에 대한 2계 미분 방정식

$$y''(t) = f(t, y, y')$$

가 주어졌을 때, 경계 조건은 $y(t_0) = y_0, y'(t_0) = y'_0$와 같이 $t=t_0$에서의 함수값과 미분값이었습니다. 고전 역학에서의 문제는 대부분 현재(즉, $t=t_0$)의 상태(위치와 속도)를 알고 있을 때, 미래($t \gt t_0)$에서의 상태를 구하는 것이 목적이기 때문에, 초기 조건이 주어지는 것이 매우 자연스럽습니다. 오일러 방법이나 룽게-쿠타 방법을 이용하여 수치적인 해를 구할 때에도 $y_0, y'_0$을 알고 있는 상황에서 $y_n = g(t_{n-1}, y_{n-1}, y'_{n-1}), y'_n = h(t_{n-1}, y_{n-1}, y'_{n-1})$ 와 같이 $n$이 증가하는 방향으로 함수값을 얻을 수 있는 점화식을 이용하였습니다. 

 

그러나, 일반적인 문제에서의 경계 조건은 초기값 문제가 아닐 수도 있습니다. 예를 든다면, $t=t_0$에서의 $y_0, y'_0$이 주어지는 대신, $t_0, t_N$에서의 함수값 $y_0, y_N$가 주어지는 경우를 생각할 수 있습니다. 시간 $t = t_N$은 정의역의 최종값이기 때문에, (이 포스팅에서는) 이와 같은 경계조건을 "최종값 문제"라고 부루도록 하겠습니다. 최종값 문제는 "최종 상태에서 원하는 결론을 얻기 위한 방정식"이라고 볼 수 있습니다. 

 

미분 방정식의 문제를 고전 역학의 문제로 생각하여 이해를 쉽게 해 보도록 하겠습니다. (사실, 글 이곳 저곳에서 미분 방정식을 푸는 수학적인 문제와 고전 역학 문제를 푸는 물리 문제를 섞어어 이야기 하곤 했습니다. 이처럼 고전 역학, 더 넓게는 물리학과 미분 방정식을 관련이 깊습니다) 최종값 문제는 "특정한 포텐셜 필드 하에서 입자의 운동의 초기 위치와 최종 위치가 주어졌을 때, 이를 만족하는 입자의 궤적을 구하는 문제"라고 할 수 있습니다. 예를들어서, 농구를 할 때 공을 던져서 골대에 넣으려고 한다면, $t_0$은 내가 공을 던지는 시점 $y_0$은 공을 던지는 시점에서의 공의 위치, $t_N$은 공이 골대에 들어가는 순간의 시점, $y_N$은 골이 골대에 들어가는 시점의 공의 위치(즉, 골대 안)이 됩니다. $t_N$은 손에서 떨어진 공이 골대에 들어가기 까지 걸리는 시간인데, 현실 세계에서는 약 $1~3$초 정도 걸리는 것 같습니다. 공을 던지는 초기 속도에 따라서 $t_N$초 후의 공의 위치 $y(t_N)$이 결정 됩니다. 고을 너무 약하게(느리게) 던진다면, 공은 골대의 앞 부분을 맞고 튕겨저 나갈 것이며, 반대로 너무 강하게(빠르게) 던진다면, 공은 골대의 뒷 부분 혹은 백보드를 맞고 튕겨져 나갈 것 입니다. 공이 골대로 들어가기 위한 정확한 속도로 던져야 골인이 됩니다. 

 

위에서 설명한 일상 생활에서의 경험을 수학적으로 옮기면, 최종값 조건이 주어진 미분 방정식을 초기값이 주어진 미분 방정식의 문제로 바꿔 풀 수 있습니다. 원래 문제의 조건에서는 $t_0, t_N$에서의 함수값 $y_0, y_N$이 주어졌지만, 특정한 $y'(t_0)$을 가정하고, $t_0$에서의 함수값과 미분값 $y_0, y'_0$이 주어진 동일한 미분 방정식을 생각할 수 있습니다. 이렇게 초기값 문제인 미분 방정식을 풀고 $t_N$에서의 함수값 $y(t_N)$이 초기에 주어진 최종값 문제의 조건 $y_N$과 같다면, 가정한 $y'(t_0)$이 실제로 $t_0$의 $y(t_0)$가 되는 것이고, 이때 얻었던 해 $y(t)$가 실제 최종값 문제의 해가 되는 것 입니다. 농구 문제에서 공을 던지는 속도를 잘 맞춰야 골인이 되는 것은 $t_0$에서의 미분값 $y'(t_0)$를 잘 가정해야 $t_N$에서의 함수값 $y(t_N)$이 $y_N$이 되는 것과 같습니다. 

 

다시 한 번 깔끔하게 요약하면, $t_0, t_N$에서의 함수값 $y_0, y_N$이 주어진 경우,

 

(1) $t_0$에서의 미분값 $y'(t_0)$를 가정한다.

(2) (1)에서 가정한대로 $t_0$에서 함수값과 미분값 $(y_0, y'(t_0))$이 주어진 초기값 문제를 풀어 $y(t)$를 구한다.

(3) (2)에서 구한 해에서 $t_N$에서의 함수값 $y(t_N)$을 구한다.

(4) (3)에서 구한 $y(t_N)$이 $y_N$과 일치하는지 비교한다. 만일 일치한다면, (2)에서 구한 $y(t)$가 해가 된다. 일치하지 않는다면, (1)로 돌아가서 다른 $y'(t_0)$에 대해서 과정을 반복한다.

 

가 됩니다. 실제로 미분 방정식을 풀어야 하는 부분은 (2)로 (2)에서의 문제는 초기값 문제이기 때문에, 이전에 설명한 오일러 방법이나 룽게-쿠타 방법을 통해서 문제를 풀 수 있습니다. (2)의 과정을 문제 없이 풀 수 있다면, 남은 문제는 (4)의 과정에서 발생합니다. (1)에서 가정한 초기값이 원하는 최종값을 주지 못한다면, 초기값을 새로운 값으로 바꾸어야 하는데, 이 과정이 새로운 문제가 됩니다. 차차 설명을 하겠지만, (4)의 과정은 수치적인 방법으로 방정식(미분 방정식이 아닌 일반적인 실수의 방정식)을 푸는 것이 됩니다. 즉, 최종값의 조건이 부여된 미분 방정식은 초기값의 조건을 맞추는 방정식의 문제가 되는 것 입니다. 

 

위에서 설명한 방식으로 최종값 조건이 주어진 미분 방정식을 푸는 것을 슈팅 방법(shooting method) 이라고 합니다. 이 방법을 왜 "슈팅" 방법이라고 부르는지는 앞서 설명한 농구에서 골을 넣는 것과의 비교를 생각하면 쉽게 이해할 수 있습니다. 

 

예시 :  원하는 시간에 원하는 위치로 입자를 옮기는 궤적 찾기

 

단순한 수학적 과정만 설명을 하면 이해가 어려울 수 있으니, 간단하면서도 개념을 쉽게 설명할 수 있는 예시를 들어 보겠습니다. 아래와 같은 미분 방정식의 해를 수치적인 방법으로 구해 보도록 하겠습니다. 사실, 매우 간단한 문제라서 해석적인 방법으로 풀이가 가능하지만, 수치적인 방법으로 문제를 해결 하고, 해석적인 방법과 비교해 보도록 하겠습니다. 

$$y''(t) + y(t) = 0, t \in [t_0 = 0, t_N = 1], y(t_0) = y_0 = 0, y(t_N) = y_N = 1$$

단순히 미분 방정식이라고 생각해도 좋고, 고전 역학계와 비교해서 생각한다면, $y''(t) + y(t) = 0$ 는 $y(t)$를 시간 $t$에서 1차원 자유도를 갖는 입자의 위치라고 할 때, 단조화 진동자의 미분 방정식이 됩니다. 편의상 질량, 용수철 상수 같은 물리 상수를 모두 $1$로 둔 것 입니다. $y(t_0) = y_0 = 0, y(t_N) = y_N = 1$는 초기($t_0$) 위치는 $0$, 최종($t_N$) 위치는 $1$이 되는 경계 조건입니다. 원점에서 운동을 시작한 입자가 1초 후에 $y=1$의 위치에 있도록 하는 해를 찾는 것이 이 미분 방정식 문제 입니다. 

 

앞에서 설명한 대로, 초기 속도(미분값) $y'(t_0) = y'_0$를 적절하게 선택하여 최종 위치가 1이 되도록 만들어 주면 됩니다. 예를들어서 $y'_0=1$ 으로 놓고, 룽게-쿠타 방법을 이용하여 $t \gt 0$ 에서의 $y(t)$를 구하면,

결과는 위 그래프와 같습니다. $t=1$에서의 $y(t=1) = 0.84$입니다. 경계조건 $y(1)=1$를 만족시키기 못 합니다. 그래프의 개형을 봐도 알 수 있고, 혹은 단조화 진동자의 물리를 생각해도 알 수 있는데, 초기 속도가 작았기 때문에 $t=1$에서의 함수값이 1보다 작았습니다. 따라서, 초기 속도(미분값)을 1 보다 더 큰 수로 해야합니다. 

$y'_0=1.5$으로 초기 미분값(속도)를 조금 더 크게 했습니다. 역시 마찬가지로 룽게-쿠타 방법으로 $t=1$에서의 함수값을 구하면 $1.26$을 얻습니다. 이번에는 목표값인 $1$보다 조금 큰 값을 얻었는데, $y'_0$값이 컸기 때문입니다. 따라서 $y(1)=1$이 되는 해를 찾으려면 $1 \lt y'_0 \lt 1.5$가 되어야 합니다. 

$y'_0=1.19$로 놓고 미분 방정식을 풀어서 $t=1$에서의 함수값을 구하면 우리가 원했던 $y(1) = 1$을 얻을 수 있습니다! "슈팅이 성공" 성공 했다고 볼 수 있습니다. 해석적인 방법이 아닌 수치적인 방법으로 해를 구한 것 이기 때문에, 이 해에 대한 수식은 따로 없고, 위 그래프의 녹색 그래프가 해 그 자체 입니다.

 

앞서 언급했던 것과 같이 위 문제는 해석적인 풀이가 가능한데, 익히 잘 알려진대로 $y(t) = A \sin t$ 가 됩니다. $y(1)=1$이 되기 위해서는 $A = \frac{1}{\sin 1}$이 ㄷ되어야 하고, 따라서 $A = 1.1883$이 됩니다. $y'(t) = A \cos t$이고 $y'(0)$을 구하면 그 값은 $1.1883$이 되는데, 수치적인 방법으로 구한 해 $1.19$와 매우 유사합니다. 위에서 푼 수치적인 방법론이 제대로 작동하였음을 해석적인 해와 비교를 통해서 확인하였습니다.

 

위 과정에서 해(그래프의 녹색 선)을 구한 과정을 다시 표현하면,

 

(1) $y'_0$에 대한 $y(1)$을 수치적인 방법으로 계산. 즉 $y(1) = f(y'_0)$ 을 만족시키는 함수 $f$를 수치적인 방법으로 계산한 것과 동일함

(2) $f(y'_0) = 1$ 이 되게 하는 $y'_0$을 찾기. 즉 방정식 문제를 풀기

 

입니다. 최종적으로 풀어야 하는 것은 (2)의 과정인데, (1)의 과정은 (2)의 과정에서 풀어야할 방정식 $f(y'_0)=1$의 함수 $f$를 생성하는 과정이었습니다. 앞에서 여러번 언급 한 것 처럼, 미분 방정식을 푸는 문제는 최종값 조건을 만족하게 하는 방정식을 푸는 것으로 바뀌었고, 방정식의 해를 찾는 일반적인 방법을 통해서 원하는 $y'_0$을 찾는 과정으로 문제를 풀 수 있습니다. 방정식의 해 $y'_0 = 1.19$를 찾을 때에는 다양한 방법을 적용할 수 있는데, 이번 포스팅에서는 이를 설명하고 적용하는 것이 주제가 아니라서 적당히 $y'_0$값을 바꾸어 가면서 $f(y'_0)=1$를 만족하는 $y'_0$를 찾았습니다. 

실제로 $f(y'_0) = y(1)$의 그래프를 그리면 위 그래프의 파란색 함수와 같습니다. 물리적 직관으로 $y'_0$가 커지면 $f(y'_0)$ 역시 커진다는 것을 알 수 있고, 빨간색으로 그린 $f(y'_0) = y(1) = 1$와의 교점이 우리가 찾는 $y'_0$값이 되는데, 약 $1.19$가 됨을 얻을 수 있습니다. 

 

슈팅 방법을 통해 문제를 한 번 풀어 봤습니다. 지금은 매우 간단한 문제를 풀었는데, 개념적으로는 위 설명이 거의 전부 입니다. 복잡한 미분 방정식이 주어졌다고 하더라도, 수치적으로 문제를 푸는데는 별다른 어려움이 없기 때문에(프로그램 코드에서 수식만 조금 바꾸면 됩니다) 같은 방식으로 풀 수 있습니다. 미분 방정식의 문제에 따라서 위에서 $f(y'_0)$이라고 쓴 함수가 달라지게 됩니다. 

 

조금 더 복잡한 문제 풀어 보기 

 

이번에는 조금 더 복잡한 문제를 풀어 보도록 하겠습니다. 

$$y''(t) = -y^3 + y', y(0) = 0, y(\pi) = 0$$

문제 입이다. $y^3$항이 있기 때문에, 비선형 미분 방정식입니다. 해석적인 방법으로는 풀기가 매우 어려운 문제 입니다. (해석적인 방법으로 풀 수 있는지 조차 잘 모르겠습니다) 고전 역학 문제로 생각한다면, 위치에 세 제곱에 비례하는 복원력이 있고 ($-y^3$이 여기에 해당됩니다), 속도에 비례하는 힘이 속도의 방향으로 작용($y'$이 여기에 해당 됩니다) 하는 시스템에서 입자의 움직임 입니다. 만일 $y'$ 항이 없다고 한다면 단순히 주기적인 진동으로 생각할 수 있지만, $y'$항이 있기 때문에, $y(t)$를 예상하기가 조금 어렵습니다.

초기 조건이 $y(0) = 0, y'(0) = 0.1$ 일때의 해 $y(t)$를 구한 그래프 입니다. 정지 상태에 있다가 속도를 $+$ 방향으로 약간 가해준 입자의 움직임이라고 생각할 수 있는데, 전체적으로 진동 운동을 하지만 ($y^3$항에 의해서), 진동은 규칙적이지 않습니다. 진폭이 계속해서 바뀌고, 진동의 주기는 점점 짧아 집니다. $y'$힘은 속도의 방향과 같은 방향으로 작용하는데, 따라서 $y'$항이 계속해서 시스템에 에너지를 공급해 주고 있습니다. 따라서 진폭이 계속해서 증가하고 있는 것 입니다. 시스템의 에너지가 계속 증가하니, 입자의 속도 역시 점점 빨라져서 그래프는 점점 더 빠르게 진동하게 됩니다. 

 

우리 문제에서는 $y'(0)$를 적절하게 조절하여 $y(\pi) = 0$를 만들어야 합니다. 서로 다른 $y'(0)$에 대해서 $t \in [0, \pi]$ 영역에서 $y(t)$를 그려보면 아래와 같습니다. 

$y'(0) \in [0.0, 1.0]$ 에서 $0.1$ 단위로 증가 시켜 가면서 $y(t)$를 그린 결과 입니다. 당연하게도 $y'(0)= 0$인 경우, $y(t)= 0$ 이 되고, 우리가 원하는 최종값 조건을 만족시키긴 하지만, 이 해는 우리가 관심있어 하는 해는 아니기 때문에 논의에서 생략하도록 하겠습니다. 

 

그래프에서 보면 밝은 하늘색인 $y'(0) = 0.3$ 의 그래프가 $t=\pi$에서 $0.07$값을 가지는 것을 확인할 수 있습니다. $0.1$ 단위로 $y'(0)$값을 바꾼 그래프 중에서는 $y(\pi) = 0$에 가장 가까운 그래프 입니다. 하지만 그 값은 $0.07$ 이라서 조금 오차가 있는데, 따라서 $y'(0) = 0.31$ 정도가 정확히 최종값 경계 조건을 만족시킬 것 같습니다. 

위 그래프는 $y'(0)$의 값에 따른 $y(t = \pi)$값의 그래프 입니다. 이 파란색 그래프와 $x$축이 만날때(빨간색 점) 의 $y'(0)$값이 우리가 찾으려고 하는 $y'(0)$인데, 그 값은 앞에서 예측 한 값과 유사한 $0.308$ 이 됩니다. 그래프의 개형을 봐서는 $y'(0)$ 가 1보다 큰 어떤 값에서 다시 한 번 $x$축과 만날 것 같은데, 그 값을 $y'(0)$으로 사용한다면, 그 때의 해 $y(t)$ 역시 경계 조건을 만족 시킵니다. 

이를 확인하기 위해서 위 그래프의 $y'(0)$의 범위를 더 넓게 그렸습니다. $y'(0) = 1.25$ 부근에서 $x$축과 만나게 됩니다. 즉, $y'(0) = 1.25$ 초기 조건도 $y(\pi) = 0$의 최종값 경계 조건을 만족하게 됩니다. 물론 같은 맥락으로 $y'(0)$가 더 큰 값에서 경계 조건을 찾는 $y'(0)$를 찾을 수도 있습니다. 처음에 찾은 $y'(0) = 0.308$ 해는 입자가 $1/2$번 왕복 운동을 하게 하는 조건이 되고, 두번째  찾은 $y'(0) = 1.25$는 입자가 $1$번 완전히 왕복 운동을 하는 경계 조건 입니다. 초기 속도에 따라서 왕복 운동의 주기가 다르기 때문에 생기는 재미있는 현상입니다. 단조화 진동자에서는 초기 속도와 상관 없이 주기 운동의 주기가 같기 때문에 이와 같은 현상은 생기지 않습니다. 슈팅 방법을 이용하여 2계 미분 방정식을 풀 때는, 고전 역학의 직관을 이용할 수 있어서 매우 좋습니다. 

 

슈팅 방법을 이용하여 (주로) 물리학에서 나오는 미분 방정식 문제를 어떻게 풀 수 있고, 그 해를 어떻게 해석할 수 있는지에 대해서 다음 포스팅을 통해서 더 공부해 보도록 하겠습니다.