📂수치해석

사다리꼴 메소드

정의 ¹

$D \subset \mathbb{R}^2$ 에서 정의된 연속함수 $f$ 에 대해 초기값 문제 $\begin{cases} y ' = f(x,y) \\ y( x_{0} ) = Y_{0} \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} + {{h} \over {2}} [ f ( x_{n} , y_{n} ) + f ( x_{n+1} , y_{n+1} ) ] $$

설명

예측자-수정자 알고리즘

오일러 메소드보다 데이터를 많이 쓴 미드포인트 메소드가 더 좋은 퍼포먼스를 보이는 것처럼 사다리꼴 메소드는 계산을 더 많이 한만큼 퍼포먼스가 향상되었을거라고 기대할 수 있다. 게다가 일단은 원스텝 메소드다. 문제는 이것이 임플리시트 메소드라는 것인데, 바로 풀어낼 수가 없으므로 방정식을 푸는 메소드를 반복하면서 풀어내는 것이 보통이다. 당연하지만 계산량이 상당히 많아지며, 원스텝 메소드임에도 불구하고 속도면에서 불리한 부분이 있다.

우선은 $y_{n+1}^{(0)}$ 를 제법 $Y_{n+1}$ 에 잘 맞도록 ‘찍어서’ $$ y_{n+1}^{(1)} = y_{n-1} + {{h} \over {2}} [ f ( x_{n} , y_{n} ) + f ( x_{n+1} , y_{n+1}^{(0)} ) ] $$ 을 얻는다. 잘 맞도록 ‘찍는다’는 건 오일러 메소드 같은 것으로 한 스텝 앞의 얼추 괜찮은 추정값 $y_{n+1}^{(0)}$ 을 얻는 것을 말한다. 그리고 충분히 정확할 때까지 계산을 반복해 $$ y_{n+1}^{(j+1)} = y_{n-1} + {{h} \over {2}} [ f ( x_{n} , y_{n} ) + f ( x_{n+1} , y_{n+1}^{(j)} ) ] $$ 를 찾아내는 식이다. $y_{n+1}^{(j)}$ 에서 $j$ 는 이터레이션을 돌린 횟수를 나타내며, 반복할수록 정확해질것이다. $$ y_{n+1} = y_{n-1} + {{h} \over {2}} [ f ( x_{n} , y_{n} ) + f ( x_{n+1} , y_{n+1} ) ] $$ 의 양변에서 $$ y_{n+1}^{(j+1)} = y_{n-1} + {{h} \over {2}} [ f ( x_{n} , y_{n} ) + f ( x_{n+1} , y_{n+1}^{(j)} ) ] $$ 를 빼면 $$ y_{n+1} - y_{n+1}^{(j+1)} = {{h} \over {2}} [ f ( x_{n+1} , y_{n+1} ) - f ( x_{n+1} , y_{n+1}^{ (j) } ) ] $$ 인데, 여기서 립시츠 조건을 가정한다면 $$ | y_{n+1} - y_{n+1}^{(j+1)} | = {{hK} \over {2}} | y_{n+1} - y_{n+1}^{(j)} | $$ 을 얻는다. 즉 $\displaystyle {{ h K } \over {2}}<1$ 이 되어야 수렴을 할텐데, $K$ 가 크다면 $h$ 는 어지간히도 작게 주어져야 한다.

