[수치해석] Error characteristics/Modified wave number

1. Order of accuracy / Leading order of error

 

지금까지 numerical method를 이용해 interpolation하거나 미분을 할 때 어떻게 구하는지, 또 그렇게 했을 때 error가 어떻게 되는지 알아봤습니다.

 

Order of accuracy, Leading order of error는 수치해석에서 가장 신경쓰는 요소 중 하나입니다. error의 order는 클수록 '대체로' 좋습니다.

 

$\theta (h^{p}) \rightarrow p \log \theta(h)$

p는 곧 $\log E$와 $\log (h)$의 기울기이기 때문에 p가 클수록 에러가 작아지면 더 급격하게 에러가 작아집니다.

 

만약 numerical method를 이용해 derivative를 구하려고 할 때 order of accuracy의 양상이 내가 구한 것과 다르다면 error가 inconsistent하다는 것이고, 이는 곧 numerical method를 적용할 만큼 데이터가 충분히 smooth하지 않다는 것입니다.

예를 들어 4th order of accuracy인 수치해석을 적용하는데 $\Delta$가 1/2배로 줄어들 때마다 error가 1/16배로 줄어들지 않는다면 inconsisten하다고 볼 수 있습니다.

 

따라서 수치해석 기법을 적용할 때, error가 consistent하게 적용되는지, 또 error가 발산하지 않고 수렴하는지 확인해야 합니다. 이 과정을 verification이라고 합니다. (다른 grid size에서 풀었을 때 error가 어떻게 되는지 확인한다.)

 

실험데이터와 수치해석 결과를 비교하는 것은 validation이라고 합니다. 이 두 개를 V&V라고 합니다.

 

 

---

 

지금까지 error에 대해서 다루었고, error의 order는 error가 줄어드는 속도일 뿐 그 자체로 order가 작다고 해서 항상 좋다고 말할 수는 없다고 했습니다.

 

그래서 이번에는 error의 characteristics과 modified wave number에 대해서 알아보려고 합니다.

 

 

2. Error Characteristics

 

error는 두 가지 특징으로 나눌 수 있습니다.

1) Dissipation

2) Dispersion

 

dissipation error는 실제 값보다 더 작게 만드는 error를 의미하고

dispersion error는 실제 값보다 더 진동하게 만드는 error를 의미합니다.

 

이 그림에서 보면 빨간색이 실제 값이고, estimate한 값이 파란색과 초록색일 때, 파란색은 빨간색에 비해 값이 작아진 것이고 초록색은 빨간색에 비해 진동하는 모습을 확인할 수 있습니다.

정확한 그림은 아니지만 대략적으로 dissipation error과 dispersion error를 시각적으로 확인하기 위해 첨부합니다.

 

 

error의 characteristics에 대해 알고 있으니, 더 자세한 error에 대한 이해를 위해 modified wave number에 대해 알아보도록 하겠습니다.

 

---

3. Modified Wave Number

 

다음과 같은 충분히 smooth한 함수 $f(x)$가 있다고 하겠습니다.

이 함수는 사실 여러 개의 wave로 분해할 수 있습니다.

 

 

wave number를 k라고 하고, k는 한 wave의 주기 length의 역수에 비례합니다. length가 길수록 k가 작습니다. 즉, 어떤 함수를 여러 개의 wave로 분해해서 볼 때 k가 작은 순서대로 보면, 파장이 가장 긴 순서대로 볼 수 있습니다.

위의 그림에서도 맨 위에 있는 wave number $k_{1}$이 가장 파장이 긴 경우입니다.

 

만약 위의 wave를 공간이 아닌 시간 domain에서 생각한다면 k는 frequency와 같은 개념일 것입니다. 왜냐하면 한 주기에 해당하는 시간의 역수에 비례하므로, $\frac{1}{T}=f$입니다.

 

여기에서는 공간에 대해서 생각하겠습니다.

 

$a_{1}\sin k_{1}x, a_{2}\sin k_{2}x, a_{3}\sin k_{3},\cdots$

 

 

