数列に対する加速法

記事の内容


今回は少し地味ですが, 数列に対する加速法をとりあげます. 代表的な加速法であるRichardsonの加速法とAitkenの加速法を実験します.

加速法の概要

概要

加速法は, 数列や級数の極限への収束を早くする方法です. Newton法などの反復法や, 数値積分のような級数計算において利用されます. また, 和算においては古くから知られていことでも有名です. 今回はこの加速法の中でも代表的なものを2つご紹介します. この記事での議論は厳密ではないので, 詳しくは[1]を参照してください.

Richardosnの加速法

ここでは, 1つ目の加速法, Richardosnの加速法について基本的なアイデアを紹介し, 実験します.

アイデア

Richardsonの加速法の基本的なアイデアを述べます. 収束する数列\(\{ x_n\}_{n=0}^{\infty}\)を考えます. この数列の極限値を\(x^*\)とします. この数列が以下のように, 近似的に, かつ漸近的に"極限値"と"減衰する成分"に分解できるとします.

\begin{equation} x_n \simeq x^* + C\lambda_1^n,\quad |\lambda_1|< 1 \end{equation}

数列の収束が遅い場合, この減衰項の影響を考慮した方が良さそうです. まず, この減衰項のオーダーを見積もるため, \(\lambda_1\)を求めます. 十分大きいnに対して, \begin{equation} \lambda_1 = \frac{\lambda_1^{n+1}}{\lambda_1^n} \simeq \frac{x_{n+1}-x^*}{x_{n}-x^*} \end{equation} として求められます. そこで, 次のように定めます.

\begin{equation} \lambda = \underset{n\to\infty}{\lim} \frac{x_{n+1}-x^*}{x_{n}-x^*} \end{equation}

次に, この\(\lambda\)を用いて, 元の数列\(\{ x_n\}_{n=0}^{\infty}\)よりも早く収束する数列を構成します. 次のように定めます.

\begin{equation} x_{n+1}^{(1)} = \frac{x_{n+1}-\lambda x_n}{1-\lambda} \end{equation}

同様に, \(\{ x_n\}_{n=0}^{\infty}\)を\(\{ x_n^{(1)}\}_{n=0}^{\infty}\)に置き換えることでさらに早く収束する数列\(\{ x_n^{(2)}\}_{n=0}^{\infty}\)が構成でき, その後も引き続き構成できます.

さて, 上式のような数列はなぜ収束が早くなる(と期待できる)のでしょうか?次のように分解して考えます.

\begin{equation} x_{n+1}^{(1)} = \frac{x_{n+1}-\lambda x_n}{1-\lambda} = \frac{-\lambda}{1-\lambda}x_n + \frac{1}{1-\lambda}x_{n+1} \end{equation}

これは点\(x_n\)と点\(x_{n+1}\)の内分点, または外分点を表します. 内分・外分の比は\(\frac{1}{1-\lambda}:\frac{-\lambda}{1-\lambda}\)です. 次に, \(\lambda\)の符号で場合分けして考えます. それぞれの場合について, 新しい数列の想定される挙動を考察します.

まず, \(\lambda > 0\)の場合を考えます. このとき, 元の数列\(\{ x_n\}_{n=0}^{\infty}\)は極限に(ほとんど)単調に近づきます. また, \(\lambda\)の絶対値が大きいほど元の数列の収束は遅いでしょう. 外分点\(x_n^{(1)} \)は\(x_n\)からは遠く, \(x_{n+1}\)に近い点です. また, \(\lambda\)の絶対値が大きいほど, \(x_{n+1}\)からはより遠ざかります. したがって, 新しい数列\(\{ x_n^{(1)}\}_{n=0}^{\infty}\)は, 元の数列の極限\(x^*\)により近づきます.

逆に\(\lambda < 0\)の場合を考えます. このとき, 数列\(\{ x_n\}_{n=0}^{\infty}\)は極限に振動しながら近づきます. また, \(\lambda\)の絶対値が大きいほど元の数列の収束は遅いでしょう. 内分点\(x_n^{(1)} \)は\(x_n\)と\(x_{n+1}\)の間の点です. また, \(\lambda\)の絶対値が大きいほど, \(x_{n}\)に近くなり, \(x_{n+1}\)から遠ざかります. したがって, 極限により近づきます.