또한 사다리꼴 메소드가 $$ Y_{n+1} = Y_{n-1} + {{h} \over {2}} [ f ( x_{n} , Y_{n} ) + f ( x_{n+1} , Y_{n+1} ) ] - {{h^3 } \over {12}} Y^{(3)} ( \xi_{n}) $$ 에서 유도되어 $O( h^3 )$ 에 의존하는 절삭 오차를 가지므로, $|y_{n+1} - y_{n+1}^{(j)} |$ 는 못해도 $O(h^4 )$ 보다는 작아줘야 한다.설명을 읽어보면 알겠지만 위와 같은 설명대로 코드를 짠다면 본질적으로는 오일러 메소드와 다르지 않고, 실제로도 오일러 메소드를 포함하고 있다. 다만 그것을 사다리꼴 메소드를 통해서 보정한다고 보는 것이 적절하다. 그래서 이런 알고리즘을 프레딕터-커렉터 알고리즘^{predictor-Corrector Algorithm}이라 부른다. 여기서 예측자Predictor의 역할은 오일러 메소드가, 수정자Corrector의 역할은 사다리꼴 메소드가 맡고 있다. 적절한 값을 예측하고 그것을 수정하는 식으로 수치적인 해를 찾는 것이다.

구현

위 스크린샷은 초기값 문제 $\begin{cases} \displaystyle y ' = {{1} \over {1 + x^2 }} - 2y^2 \\ y(0) = 0 \end{cases}$ 를 사다리꼴 메소드와 오일러 메소드로 풀고 트루 솔루션 $\displaystyle Y = {{x } \over {1 + x^2}}$ 으로 오차를 비교한 결과다. 당연하지만 사다리꼴 메소드의 오차가 훨씬 작다.

위 스크린샷은 초기값 문제 $\begin{cases} \displaystyle y ' = x - y^2 \\ y(0) = 0 \end{cases}$ 사다리꼴 메소드와 미드포인트 메소드로 풀고 그 수치해를 비교한 결과다. 미드포인트 메소드는 뒤로 가면서 값이 요동치는 걸 보아 패러사이틱 솔루션이 있는데, 사다리꼴 메소드는 안정적으로 문제를 잘 풀어주고 있다.

아래는 R로 작성한 코드다. j는 이터레이션을 반복할 횟수로써, 따로 입력하지 않으면 오일러 메소드로 한번 찍기만 하고 넘어간다.

Euler<-function(f,Y_0,a,b,h=10^(-3))
{
  Y <- t(Y_0)
  node <- seq(a,b,by=h)
  
  for(x in node)
  {
    Y<-rbind(Y,Y[length(Y[,1]),]+h*f(x,Y[length(Y[,1]),]))
  }
  
  return(Y)
}
 
Midpoint<-function(f,Y_0,a,b,h=10^(-3))
{
  Y <- t(Y_0)
  Y<-rbind(Y,Y[1,]+h*f(a,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)
}
 
Trapezoidal<-function(f,Y_0,a,b,h=10^(-3),j=0)
{
  Y <- t(Y_0)
  Y<-rbind(Y,Y[1,]+h*f(a,Y[1,]))
  node <- seq(a,b,by=h)
  
  for(x in node[-1])
  {
    Y_guess<-Y[length(Y[,1]),] + h* f(x,Y[length(Y[,1]),] )
    if(j>0)
    {
      for(temp in 1:j)
      {
        Y_guess<-Y[length(Y[,1]),]+(h/2)*( f(x,Y[length(Y[,1]),] )+ f(x+h,Y_guess) )
      }
    } else {"Wrong Iteration Number"}
    Y<-rbind(Y,Y[length(Y[,1]),]+(h/2)*( f(x,Y[length(Y[,1]),] )+ f(x+h,Y_guess) ))
  }
  
  return(Y)
}
 
Y<-function(x) {x / (1 + x^2)}
 
f<-function(x,y) {1/(1+x^2) + - 2*(y^(2))}
out<-Trapezoidal(f,seq(0,0,len=1),0,2,h=0.1)
abs(Y(seq(0,2.1,by=0.1)) - out[,1])
out<-Euler(f,seq(0,0,len=1),0,2,h=0.1)
abs(Y(seq(0,2.1,by=0.1)) - out[,1])
 
g<-function(x,y) {x-y^2}
out<-Trapezoidal(g,seq(0,0,len=1),0,3.25,h=0.25)
out[,1]
out<-Midpoint(g,seq(0,0,len=1),0,3.25,h=0.25)
out[,1]