📂数値解析

台形法

定義 ¹

$D \subset \mathbb{R}^2$で定義された連続関数に対する初期値問題$\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} ) ] $$

説明

予測子-修正子アルゴリズム

Eulerメソッドよりもデータを多く使ったMidpointメソッドがより良いパフォーマンスを示すように、台形メソッドは計算を行う分だけパフォーマンスが向上することが期待できる。ただし、基本的にはワンステップメソッドである。問題は、これがimplicit methodであることだが、直接解けないので、方程式を解くメソッドを繰り返し使いながら解くことが一般的である。当然ながら計算量はかなり多くなり、ワンステップメソッドであっても速度面で不利な点がある。

まずは、$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)} ) ] $$ を得る。‘当てる’とは、Eulerメソッドのようなものを使って一ステップ先のまあまあ良い推定値$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 )$よりも小さくしなければならない。説明を読めばわかるが、上述のようにコードを書けば、本質的にはEulerメソッドと差はなく、実際にそれを含んでいる。ただし、それを台形メソッドを通して修正するという考え方が適切である。そのため、このアルゴリズムは予測子-修正子アルゴリズムと呼ばれる。ここで、予測子Predictorの役割はEulerメソッドが、修正子Correctorの役割は台形メソッドが担っている。適切な値を予測し、それを修正して数値解を見つけるのが目標である。

実装

上のスクリーンショットは、初期値問題$\begin{cases} \displaystyle y ' = {{1} \over {1 + x^2 }} - 2y^2 \\ y(0) = 0 \end{cases}$を台形メソッドとEulerメソッドで解き、その誤差をtrue solution$\displaystyle Y = {{x } \over {1 + x^2}}$と比較した結果である。台形メソッドの誤差の方がはるかに小さいのは当然である。

上のスクリーンショットは、初期値問題$\begin{cases} \displaystyle y ' = x - y^2 \\ y(0) = 0 \end{cases}$を台形メソッドとMidpointメソッドで解き、その数値解を比較した結果である。Midpointメソッドは進むにつれて値が揺れるのが見られ、パラサイティックソリューションがあることを示しているが、台形メソッドは安定して問題をよく解いている。

以下はRで書かれたコードである。jは、別途入力しない場合はEulerメソッドで一度だけ予測して進む、繰り返すイテレーションの回数を示す。

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]