実験

Richardosnの加速法を適用して, 実験してみます. ここでは例を2つほど示します. 元の数列を\(\{ x_n^{(0)}\}_{n=0}^{\infty}\), \(\{ x_n^{(m-1)}\}_{n=0}^{\infty}\)を加速させた数列を\(\{ x_n^{(m)}\}_{n=0}^{\infty}\)で表します. 今回の実験では, 4段階加速させました.

例1


まず, 1つ目の例です. 以下の数列を考えます.

\begin{equation} x_n = \frac{1}{n}\left( \frac{1}{2} \right)^n,\quad x_0=1 \end{equation}

この数列の極限は以下の通りです.

\begin{equation} \underset{n\to \infty}{\lim}x_n = 0 \end{equation}

また, \(\lambda\)の値を求めると, 以下のようになります.

\begin{equation} \lambda = \underset{n\to \infty}{\lim}\frac{x_{n+1}}{x_n} = \frac{1}{2} \end{equation}

下図に, Richardosnの加速法を適用した誤差を対数スケールでプロットしました. 狙い通り, 加速させるほど極限値への収束が早くなります.

【コードの実行結果】

例2


次に, 2つ目の例です. この例では, \(|\lambda|=1\)となります. 先ほどまでは\(|\lambda|< 1\)の場合を考えていました. =が成り立つ境界的な状況でも加速できるか実験してみます. 以下の数列を考えます.

\begin{equation} x_{n+1} = x_n - \frac{1}{2}\left( x_n^2-4\right),\quad x_0=1 \end{equation}

この数列の極限は以下の通りです.

\begin{equation} \underset{n\to \infty}{\lim}x_n = 2 \end{equation}

また, \(\lambda\)の値を求めると, 以下のようになります.

\begin{equation} \lambda = \underset{n\to \infty}{\lim}\frac{x_{n+1}-2}{x_n-2} = -1 \end{equation}

下図に, Richardosnの加速法を適用した誤差をプロットしました. 元の数列の収束はかなり遅いです. どうやら加速もうまくいってそうです. 他の境界的な例でもうまくいくのか気になりますね...(この例では偶然うまくいっているだけかも).

【コードの実行結果】

Aitkenの加速法

ここでは, 2つ目の加速法, Aitkenの加速法について基本的なアイデアを紹介し, 実験します.

アイデア

前節のRichardsonの加速法では, \(\lambda\)の値が既知であることが前提です. やってみると分かるのですが, 求めるのが意外と大変です. Aitkenの加速法は\(\lambda\)を推定し, 似たようなことを行います. 前節の\(\lambda\)の推定値として, 次の\(\hat{\lambda}_n\)を用います.

\begin{equation} \hat{\lambda}_n = \frac{x_{n+1}-x_{n}}{x_{n}-x_{n-1}} \end{equation}

実際, 次のように変形すると推定値として妥当であることが分かります.

\begin{equation} \hat{\lambda}_n = \frac{(x^* + C\lambda^{n+1})-(x^* + C\lambda^{n})}{(x^* + C\lambda^{n})-(x^* + C\lambda^{n-1})} = \lambda \end{equation}

これを推定値として, 次のように新たな数列を構成します.

\begin{equation} x_{n+1}^{(1)} = \frac{x_{n+1}-\hat{\lambda}_nx_n}{1-\hat{\lambda}_n} \end{equation}

この式自体は, Richardsonの加速法とほぼ同じです. さらに変形すると, 次のようになります.

\begin{equation} x_{n+1}^{(1)} = x_{n+1} - \frac{(x_{n+1}-x_{n})^2}{x_{n+1}-2x_{n}+x_{n-1}} \end{equation}

これで, \(\lambda\)の情報がなくても加速させた数列が構成できます. また, 何段階にも加速できます.

実験

Aitkenの加速法を適用して, 実験してみます. ここでは例を3つほど示します. 元の数列を\(\{ x_n^{(0)}\}_{n=0}^{\infty}\), \(\{ x_n^{(m-1)}\}_{n=0}^{\infty}\)を加速させた数列を\(\{ x_n^{(m)}\}_{n=0}^{\infty}\)で表します. 今回の実験では, 4段階加速させました.

