이런저런 생각을 하다가 미분방정식의 최종 시간에서의 값이 정해질 때 역으로 풀면 (즉, 시간을 거꾸로 해서 풀면) 풀릴 지에 대한 궁금증이 생겨서 테스트한 것이다. % Solve differential equation using MATLAB odesolver tspan = [0 5]; y0 = 10; [t,y] = ode45(@(t,y) -2*y, tspan, y0); figure; plot(t,y,'-o'); dt = 0.01 ; N = tspan(2) / dt ; x = zeros(N+1,1) ; f = @(t,y) -2*y ; x(1) = y(end) ; % solve reversely for i =1:N t_tmp = 5 - dt * i ; x(i+1) = x(i) - dt * f(t_tmp,x(..

Introduction 설명을 시작하기 앞서 우리가 구하고자 하는 Ordinary Differential Equation(ODE)는 다음과 같습니다. $y' = f(y,t)$ 저번 글에서 작성했던 RK4 식은 아래와 같습니다. $y^{n+1}=y^{n}+\frac{1}{6} k_{1}+\frac{1}{3}\left(k_{2}+k_{3}\right)+\frac{1}{6} k_{4}$ $k_{1}=\Delta t f\left(y^{n}, t^{n}\right)$ $k_{2}=\Delta t f\left(y^{n}+\frac{1}{2} k_{1}, t^{n}+\frac{\Delta t}{2}\right)$ $k_{3}=\Delta t f\left(y^{n}+\frac{1}{2} k_{2}, t^{n}+\fra..

ODE에 대한 Introduction에서 적분을 할 때 정해놓은 $\Delta t$가 아니라 그 안에서 여러 스텝을 거쳐서 적분을 하는 방법은 substep method이고, $t^{n+1}$라는 미래 시간을 계산하기 위해 현재와 그 이전 스텝을 쓰는 것을 multi step method라고 했습니다. 이번 글과 다음 글에서 소개할 방법은 multi step method의 범주에 드는 방법들입니다. 위의 그림처럼 $\Delta t$ 간격의 시간에서 $t^{n-1}$ 등 과거 스텝도 사용합니다. 이 방법에는 몇 가지 단점이 있습니다.첫 번째, timestep size가 고정됩니다.두 번째, self starting이 불가능합니다. 즉, 초기 조건을 알아도 그 초기 조건을 현재로 생각했을 때 이전 스텝도 알..

1. Introduction Runge Kutta 방법(이하 RK 방법)은 explicit method 중에서 nonlinear에 적합하기 때문에 많이 쓰입니다. substep method라고 해서, RK방법은 한 타임스텝을 한 번에 계산하는 것이 아니라 여러 스텝으로 나눠서 진행하는 방법입니다. Advantage 1) Stability가 좋다 (일반 explicit method에 비해 stability 영역이 넓음) 2) time step size를 바꿔가면서 사용할 수 있다 (이 글에서는 균일한 간격을 사용하지만 달라도 가능하다) 3) self-starting (이후에 나오는 multistep method는 초기 조건만으로 시작할 수가 없지만 이 방법은 가능) RK방법도 2차, 3차, 4차 등 차수를..

글에서 수치해석 방법을 줄여서 씁니다. EE : Explicit Euler IE : Implicit Euler TR : Trapezoidal method Explicit method가 implicit method 대신 쓰이는 이유 지금까지 Explicit Euler, Implicit Euler, Trapezoidal method를 살펴보았는데 trapezoidal method는 accuracy에서나, stability에서나 좋은 점이 많습니다. EE, IE 방법보다 정확도도 1차 높고 stability 영역도 exact solution과 동일합니다. 하지만 그럼에도 불구하고 Trapezoidal method를 쓰지 않고 Predictor-Corrector method(PC), 또는 Runge-Kutta ..

