📂数値解析

数値積分

定義 ¹

$f : [a,b] \to \mathbb{R}$が$[a,b]$で積分可能であり、$[a,b]$を$a = x_{0} < \cdots < x_{n} = b$のようなノードポイントで分割したとする。

1. 積分オペレータ$I$を$\displaystyle I(f) := \int_{a}^{b} f(x) dx$のように定義する。
2. 積分オペレータ$I_{n}$を$\displaystyle I_{n} (f) := \sum_{k=1}^{n} \int_{x_{k-1}}^{x_{k}} f(x) dx$のように定義する。
3. エラー$E_{n}$を$E_{n} (f) := I (f) - I_{n} ( f )$のように定義する。
4. $\displaystyle \lim_{n \to \infty} {{\tilde{E}_{n} (f) } \over { E_{n} (f) }} = 1$を満たす$\tilde{E}_{n}$を$E_{n}$のアシンプトティックエラーという。

説明

定積分を計算することは数値解析の目標の一つであることに誰もが異論はないだろう。実際、区分求積法のアイデアを思い起こせば、空間を可能な限り細かく分割すればするほど、積分値をよく近似できるだろうと予測できる。だから、数値解析で興味を持つ数値的積分とは、具体的にその誤差がどれくらいか、閉区間でない場合はどのように積分すべきかなどであろう。

例示

台形則

数値積分の例として、$f \in C^2 [a,b]$に対する台形則を見てみよう。台形則は$I(f)$に対して$f$を区間$[a,b]$でリニアインターポレーションしてその関数の積分値として$I(f)$を近似する方法である。

$$ I_{1}^{1} (f) := \left( {{ b - a } \over { 2 }} \right) [ f(a) + f(b) ] $$ この場合、実際の$I(f)$と$I_{1}^{1} (f)$の誤差$E_{1}^{1} (f)$はある$\xi \in [a,b]$に対して以下のように計算される。

ポリノミアルインターポレーション:

• [4] 実際の関数との誤差：$(n+1)$回微分可能な$f : \mathbb{R} \to \mathbb{R}$とある$\xi \in \mathscr{H} \left\{ x_{0} , \cdots , x_{n} \right\}$に対して、$f$のポリノミアルインターポレーション$p_{n}$はある$t \in \mathbb{R}$に対して$\displaystyle f(t) - p_{n} (t) = {{ (t - x_{0}) \cdots (t - x_{n}) } \over { (n+1)! }} f^{(n+1)} ( \xi )$を満たす。

$$ \begin{align*} \displaystyle E_{1}^{1} (f) :=& I(f) - I_{1}^{1} (f) \\ =& \int_{a}^{b} \left[ f(x) - {{ f(b) ( x - a ) - f(a) (x - b) } \over { b - a }} \right] dx \\ =& \int_{a}^{b} \left[ f(x) - p_{1} (x) \right] dx \\ =& {{1} \over {2}} f '' ( \xi ) \int_{a}^{b} (x-a) (x-b) dx \\ =& \left[ {{1} \over {2}} f '' ( \xi ) \right] \left[ - {{1} \over {6}} (b-a)^{3} \right] \\ =& - {{1} \over {12}} (b-a)^{3} f '' ( \xi ) \end{align*} $$

合成台形則