例3


以下の数列を考えます. 前節の1つ目の例と同じ数列です.

\begin{equation} x_n = \frac{1}{n}\left( \frac{1}{2} \right)^n,\quad x_0=1 \end{equation}

下図に, Aitkenの加速法を適用した誤差をプロットしました.

【コードの実行結果】

例4


以下の数列を考えます.

\begin{equation} x_n = \left( 1+ \frac{1}{n} \right)^n,\quad x_0=1 \end{equation}

自然対数の底の定義式なので, 極限は以下の通りです.

\begin{equation} \underset{n\to \infty}{\lim}x_n = e \end{equation}

下図に, Aitkenの加速法を適用した誤差をプロットしました.

【コードの実行結果】

例5


最後の例です. ここでは, 以下の部分和をもつ無限級数の収束を加速させます.

\begin{equation} S_n = \sum_{k=1}^{n} \frac{1}{k^2} \end{equation}

よく知られている通り, この級数の極限は以下の通りです.

\begin{equation} \underset{n\to \infty}{\lim}S_n = \frac{\pi^2}{6} \end{equation}

下図に, Aitkenの加速法を適用した誤差をプロットしました.

【コードの実行結果】

コード

【Juliaコード1; インポート】
using Plots
pyplot()
【Juliaコード2; 関数の定義】
#数列を計算する関数
function calc_xs(func, nmax, x₀)
    xs = zeros(nmax)
    xs[1] = x₀
    for n in 2:nmax
        xs[n] = func(n, xs[n-1])
    end
    return xs
end

#Richardsonの加速法
function updateR(a,b,λ)
    return (a-λ*b)/(1-λ)
end

#Richardsonの加速法
function Richardson(xs, λ)
    N = length(xs)
    xs_new = zeros(N-1)
    for n in 2:N
        xs_new[n-1] = updateR(xs[n], xs[n-1], λ)
    end
    return xs_new
end

#Aitkenの加速法
function updateA(a,b,c)
    return a-(a-b)^2/(a-2*b+c)
end

#Aitkenの加速法
function Aitken(xs)
    N = length(xs)
    xs_new = zeros(N-2)
    for n in 3:N
        xs_new[n-2] = updateA(xs[n], xs[n-1], xs[n-2])
    end
    return xs_new
end
【Juliaコード3; Richardosn実験1】
#数列の計算
nmax = 30
λ = 1/2
xs0 = calc_xs(((n,x)->1/n/2^n), nmax, 1)
xs1 = Richardson(xs0, λ)
xs2 = Richardson(xs1, λ)
xs3 = Richardson(xs2, λ)
xs4 = Richardson(xs3, λ)

#誤差の計算
limx = 0.0
err0 = abs.(limx.-xs0)
err1 = abs.(limx.-xs1)
err2 = abs.(limx.-xs2)
err3 = abs.(limx.-xs3)
err4 = abs.(limx.-xs4)

#プロット
fig1 = plot(title="Richardson", legend=:bottomleft, xticks=1:nmax, yscale=:log10, xlabel="n", ylabel="log10|err|")
plot!(1:nmax, err0, label="x^0", marker=:auto)
plot!(2:nmax, err1, label="x^1", marker=:auto)
plot!(3:nmax, err2, label="x^2", marker=:auto)
plot!(4:nmax, err3, label="x^3", marker=:auto)
plot!(5:nmax, err4, label="x^4", marker=:auto)
savefig(fig1, "figs-accel/fig1.png")
【Juliaコード4; Richardosn実験2】
#数列の計算
nmax = 30
λ = -1
xs0 = calc_xs(((n,x)->x-(x^2-4)/2), nmax, 1)
xs1 = Richardson(xs0, λ)
xs2 = Richardson(xs1, λ)
xs3 = Richardson(xs2, λ)
xs4 = Richardson(xs3, λ)

#誤差の計算
limx = 2.0
err0 = abs.(limx.-xs0)
err1 = abs.(limx.-xs1)
err2 = abs.(limx.-xs2)
err3 = abs.(limx.-xs3)
err4 = abs.(limx.-xs4)

