[수치해석] System of ODEs - Initial value problem

 

이 글에서는 여러 개의 state가 있는 ODE를 푸는 과정에 대해 다룬다.

 

 이 글만 봐서는 이해가 안되는 맥락의 내용이 많은데 이런 것은 앞에서 다룬 여러 방법들을 살펴봐야한다.

 

https://normal-engineer.tistory.com/203

[수치해석] Numerical solution of ODE (7) Leap frog method

 

 

1. Formulation

가장 일반적인 system of ODEs는 다음과 같이 state $y$가 벡터인 형태로 나타난다.

 

$\mathbf{y}^{'}=\mathbf{f}( \mathbf{y},t) $

 

ODEs 중에서도 선형인 경우에는 쉽게 분석을 할 수 있다.

$\frac{dy}{dt}=Ay=S^{-1}\Lambda S y$

$A$ matrix는 eigen decomposition에 따라 $S^{-1}\Lambda S$로 바꿀 수 있다.

 

위 식을 다시 정리하면,

$\frac{d Sy}{dt}=\Lambda (Sy)$

$\frac{dz}{dt}=\Lambda z$

여기서 $z$는 벡터이다. 따라서

각 eigenvector에 대하여

$z^{'}_{i}=\lambda_i z_i$로 표현할 수 있다.

각 $z_i$에 대하여 앞서 다뤘던 1차원 미분방정식과 동일하게 분석할 수 있다.

 

$z^{n+1}-z^{n}=\Lambda \Delta t z^{n}$

$z_{i}^{n+1}=(1+\lambda_i \Delta t)z_{i}^n$이므로

$z_{i}^{n}=(1+\lambda_i \Delta t)^n z_i^0=\sigma_i^n z_i^0$

 

따라서 $\sigma_i$에 따라 System of ODEs의 stability가 결정된다.

 

$\lambda_i = \lambda_{iR}+i\lambda_{iI}$로 complex value라고 할 때

$\| (1+\lambda_R \Delta t) + i (\lambda_I \Delta t) \| \leq 1$를 만족시켜야 numerical solution이 stable하게 나온다.

 

 

위 그림에서 x축은 $\lambda_R \Delta t $, y축은 $i\lambda_{I}\Delta t$를 의미한다.

 

따라서 $\Delta t$ (time step)은 $2/|\lambda_{i,max}|$를 만족시켜야 한다.

여기에서도 condition number 개념을 적용할 수 있다.

$K(A)=\frac{ | \lambda_{max} |}{ |\lambda_{min}| } $

matrix가 커지만 $K$도 커지는 경향이 있다. 즉, ODE 갯수가 많아지면 condition number가 커진다. stiffness가 커지면 $\Delta t$가 작아져야 한다.

 

예를 들어

$u^{'}=998u+1998v$

$v^{'}=-999u -1999v$

$u(0)=1, v(0)=1$

$\lambda_1 = -1, \; \lambda_2=-1000$

 

따라서 $u=4e^{-t}-3e^{-1000t}, v=-2e^{-t}+3e^{-1000t}$로 나온다. 즉, 응답 속도가 매우 다른 두 개의 term이 state들에 들어오기 때문에 그만큼 $\Delta t$가 매우 작아져야한다.

 

2. Nonlinear problem

 

Nonlinear Problem일 경우에는 다른 분석을 적용해야한다. 위와 같은 선형분석이 안되기 때문이다.

 

2-1) Linearization

 

$\mathbf{y}^{'}=\mathbf{f}( \mathbf{y},t) $ (nonlinear)

이 식을 implicit Euler를 적용한다고 해보자.

 

$y^{n+1}=y^n+f(y^{n+1},t^{n+1})$

$f$가 비선형이기 때문에 

 

$f(y^{n+1})=f(y^n)+(y^{n+1}-y^{n}) \left. \frac{\partial f}{\partial y} \right|_{y=y^{n}} + \frac{ (y^{n+1}-y^{n})^2 }{ 2! } \left. \frac{\partial^2 f}{ \partial y^2} \right|_{y=y^n}$

 

Trapezoidal 방법으로 $y^{n+1}$를 표현하면 다음과 같다.

 

$y^{n+1}=y^n +\frac{\Delta t}{2} \left[ f(y^{n+1},t^{n+1})+ f(y^n,t^n) \right]$

여기에 있는 $f(y^{n+1},t^{n+1})$를 taylor expansion한 값으로 대체한다.

 

$$ y^{n+1}=y^n + \frac{\Delta t}{2}  \left[ f(y^n,t^{n+1})+(y^{n+1}-y^{n}) \left. \frac{\partial f}{\partial y} \right|_{(y^n, t^{n+1})}+f(y^n, t^n)  \right] +\theta(\Delta t^3)$$

 

이제 우변, 좌변을 $y^{n+1}$과 $y^{n}$에 대해서 정리할 수 있게 된다. 대신 시간 $t^n, t^{n+1}$을 잘 구별해야 한다.

 

$$\left( 1-\frac{\Delta t}{2} \left. \frac{\partial f}{\partial y} \right|_{y^n, t^{n+1}}\right) y^{n+1} =y^n+\frac{\Delta t}{2} \left[ f(y^n, t^{n+1})+f(y^n, t^n) \right] -y^n \frac{\Delta t}{2} \left. \frac{\partial f}{\partial y} \right|_{(y^n, t^{n+1})} $$

 

 

이를 n개 state로 확장하면 다음과 같다.

 

$$ \frac{d u}{d t}=f\left(u_1, u_2, \ldots, u_m\right)=f(u) $$

state가 n개라고 하자.

 

앞선 방식으로 trapezoid 방법으로 이산화한다.

 

$$ u^{n+1}=u^n+\frac{\Delta t}{2}\left[f\left(u^{n+1}\right)+f\left(u^n\right)\right]  $$

$$ f_i\left(u^{n+1}\right)=f_i\left(u^n\right)+\left.\left(u_1^{n+1}-u_1^n\right) \frac{\partial f_i}{\partial u_1}\right|^n+\left.\left(u_2^{n+1}-u_2^n\right) \frac{\partial f_i}{\partial u_2}\right|^n+\cdots +\theta(\Delta t^2)$$

 

우변을 정리하면 $f_{i}(u^n)+\sum_{j=1}^m\left(u_j^{n+1}-u_j^n\right) \left. \frac{\partial f}{\partial u_j}\right|^n+\theta\left(\Delta t^2\right)$

 

이를 systematic하게 표현하면 다음과 같은 matrix 형태로 표현할 수 있다.
$$ f\left(u^{n+1}\right)=f\left(u^n\right)+A\left(u^{n+1}-u^n\right)+\theta\left(\Delta t^2\right)$$

where $A=\left[\begin{array}{lll}\frac{\partial f_1}{\partial u_1} & \frac{\partial f_1}{\partial u_2} & \cdots & \frac{\partial f_1}{\partial u_m}  \\ &  & & \\ \frac{\partial f_n}{\partial u_1} & & & \frac{\partial f_n}{\partial u_m}\end{array}\right]^n=$ Jacobian matrix

 

 

$$ u^{n+1}=u^n +\frac{\Delta t}{2} \left[ f(u^n) +A(u^{n+1}-u^n) +f(u^n) \right] +\theta (\Delta t^3)$$

$$\left( I -\frac{\Delta t}{2} A\right) u^{n+1} = \left( I+\frac{\Delta t}{2} A\right ) u^n + \Delta t f(u^n)$$

 

$$\left( I-\frac{\Delta t}{2} A\right)$$는 full matrix이다.

$\Delta t$가 크면 inversion 해야하는데 full matrix의 inverse matrix를 구하는 데에는 $\mathcal{O} (n^3)$ 시간이 걸리기 때문에 차라리 newton iteration을 사용하는 것이 나을 수 있다.

$\Delta t$가 작으면 inversion이 필요없을 수 있다.

 

2-2) Newton iteration

$f(x)$를 residual이라고 할 때 우리는 $f(x)=0$로 만들고 싶어한다. newton iteration은 수치적인 방법으로 $f$의 root를 구하기 위한 방법이다.

 

$$f^{\prime}\left(x^k\right)=\frac{f\left(x^{k+1}\right)-f\left(x^k\right)}{x^{k+1}-x^k}=x^{k}-\frac{f(x^k)}{f^{'}(x^k)}$$

 

위 식에서 $f(x^{k+1})=0$로 만들고자 한다고 하자.

$$f^{\prime}\left(x^k\right)=x^{k}-\frac{f(x^k)}{f^{'}(x^k)}$$로 계산해서 $x^{k}$와 $x^{k+1}$의 차이가 거의 없으면 반복을 중단한다.

 

다시 원래 미분방정식을 푸는 문제로 돌아오자.

$y^{'}=g(y,t)$에서 iterative하게 풀고자 할 때

$y^{n+1}=y^n+\Delta t g\left(y^{n+1}, t\right)$

$f(y)=y^{n+1}-y^n-\Delta t g\left(y^{n+1}, t\right)$가 0이 되도록 newton method를 적용해서 $y^{n+1}$을 구할 수 있다.

 

$$\underline{u}^{n+1}=\underline{u}^n+\frac{\Delta t}{2}\left[f\left(\underline{u}^{n+1}\right)+f\left(\underline{u}^n\right)\right] $$

 

$$F\left(u^{n+1}\right)= \underline{u}^{n+1}-\underline{u}^n-\frac{\Delta t}{2}\left[f\left(\underline{u}^{n+1}\right)+f\left(\underline{u}^n\right)\right]=0$$

 

$$\left.\frac{\partial F_i}{\partial u_j}\right|^k=\frac{F\left(\underline{u}^{k+1}\right)-\underline{F}\left(\underline{u}^k\right)}{ \underline{u}^{k+1}-\underline{u}^k}$$

$$A\underline{u}^{k+1}=A\underline{u}^k -\underline{F}(\underline{u}^{k})$$