今、$[a,b]$を$n$個の区間に分割し、各$[x_{k-1} , x_{k} ]$で台形則で積分を近似すると考えよう。この方法を合成台形則という。 $$ I_{n}^{1} (f) := \sum_{k=1}^{n} {{ h } \over { 2 }} \left( f (x_{k-1}) + f( x_{k} ) \right) $$ では、実際の$I(f)$と$I_{n}^{1} (f)$の誤差はある$\xi_{k} \in [x_{k-1}, x_{k} ]$に対して以下のように計算される。 $$ \begin{align*} \displaystyle E_{n}^{1} (f) =& I (f) - I_{n}^{1} (f) \\ =& \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ h^3 n } \over { 12 }} \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right) \end{align*} $$ ここで、$\displaystyle \left( {{ 1 } \over { n }} \sum_{k=1}^{n} f '' ( \xi_{k} ) \right)$は$f’’ ( \xi_{k} )$の平均なので、中間値の定理により、ある$\xi \in [a,b]$が存在して$\displaystyle \min_{x \in [a,b]} f ''(x) \le f ''( \xi ) \le \max_{x \in [a,b]} f ''(x)$を満たすべきである。したがって、 $$ E_{n}^{1} (f) = - {{ (b-a) h^2 } \over {12}} f '' (\xi ) $$ このように誤差を直接式で見ることで、積分区間$[a,b]$が広がるほど誤差が大きくなり、$h$が小さいほど正確になることがわかる。しかし、$f’’ ( \xi ) $のために具体的にどれほど間違っているかはわからない。これを克服するために以下のようにアシンプトティックエラーを求める。 $$ \begin{align*} \lim_{n \to \infty} {{ E_{n}^{1} (f) } \over { h^2 }} =& \lim_{n \to \infty} {{1} \over {h^2}} \sum_{k=1}^{n} \left( - {{ h^3 } \over { 12 }} f '' ( \xi_{k} ) \right) \\ =& - {{ 1 } \over { 12 }} \lim_{ n \to \infty} \sum_{k=1}^{n} h f '' ( \xi_{k} ) \\ =& - {{ 1 } \over { 12 }} \int_{a}^{b} f ''(x) dx \\ =& - {{ 1 } \over { 12 }} [ f '(b) - f '(a) ] \end{align*} $$ 誤差$E_{n}^{1} (f)$だけを知っていても定理によれば$f’’ ( \xi )$が存在することはわかるが、具体的な数値をつけることはできなかった。しかし、今、 $$ \tilde{E}_{n}^{1} (f) := \lim_{n \to \infty} E_{n}^{1} (f) $$ とした場合、十分に大きい$n$に対して $$ E_{n}^{1} (f) \approx - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ] = \tilde{E}_{n}^{1} (f) $$ であることがわかり、$f ' (a) , f '(b)$の計算は比較的簡単なので、比較的簡単かつ正確に誤差を把握することができる。しかも、ここまで来れば誤差がその程度であると終わる理由もない。計算をする前からこのように誤差が出ると予測できるならば、始めから誤差を補正することもできる。 $$ I(f) - I_{n}^{1} (f) \approx \tilde{E}_{n}^{1} (f) \implies I(f) \approx I_{n}^{1} (f) + \tilde{E}_{n}^{1} (f) $$

修正台形則

$I(f)$を計算するオペレータとして以下のように$CT_{n}$を定義しよう。 $$ CT_{n} (f) := h \left[ {{1} \over {2}} f (x_{0} ) + f (x_{1} ) + \cdots + f (x_{n-1} ) + {{1} \over {2}} f (x_{n} ) \right] - {{ h^2 } \over { 12 }} [ f '(b) - f '(a) ] $$

この方法を修正台形則と言う。[ 注記: 上で台形則が改良される過程を見ていれば、わざわざこれらを明確に区別する必要はない。単に台形則と言っても意味は通じる。 ]

実装

実際に$\displaystyle \int_{0}^{\pi} e^{x} \cos (x) dx = -{{e^{\pi} + 1} \over {2}} \approx - 12.0703$を$I_{n}^{1}$と$CT_{n}$で計算してみた結果は以下の通り。

最初の列は$n$を示し、2番目の列は$\left| I(f) - I_{n}^{1} (f) \right|$、3番目の列は$\left| I(f) - CT (f) \right|$を示している。合成台形則はノードを2倍に増やすたびに誤差が$\displaystyle {{1} \over {4}}$倍に減少し、修正台形則は$\displaystyle {{1} \over {16}}$倍に減少していることが確認できる。数値だけを見ても修正台形則がはるかに優れていることは明らかだが、特に速度が非常に速い点が良い。

以下はPythonで書かれた例示コードである。

import numpy as np
def f(x) :
    return np.exp(x)*np.cos(x)
 
def d(f,x,tol=10**(-8)) :
    h=1
    d1=1
    d2=0
 
    while(abs(d1-d2)>tol) :
        d1=d2
        d2=(f(x+h)-f(x))/h
        h=h/2
 
    return(d2)
 
def I(f,a,b,n) :
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f(x)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x))
 
def CT(f,a,b,n) :
    h=(b-a)/n
    x=np.linspace(a,b,n+1)
    x=f(x)
    x[0]=x[0]/2
    x[-1]=x[-1]/2
    return(h*sum(x)-(h**2)*(d(f,b)-d(f,a))/12)
 
 
TV=-(np.exp(np.pi)+1)/2
print("True Value : %2.4f" % TV)
 
print("|   n |  I(f) Err | CT(f) Err |")
print("="*30)
for n in range(1,10) :
    print("| %3d |  %8.2e |  %8.2e |" % (2**n,abs(TV-I(f,0,np.pi,2**n)),abs(TV-CT(f,0,np.pi,2**n))))