#プロット
fig2 = plot(title="Richardson", legend=:bottomleft, xticks=1:nmax, yscale=:log10, xlabel="n", ylabel="log10|err|")
plot!(1:nmax, err0, label="x^0", marker=:auto)
plot!(2:nmax, err1, label="x^1", marker=:auto)
plot!(3:nmax, err2, label="x^2", marker=:auto)
plot!(4:nmax, err3, label="x^3", marker=:auto)
plot!(5:nmax, err4, label="x^4", marker=:auto)
savefig(fig2, "figs-accel/fig2.png")
【Juliaコード5; Aitken実験1】
#数列の計算
nmax = 30
xs0 = calc_xs(((n,x)->1/n/2^n), nmax, 1)
xs1 = Aitken(xs0)
xs2 = Aitken(xs1)
xs3 = Aitken(xs2)
xs4 = Aitken(xs3)

#誤差の計算
limx = 0.0
err0 = abs.(limx.-xs0)
err1 = abs.(limx.-xs1)
err2 = abs.(limx.-xs2)
err3 = abs.(limx.-xs3)
err4 = abs.(limx.-xs4)

#プロット
fig3 = plot(title="Aitken", legend=:bottomleft, xticks=1:nmax, yscale=:log10, xlabel="n", ylabel="|err|")
plot!(1:nmax, err0, label="x^0", marker=:auto)
plot!(3:nmax, err1, label="x^1", marker=:auto)
plot!(5:nmax, err2, label="x^2", marker=:auto)
plot!(7:nmax, err3, label="x^3", marker=:auto)
plot!(9:nmax, err4, label="x^4", marker=:auto)
savefig(fig3, "figs-accel/fig3.png")
【Juliaコード6; Aitken実験2】
#数列の計算
nmax = 30
xs0 = calc_xs((n,x)->((1+1/n)^n), nmax, 1)
xs1 = Aitken(xs0)
xs2 = Aitken(xs1)
xs3 = Aitken(xs2)
xs4 = Aitken(xs3)

#誤差の計算
limx = exp(1)
err0 = abs.(limx.-xs0)
err1 = abs.(limx.-xs1)
err2 = abs.(limx.-xs2)
err3 = abs.(limx.-xs3)
err4 = abs.(limx.-xs4)

#プロット
fig4 = plot(title="Aitken", legend=:bottomleft, xticks=1:nmax, yscale=:log10, xlabel="n", ylabel="|err|")
plot!(1:nmax, err0, label="x^0", marker=:auto)
plot!(3:nmax, err1, label="x^1", marker=:auto)
plot!(5:nmax, err2, label="x^2", marker=:auto)
plot!(7:nmax, err3, label="x^3", marker=:auto)
plot!(9:nmax, err4, label="x^4", marker=:auto)
savefig(fig4, "figs-accel/fig4.png")
【Juliaコード7; Aitken実験3】
#数列の計算
nmax = 30
xs0 = calc_xs((n,x)->(x+1/n^2), nmax, 1)
xs1 = Aitken(xs0)
xs2 = Aitken(xs1)
xs3 = Aitken(xs2)
xs4 = Aitken(xs3)

#誤差の計算
limx = π^2/6
err0 = abs.(limx.-xs0)
err1 = abs.(limx.-xs1)
err2 = abs.(limx.-xs2)
err3 = abs.(limx.-xs3)
err4 = abs.(limx.-xs4)

#プロット
fig5 = plot(title="Aitken", legend=:bottomleft, xticks=1:nmax, yscale=:log10, xlabel="n", ylabel="|err|")
plot!(1:nmax, err0, label="x^0", marker=:auto)
plot!(3:nmax, err1, label="x^1", marker=:auto)
plot!(5:nmax, err2, label="x^2", marker=:auto)
plot!(7:nmax, err3, label="x^3", marker=:auto)
plot!(9:nmax, err4, label="x^4", marker=:auto)
savefig(fig5, "figs-accel/fig5.png")
参考文献

      [1]杉原正顕, 室田一雄, 数値計算法の数理, 岩波書店, 1998
      [2]有本卓, 数値解析(1), コロナ社, 1981
      [3]藤野清次, 数値計算の基礎―数値解法を中心に, サイエンス社, 1997
      [4]伊理正夫, 藤野和建, 数値計算の常識, 共立出版, 2019