📂数値解析

台形法

定義¹

$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]