우리가 $f^{'}(x)$를 Numerical method로 구할 때 그 결과를 보면

 

$\frac{\partial f}{\partial x}+\text{Leading Order of Error}$

 

뒤에 error가 존재합니다. 그런데 Numerical method는 discrete하게 $\Delta x$마다 포인트를 나눠서 계산하는 것이기 때문에, k가 큰 성분(즉, 파장이 매우 짧은 성분)은 같은 $\Delta x$에 대한 resolution이 충분하지 않습니다.

파장이 짧은 성분에 대해서는 계산을 해내기 어렵다는 뜻입니다.

 

따라서 numercial error는 k에 dependent하고, 우리가 error의 특성을 알기 위해 wave number를 알아야 하는 이유가 바로 여기에 있습니다.

 

* error의 order는 global한 특성이며, order가 높다고 반드시 좋은 것이 아니기 때문에 wave number를 파악할 필요가 있습니다.

 

---

수치해석을 할 때 exact wave number(analytic wave number)를 구할 수 없어서 그와 유사한 wave number를 구할 수 있습니다. 이것이 바로 modified wave number입니다.

 

modified wave number를 알면 accuracy와 error characterstics을 알 수 있습니다.

 

직접 modified wave number을 구해서 알아보겠습니다.

 

a harmonic function with period L가 있다고 하겠습니다.

$f(x)=e^{ikx}=\cos kx +i\sin kx \rightarrow \text{  (oscillatory function)}$

여기서 k는 wave number입니다. (spatial domain에서는 wave number, time domain에서는 frequency)

 

N개의 grid로 나누었을 때 maximum wave number는 정해져있습니다. (구별할 수 있는 최대의 wave number)

 

wave number $k=\frac{2 \pi}{L}n$ $(n=0,1,2,\ldots,\frac{N}{2})$

L은 우리가 고려하는 domain size를 말합니다.

maximum wave number는 $k=\frac{2\pi}{L}\frac{N}{2}=\frac{\pi N}{L}$

 

maximum wave number은 일종의 nyquist frequency와 개념적으로 비슷해보입니다.

maximum wave number : N개의 grid일 때 resolve할 수 있는 wave의 최대 wave number

 

그리고 f의 미분값은 continuous할 때 $f^{'}(x)=ike^{ikx}$로 exact differentiation을 알고 있습니다.

 

Finite difference for $f^{'}$를 구해보겠습니다.

 

2nd order central diff

 

$f_{j}^{'}=\frac{f_{j+1}-f_{j-1}}{2h}+\theta(h^{2})$

 

$=\frac{e^{ik(x_{j}+h)}-e^{ik(x_{j}-h)}}{2h}+\theta(h^{2})$

 

$=e^{ikx_{j}}\left [ \frac{e^{ikh}-e^{-ikh}}{2h} \right ]+\theta(h^{2})$

 

$=e^{ikx_{j}}\left [ \frac{i\sin kh}{h} \right ]+\theta(h^{2})$

 

$=i\frac{\sin kh}{h}e^{ikx_{j}}+\theta(h^{2})$

 

exact differential과 비교해보면 $f^{'}(x_{j})=ike^{ikx_{j}}$이므로

$\frac{i\sin kh}{h}\leftrightarrow ik$

두 값이 대응됩니다.

 

$\lim_{h\rightarrow 0}\frac{i\sin kh}{kh}\cdot k = k$입니다.

 

 

 

위 그림을 보면 exact wave number $k$와 modified wave number $k_{2CD}^{'}=\frac{\sin kh}{h}$ 둘의 차이를 보여줍니다. ($k_{2CD}^{'}$는 2nd order central일 때 modified wave number)

 

wave number가 작을 때는 $k$와 $k_{2CD}^{'}$의 차이가 작지만, wave number가 크면 오차가 크다.

 

 

이 오차는 dispersion error로 반영된다.

 

 

지금 본 것은 2nd order central diff인데 1st order forward&backward diff의 modified wave number은?

 

$f_{j}^{'}=\frac{f_{j}-f_{j-1}}{h}+\theta (h^{1}) $

 

$\cong \frac{e^{ikx_{j}}-e^{ik(x_{j}-h)}}{h}$

 

$=e^{ikx_{j}}\frac{1-e^{-ikh}}{h}$

 

$=e^{ikx_{j}} \left ( \frac{1-\cos kh+ i \sin kh}{h} \right )$

 

exact wave number는 $f^{'}(x_{j})=ike^{ikx_{j}}$를 통해서 $k$임을 알고 있습니다.

 

$\frac{1-\cos kh + i\sin kh}{h}\leftrightarrow ik$

이 두 값이 서로 대응됩니다.

 

$k_{1BD}^{'}=\frac{-i(1-\cos kh)}{h}+\frac{\sin kh}{h}=\frac{\sin kh}{h}+\frac{i(\cos kh-1)}{h}$

 

참고) $k_{2CD}^{'}=\frac{\sin kh}{h}$

 

1st order backward diff의 경우에 modified wave number는 dispersion error와 dissipation error 모두 존재합니다. 또한 $\frac{\sin kh}{h}$와 $k$의 차이를 비교하는 것 뿐 아니라 허수 부분인 $\frac{i(\cos kh-1)}{h}$이 추가되었는데, 이는 순수하게 다 error로 남습니다.

 

1st order backward diff는 dissipation error의 영향이 더 큽니다 아래는 exact wave number와 $\frac{(\cos kh-1)}{h}$의 차이를 나타낸 것입니다.

 

 

실제 문제에 적용할 경우

 

first-order PDE

$\frac{\partial u(t,x) }{\partial t}+\frac{\partial u(t,x)}{\partial x}=0$

$u(t,x)=f(t)e^{ikx}$라고 가정하면

 

$e^{ikx}\frac{df}{dx}+ikf(t) e^{ikx}=0$

위 식에서 $ik^{'}f(t) e^{ikx}$에서 numerical 방법을 썼을 때는 modified wave number이므로 $ik^{'}f(t) e^{ikx}$로 표시해야합니다.

 

위 식을 정리하면 $\frac{df}{dt}=-ikf$

$f=f(0)e^{-ikt}$

$u(t,x)=f(0)e^{ik(x-t)}$

 

위의 식이 exact solution이지만 modified wave number하고 생각하면

 

$\frac{df}{dt}=-ik_{1BD}^{'}f$이므로 (1st order backward diff일 때)

$u(t,x)=f(0)e^{ikx}e^{-ik_{1BD}^{'}t}$ 이 결과와 exact solution인 위 식의 결과를 비교해보면 오차를 확인할 수 있습니다.

 

$k_{2CD}^{'}=\frac{\sin kh}{h}$

$k_{1BD}^{'}=\frac{\sin kh}{h}+\frac{i(\cos kh-1)}{h}$

sin과 cos에 대한 taylor expansion을 통해 각 항의 order of accuracy를 구할 수 있습니다.

 

$k_{1BD}^{'}=\frac{1}{h}\left \{ (kh)-\frac{(kh)^{3}}{3!}+\theta[(kh)^{5}] \right \}+\frac{i}{h}\left \{ 1-\frac{(kh)^2}{2!}+\theta[(kh)^{4}]-1 \right \}$

$\frac{\sin kh}{h}=k+\theta(h^{2})$

$\frac{(\cos kh-1)}{h}=0+\theta(h^{1})$

 

 

 

이렇게 modified wave number를 구하는 이유는,

1) order of accuracy 측정

2) error characteristics 파악

에 목적이 있습니다.

위에서 구한 modified wave number를 solution에 대입하겠습니다.

 

$\text{2nd order central diff : }u(t,x)=f(0)e^{ikx}\exp(-i\left [ k -\frac{k^{3}h^{2}}{6}+\cdots \right ])=f(0)e^{ik(x-t)}e^{i \left ( \frac{k^{3}h^{2}}{6}\cdots \right ) t}$

 

$f(0)e^{ik(x-t)}$가 exact solution이므로 뒤에 붙은 $e^{i \left ( \frac{k^{3}h^{2}}{6}\cdots \right ) t}$이 오차를 발생시킵니다. 이 error는 오일러 공식에 따라 $\cos \left(\frac{k^{3}h^{2}}{6}t \right)+i\sin \left(\frac{k^{3}h^{2}}{6}t \right)$와 같이 cos와 sin으로 바꿀 수 있으므로 oscillratory term이고 곧, dispersive error라고 볼 수 있습니다.

 

$\text{1st order backward diff : }u(t,x)=f(0)e^{ikx}\exp(-i\left [ k -\frac{k^{3}h^{2}}{6}+\cdots \right ])\exp(\left [ \frac{-k^{2}h}{2}t+\cdots \right ])$

$=f(0)e^{ik(x-t)}e^{i \left ( \frac{k^{3}h^{2}}{6}\cdots \right ) t}e^{\left ( -\frac{k^{2}h}{2}t\cdots \right ) t}$

1BD의 경우에는 wave number에서 실수 부분과 허수 부분 둘 다 있기 때문에 따로 나눠서 계산합니다.

 

그 결과 2CD처럼 oscillatory error도 있지만

$e^{\left ( -\frac{k^{2}h}{2}t\cdots \right ) t}$와 같이 dissipation error도 존재합니다.

 

 

 

 

그리고 주의해야할 점은

error가 발산하면 안된다는 것입니다.

 

1BD의 경우 추가적으로 $e^{i \left ( \frac{k^{3}h^{2}}{6}\cdots \right ) t}$와 $e^{\left ( -\frac{k^{2}h}{2}t\cdots \right ) t}$항이 error로 작용하는데 첫 번째 dispersion term은 크기 자체에 영향을 미치지 않습니다. 하지만 두 번째 dissipation error term은 크기에 영향을 미칩니다.

 

따라서

 

$\frac{\partial u }{\partial t} + c \frac{\partial u}{\partial x}=0 \Rightarrow x-ct=constant$

다음과 같은 Partial differential equation에서 $c>0$이면 시간이 지나면, $ct$가 증가하므로 $x-ct$를 constant로 유지하기 위해 $x$도 증가해야합니다. 즉, 시간이 지남에 따라 x축 기준 증가하는 방향으로 진행됩니다.

 

반대로 $c<0$일 때는 x축 기준 감소하는 방향으로 진행됩니다.

 

 

각각의 경우에 따라 적절한 numerical scheme이 달라집니다.

 

 

1) $c>0$일 때

 

1BD를 적용하면 위에서 했던 대로 $e^{\left ( -\frac{k^{2}h}{2}t\cdots \right ) t}$항으로 인해 decaying합니다.

그러나 1FD(1st order forward finite difference)를 적용하면, 이전에 다루었던 포스팅(normal-engineer.tistory.com/97)에서 알 수 있는 것처럼

1BD와 1FD는 error term의 부호가 반대이기 때문에  $e^{\left ( -\frac{k^{2}h}{2}t\cdots \right ) t}$항 대신  $e^{\left ( +\frac{k^{2}h}{2}t\cdots \right ) t}$으로 시간이 지남에 따라 error term이 발산하게 됩니다.

 

 

2) $c<0$일 때

 

x의 진행방향이 $c>0$일 때와 반대이기 때문에 1FD를 적용하는 것이 더 적절합니다.

 

이렇게 진행방향과 거꾸로 가는 방향의 scheme을 사용하는 것을 upwind method라고 합니다.

각 방법의 정확도와 별개로, solution이 발산하는 것은 적절하지 않기 때문에 많은 상용프로그램에서도 default로 사용하는 방법이라고 합니다.

 

 

---

modified wave number까지 얘기했으니 이제 numerical method에서 interpolation, differentiation를 배웠고 마지막으로 integration을 정리하겠습니다.