(앞선 글부터 보면 알 수 있지만 문자를 표현할 때 아래 첨자는 공간 상의 좌표, 위 첨자는 시간 상의 좌표를 의미합니다.) trapezoidal 방법은 implict method 중 하나입니다. (즉, 미분할 때 다음 값($y^{n+1},t^{n+1}$)도 가져다 씀) 이 방법은 PDE에 쓰일 때 Crank-Nicolson method라고 합니다. $y^{'}=f(y,t)$ $\frac{y^{n+1}-y^{n}}{\Delta t}=\frac{1}{2}\left[f(y^{n},t^{n})+f(y^{n+1},t^{n+1})\right]$ implicit이면 stable 영역에서 항상 stable합니다. 위 방법의 장점은 accuracy가 globally 2nd order($\Delta t^{2}$)라는 점..

이번에는 explicit euler 방법이 아닌 implicit euler 방법을 사용하고자 합니다. # Notation 참고 $x_{B}^{A}$ : $A$는 시간에서의 위치(Timesteps), $B$는 공간에서의 위치를 의미합니다. $y^{'}=f(y,t)$ model problem인 $y^{'}=\lambda y$을 implicit euler 방법으로 풀어보겠습니다. $\frac{y^{n+1}-y^{n}}{\Delta t}=f(y^{n+1},t^{n+1})$ : implicit euler이므로 현재 y,t값이 아닌 미래의 값(n+1)을 사용합니다. $\frac{y^{n+1}-y^{n}}{\Delta t}=\lambda y^{n+1} \Rightarrow y^{n+1}-\lambda \Delta t..

Explicit euler 방법에 대해서 설명하겠습니다. (이전 글에서도 다룬 적이 있어서 링크를 달아놓겠습니다.) $y^{'}=f(y,t)$이라고 할 때 시간 축에서 한 step 움직였을 때, y값을 구하고자 합니다. 이전 step의 값을 알고 있다고 가정하고, taylor expansion을 통해서 구합니다. $y^{n+1}=y(t^{n}+\Delta t)=y^{n}+\Delta t (y^{n})^{'}+\frac{\Delta t^{2}}{2!}(y^{n})^{''}+\frac{\Delta t^{3}}{3!}(y^{n})^{'''}+\cdots$ 이 식을 $y^{'}$에 대해서 정리합니다. $(y^{n})^{'}=\frac{y^{n+1}-y^{n}}{\Delta t}-\frac{\Delta t}{2!..

이번에는 ODE를 수치해석적인 방법으로 푸는 것에 대해 다루겠습니다. 1. Introduction 예를 들어 $y^{''}+wy=f(x)$ $\frac{d^{2}y}{dx^{2}}+wy=f(x)$ 적분 상수가 2개 나오므로 조건이 2개가 필요합니다. initial value problem(IVP)이면 초기조건이 주어지고 $y(0)=y_{0}$ $\left. \frac{dy}{dx} \right |_{x=0}=v_{0}$ boundary value problem(BVP)이면 $y(0)=y_{0}$ $y(L)=y_{L}$ 일단 처음에는 IVP를 먼저 풀고, 그 다음에 BVP로 확장하겠습니다. 모든 방법은 $0\leq t \leq t^{n}$에서의 solution을 알고 있다고 가정하고 그 solution을 ..

이번에는 numerical integration의 마지막, Gauss quadrature에 대해서 알아보겠습니다. 1. Definition of Gauss quadrature Gauss quadrature는 특정 weight($w_{i}$)와 함수값 ($f(x_{i})$)의 곱을 discrete하게 더했을 때 그 구간의 definite integral과 거의 approximate하게 만드는 방법입니다. 식으로 표현하면 다음과 같습니다. $\int_{a}^{b} f(x) dx \approx \sum_{i=0}^{n}f(x_{i})\cdot w_{i}\cdots\cdots (*)$ 이 때, $h\text{ (interval between the points)}, \{w_{i}\}, \{x_{i}\}$는 균일..