📂수치해석

미드포인트 메소드

메소드

$D \subset \mathbb{R}^2$ 에서 정의된 연속함수 $f$ 에 대해 초기값 문제 $\begin{cases} y ' = f(x,y) \\ ( y( x_{0} ), y (x_{1}) ) = ( Y_{0} , Y_{1} ) \end{cases}$ 가 주어져 있다. 구간 $(a,b)$ 를 $a \le x_{0} < x_{1} < \cdots < x_{n} < \cdots x_{N} \le b$ 와 같은 노드 포인트들로 쪼갰다고 하자. 특히 충분히 작은 $h > 0$ 에 대해 $x_{j} = x_{0} + j h$ 이라고 하면 초기값 $y_{0} = Y_{0}$ 에 대해 $$y_{n+1} := y_{n-1} + 2 h f ( x_{n} , y_{n} )$$

유도 ¹

$Y’(t)$ 를 $x_{n}$ 에 대해 테일러 전개하면 $$Y’(t) = y ' ( x_{n} ) + ( t - x_{n} ) Y’’ (x_{n} ) + {{( t - x_{n} )^2} \over {2}} Y^{(3)} ( \xi_{n} )$$ 를 만족하는 $\xi_{n} \in [ x_{n-1} , x_{n} ]$ 이 존재한다. 양변에 $t= x_{n-1}$ 부터 $x_{n+1}$ 까지의 정적분을 취하면 $$\int_{x_{n-1} }^{ x_{n+1} } Y’(t) dt = \int_{x_{n-1} }^{ x_{n+1} } \left[ y ' ( x_{n} ) + ( t - x_{n} ) Y’’ (x_{n} ) + {{( t - x_{n} )^2} \over {2}} Y^{(3)} ( \xi_{n} ) \right] dt$$ $x_{n-1} = x_{n} - h$ 이고 $x_{n+1} = x_{n} + h$ 이므로 $$\int_{x_{n-1} }^{ x_{n+1} } Y’(t) dt = ( x_{n} + h - (x_{n} - h) ) y ' ( x_{n} ) + {{ h^3 - ( - h^3 ) } \over {6}} Y^{(3)} ( \xi_{n} )$$ 미적분학의 기본정리에 의해 $$Y(x_{n+1}) - Y(x_{n-1})= 2h y ' ( x_{n} ) + {{ h^3 } \over {3}} Y^{(3)} ( \xi_{n} )$$ 가정에서 $Y’ = f$ 고, $\displaystyle {{ h^3 } \over {3}} Y^{(3)} ( \xi_{n} )$ 을 탈락시키면 $$y_{n+1} = y_{n-1} + 2 h f ( x_{n} , y_{n} )$$

■

오일러 메소드와의 비교

오차 분석

오일러 메소드와 비교하자면 두 개의 데이터를 사용하니까 계산은 많아지겠지만 오차가 더 작을 것이라고 짐작할 수 있다. 실제로 오일러 메소드의 오차는 조건이 많이 주어져 있어도 $O (h)$ 에 바운드되어있을 뿐이다.

반면 미드포인트 메소드는 멀티 포인트 메소드로써 그 오차를 계산하면 $$\max_{ x_{0} \le x_{n} \le b} | Y (x_{n} ) - y_{h} (x_{n} ) | \le e^{2K(b-x_{0})} \eta (h) + \left[ {{ e^{2K (b - x_{0} ) } - 1 } \over {2K}} \right] \left[ {{1} \over {3}} h^2 | Y^{(3)} |_{\infty} \right]$$ 여기서 초기오차는 $\displaystyle \eta (h) = \max \left\{ |Y_{0} - y_{0}| , | Y(x_{1}) - y_{h} (x_{1}) | \right\}$ 인데, $y_{0} =Y_{0}$ 이므로 사실 $\eta (h) = | Y(x_{1}) - y_{h} (x_{1}) |$ 이 된다. $Y$ 를 테일러 전개하면 $$Y (x_{1}) = Y(x_{0}) + h Y(x_{0}) + {{h^{2} } \over {2}} Y’’ (\zeta) = y_{0} + h f(x_{0} , y_{0}) + {{h^{2} } \over {2}} Y’’ (\zeta)$$ 이때 오일러 메소드에 따라 $y_{0} + h f(x_{0} , y_{0}) = y_{1}$ 이므로, $$Y(x_{1}) - y_{h} (x_{1}) = {{h^{2} } \over {2}} Y’’ (\zeta) = O(h^{2} )$$ 결과적으로 미드포인트 메소드는 $O (h^2)$ 에 바운드되어 오일러 메소드보다 정확도 면에서 나은 점이 있다고 말할 수 있다.

안정성 이슈

하지만 미드포인트 메소드는 치명적인 결점이 있다. 다음은 R 로 미드포인트 메소드를 구현한 코드다.

Midpoint<-function(f,Y\_0,Y\_1=NA,a,b,h=10^(-3))
{
  Y <- t(Y\_0)
  if(is.na(Y\_1))
  {
    Y<-rbind(Y,Y[1,]+h*f(a,Y[1,]))
  }
  else
  {
    Y<-rbind(Y,Y\_1)
  }
  node <- seq(a,b,by=h)
  
  for(x in node[-1])
  {
    Y<-rbind(Y,Y[length(Y[,1])-1,]+2*h*f(x,Y[length(Y[,1]),]))
  }
  
  return(Y)
}
 
f<-function(x,y) {1/(1+x^2) + - 2*(y^(2))}
out<-Midpoint(f,seq(0,0,len=1),0,2,h=0.1)
out[,1]
 
g<-function(x,y) {x-y^2}
out<-Midpoint(g,seq(0,0,len=1),0,3.25,h=0.25)
out[,1]

실행시켜보면 함수 f()와 달리 g()는 솔루션이 제대로 구해지지 않는데, 이는 미드포인트 메소드의 안정성Stability에 문제가 있음을 보여주는 예시다.