equal_l2’s blog

※記載されている内容の正確性は保証しませんが、間違いを指摘していただければ直します。

離散フーリエ変換への道 (4) 離散フーリエ変換

前:
equal-l2-works.hatenablog.com

離散時間フーリエ変換では、逆変換が離散化されていない。

f[n]を周期Nの周期列とする。
これを離散時間フーリエ変換によって周期2\piで値が\frac{2\pi}{N}ごとに存在する離散関数に写す。

離散関数をデルタ関数の和で連続関数F(\omega)として表す。
\displaystyle F(\omega) = \sum^{\infty}_{n=-\infty} c_n \delta(\omega-\frac{2\pi n}{N})
(c_nは適当な係数)

この式を離散時間フーリエ逆変換の式
\displaystyle f[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} F(\omega) e^{i \omega n} \mathrm{d}\omega
に代入すると、
\begin{align}
f[n] 
&= \frac{1}{2\pi} \int_{-\pi}^{\pi} F(\omega) e^{i \omega n} \mathrm{d}\omega \\
&= \frac{1}{2\pi} \int_{-\pi}^{\pi} \left[ \sum^{\infty}_{k=-\infty} c_k \delta(\omega-\frac{2\pi k}{N}) \right] e^{i \omega n} \mathrm{d}\omega \\
&= \frac{1}{2\pi}  \sum^{N/2}_{k=-N/2} c_k e^{i \frac{2\pi k}{N} n} \\
&= \frac{1}{2\pi}  \sum^{N-1}_{k=0} c_k e^{i \frac{2\pi k}{N} kn}
\end{align}

c_k=\frac{2\pi}{N}F[k]となるような数列F[k]を導入して、
\displaystyle f[n] = \frac{1}{N}  \sum^{N-1}_{k=0} F[k] e^{i \frac{2\pi}{N} kn}
これが離散フーリエ逆変換である。

これに対応するように離散時間フーリエ変換の式を変えてやると、
\displaystyle F[k] = \sum^{N-1}_{n=0} f[n] e^{-i \frac{2\pi}{N} kn}
これが離散フーリエ変換である。

離散フーリエ変換への道 (3) 離散時間フーリエ変換

前:
equal-l2-works.hatenablog.com

離散時間フーリエ変換

離散集合f[n]を考える。


ここで、f[n]の各要素は等間隔に増大する変数x[n]に対応するものとする。
この対応をf(x[n])=f[n]とする。
ただし、x[0]=0とする。

今、変数x[n]-x[n-1] = \Delta xとすると、
x[n]=n \Delta x
と書ける。


離散関数f(x[n])を可積分な連続関数f(x)として扱うには、各サンプルにデルタ関数を掛けて
 \displaystyle f(x) = \sum^{\infty}_{n=-\infty} f[n] \delta(x-x[n])

あとはこの関数をフーリエ変換して、

\begin{align}
F(\omega) 
&= \int_{-\infty}^{\infty} f(y) e^{-i \omega y} \ \mathrm{d}y \\
&= \int_{-\infty}^{\infty} \left[ \sum^{\infty}_{n=-\infty} f[n] \delta(x-x[n]) \right] e^{-i \omega y} \ \mathrm{d}y \\
&= \sum^{\infty}_{n=-\infty} f[n] e^{-i \omega x[n]} \\
&= \sum^{\infty}_{n=-\infty} f[n] e^{-i \omega n \Delta x}
\end{align}


ここで\omega \Delta x\omegaとして再定義すると、
\displaystyle F(\omega) = \sum^{\infty}_{n=-\infty} f[n] e^{-i \omega n}
この変換を離散時間フーリエ変換という。

離散時間フーリエ逆変換

F(\omega) = \displaystyle \mathcal{F}\left[\sum^{\infty}_{n=-\infty} f[n] \delta(x-x[n])\right]
であることから、この両辺をフーリエ逆変換すれば
\displaystyle
\sum^{\infty}_{n=-\infty} f[n] \delta(x-x[n]) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega x} \mathrm{d}\omega
となる。

この式から元の離散集合f[n]を取り出すには、x = x[n]となるような積分で集合を作ればよい。
ただし、F(\omega)は項e^{-i \omega n}を含むことから\omegaについて周期2\piの周期性を持つので、積分範囲の大きさは2\piでよい。
ゆえに、離散時間フーリエ逆変換の式は、
\displaystyle f[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} F(\omega) e^{i \omega n} \mathrm{d}\omega
となる。

離散フーリエ変換への道 (2) フーリエ変換

前:
equal-l2-works.hatenablog.com

フーリエ変換

周期2Lの関数に対する複素フーリエ級数展開の式は次のようなものであった。
\displaystyle f(x)=\sum_{n=-\infty}^{\infty} c_n e^{in(\pi/L)x}
\displaystyle c_n = \frac{1}{2L} \int_{-L}^{L} f(y) e^{-in(\pi/L)y} \ \mathrm{d}y


係数を級数の式へ代入すると、

\displaystyle f(x)=\sum_{n=-\infty}^{\infty} \left[ \frac{1}{2L} \int_{-L}^{L} f(y) e^{-in(\pi/L)y} \ \mathrm{d}y \right] e^{in(\pi/L)x}


これを非周期関数に適用できるようにするには、L \rightarrow \inftyの極限をとればよい。
 \omega_n = n(\pi/L), \Delta\omega = \omega_n - \omega_{n-1} = \pi/Lを導入すると、

\displaystyle f(x)=\sum_{n=-\infty}^{\infty} \left[ \frac{1}{2\pi} \int_{-L}^{L} f(y) e^{-i \omega_n y} \ \mathrm{d}y \right] e^{i \omega_n x} \Delta\omega

いま、リーマン積分の定義
\displaystyle
\lim_{\Delta\omega \rightarrow 0} \sum_{n=-\infty}^{\infty} g(\omega_n) \Delta \omega = \int_{-\infty}^{\infty} g(\omega) \mathrm{d} \omega
を適用して和を積分にしてやると、


\displaystyle
f(x) = \int_{-\infty}^{\infty} \left[ \frac{1}{2\pi} \int_{-\infty}^{\infty} f(y) e^{-i \omega y} \ \mathrm{d}y \right] e^{i \omega x} \mathrm{d}\omega


これは、次のような2本の式として書き直せる。

\begin{align}
F(\omega) &= \int_{-\infty}^{\infty} f(y) e^{-i \omega y} \ \mathrm{d}y \\
f(x) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega x} \mathrm{d}\omega
\end{align}
F(\omega)f(x)フーリエ変換という。
f(x)F(\omega)の逆フーリエ変換という。
フーリエ変換を行う写像\mathcal{F}、逆フーリエ変換を行う写像(\mathcal{F}の逆写像)を\mathcal{F}^{-1}と表す。
すなわち、\mathcal{F}\left[f(x)\right] = F(\omega),\mathcal{F}^{-1}\left[F(\omega)\right] = f(x)

離散フーリエ変換への道 (1) フーリエ級数展開

離散フーリエ変換を勉強することにしたので、勉強したことをメモしておく。

フーリエ級数展開とは

任意の周期関数を\sin関数と\cos関数の組み合わせで展開する方法。

周期2\piの関数のフーリエ級数展開

周期2\piの関数f(x)は次のようにフーリエ級数展開できる。

\displaystyle
f(x)=a_0+\sum_{n=1}^{\infty}(a_n\cos{nx}+b_n\sin{nx})

係数a_n,b_n \in \mathbb{R}は次のように求められる。

\begin{align}
a_n &= \frac{1}{\pi} \int_{-\pi}^{\pi} f(y)\cos{ny} \ \mathrm{d}y \\
b_n &= \frac{1}{\pi} \int_{-\pi}^{\pi} f(y)\sin{ny} \ \mathrm{d}y
\end{align}

本によっては積分区間が[0,2\pi]のものがあるが、積分区間の大きさが2\piであれば何でも良いらしい。
ネットで証明を見つけたので載せておく。

積分な周期2\piの関数g(x)について\int_{0}^{2\pi} g(x) \ \mathrm{d}x = \int_{a}^{a+2\pi} g(x) \ \mathrm{d}xを示す。

 a+b=2k\piとなるようなk \in \mathbb{Z}, \ b \in \{x\in\mathbb{R} | 0 \leq x \lt 2\pi\}を考えると、
 
\begin{align}
\int_{a}^{a+2\pi} g(x) \ \mathrm{d}x
&= \int_{a}^{a+b} g(x) \ \mathrm{d}x + \int_{a+b}^{a+2\pi} g(x) \ \mathrm{d}x \\
&= \int_{a}^{2k\pi} g(x) \ \mathrm{d}x + \int_{2k\pi}^{a+2\pi} g(x) \ \mathrm{d}x \\
&= \int_{a}^{2k\pi} g(x) \ \mathrm{d}x + \int_{2k\pi-2\pi}^{a} g(x) \ \mathrm{d}x \\
&= \int_{2k\pi-2\pi}^{2k\pi} g(x) \ \mathrm{d}x\\
&= \int_{0}^{2\pi} g(x) \ \mathrm{d}x\\
\end{align}
よって、周期2\piの関数g(x)の定積分積分区間の大きさが2\piであれば同じである。

周期2\piの関数f(x)\sin関数や\cos関数の積は、周期2\piの関数である。
よって、フーリエ係数は積分区間が2\piであれば積分区間によって変化しない。

周期2L \ (L \lt 0)の関数のフーリエ級数展開

周期2Lの関数f(x)は次のようにフーリエ級数展開できる。

\displaystyle
f(x)=a_0+\sum_{n=1}^{\infty}(a_n\cos{\frac{n\pi x}{L}}+b_n\sin{\frac{n\pi x}{L}})

係数a_n,b_n \in \mathbb{R}は次のように求められる。

\begin{align}
a_n &= \frac{1}{L} \int_{-L}^{L} f(y)\cos{\frac{n\pi y}{L}} \ \mathrm{d}y \\
b_n &= \frac{1}{L} \int_{-L}^{L} f(y)\sin{\frac{n\pi y}{L}} \ \mathrm{d}y
\end{align}

これらの式は、h(t)=f(Lt/\pi)なる周期2\piの関数のフーリエ級数展開を考えてやると導出できる。

複素フーリエ級数展開

今まで扱ってきたフーリエ級数展開に、オイラーの公式やド・モアブルの公式などをうまく作用させてやると、違う形のフーリエ級数展開が得られる。

周期2\piの関数f(x)は次のようにもフーリエ級数展開できる。
\displaystyle f(x)=\sum_{n=-\infty}^{\infty} c_n e^{inx}
係数c_n \in \mathbb{C}は次のように求められる。
\displaystyle c_n = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(y) e^{-iny} \ \mathrm{d}y

周期2Lの場合も同様に、
\displaystyle f(x)=\sum_{n=-\infty}^{\infty} c_n e^{in(\pi/L)x}
\displaystyle c_n = \frac{1}{2L} \int_{-L}^{L} f(y) e^{-in(\pi/L)y} \ \mathrm{d}y

広義Runge–Kutta法の数値安定性

数値解法の安定性とは

(大雑把に言えば)どのくらい刻み幅が大きくなると計算が明らかにおかしくなるか、ということである。

A-stabilityとは

数値解法の安定性を評価する指標の一つである。

次の常微分方程式を考える。
\displaystyle\frac{\mathrm{d} x(t)}{\mathrm{d} t}=kx(t) \ (k \in \mathbb{C}) \cdots (1)
これをx(0)=1として解析的に解けば
\displaystyle x(t)=e^{kt}

このとき、\mathrm{Re}(k)<0ならば
\displaystyle\lim_{t \to \infty} x(t) = 0 \cdots (2)
である。

ある数値解法で(1)式を解いたとき、刻み幅hにかかわらず式(2)を満たすならば、その解法はA-stableであるという。

安定領域とは

上の係数kと、刻み幅hの積hk\lim_{t \to \infty} x(t) = 0を満たすhkの領域である。
広義Runge–Kutta法においてはこの領域は容易に求められる。

広義Runge–Kutta法において次の式が成り立つ。
\displaystyle x_{n+1}=\phi(hk) \cdot x_{n}
ここで、x_{n}=x(hn)+x_0である。x_0=x(t_0)であり、t_0tの初期値である。
また、\phiは各解法で固有な関数で、安定関数とよばれる。
安定領域は、\{z \in \mathbb{C}\ |\ |\phi(z)| < 1\}なる領域である。

ある解法がA-stableであるには、安定領域が\{z \in \mathbb{C}\ |\ \mathbb{Re}(z) < 0\}なる領域を含んでいればよい。

各解法の安定領域

今までこのブログで書いてきた各解法(すべて広義Runge–Kutta法である)について、その安定領域を調べる。

陽的Euler法

陽的オイラー法の式は
\displaystyle x_{n+1} = x_n + f(t_n,x_n)h

これに、f(t_n,x_n) = kx_nを代入して、
\begin{eqnarray}
		x_{n+1}
		&=& x_n + kx_n h
		&=& (1+kh)x_n
	\end{eqnarray}

よって、安定関数は
\displaystyle\phi(z)=1+z
ゆえに、安定領域は
\displaystyle|1+z|<1
(以後の安定関数の計算も全て同様の方法である)

f:id:equal_l2:20150922093556p:plain

(縦軸:虚軸(\mathrm{Im}) 横軸:実軸(\mathrm{Re}) 以後の画像も同様)

陰的Euler法

陰的オイラー法の式は
\displaystyle x_{n+1} = x_n + f(t_{n+1},x_{n+1})h

安定関数は
\displaystyle\phi(z)=\frac{1}{1-z}
ゆえに、安定領域は
\displaystyle\left|\frac{1}{1-z}\right|<1
すなわち
\displaystyle|1-z|>1

f:id:equal_l2:20150922093612p:plain

この解法はA-stableである。

台形法

台形法の式は
\displaystyle x_{n+1} = x_n + \frac{h}{2}\left( f(t_n,x_n) + f(t_{n+1},x_{n+1}) \right)
安定関数は
\displaystyle\phi(z)=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}
ゆえに、安定領域は
\displaystyle\left|\frac{1+\frac{z}{2}}{1-\frac{z}{2}}\right|<1
変形して、
\displaystyle\mathrm{Re}(z)<0

f:id:equal_l2:20150922093624p:plain

この解法はA-stableである。

陽的中点法

陽的中点法の式は
\displaystyle x_{n+1} = x_{n} + f(t_n+\frac12 h,x_n + \frac12 h f(t_n,x_n))h
安定関数は
\displaystyle\phi(z)=\frac12 z^2 + z + 1

f:id:equal_l2:20150922093634p:plain

陰的中点法

陰的中点法の式は
\displaystyle x_{n+1} = x_{n} + f(t_n+\frac12 h,\frac12 \left( x_n + x_{n+1} \right))h
安定関数は
\displaystyle\phi(z)=\frac{1+\frac{z}{2}}{1-\frac{z}{2}}
これは台形法と同じ安定関数なので、安定領域は
\displaystyle\mathrm{Re}(z)<0
(図省略)

陽的化した陰的Euler法

陽的化した陰的オイラー法の式は
\displaystyle x_{n+1} = x_n + f(t_{n+1}, x_n + f(t_n,x_n)h)h
安定関数は
\displaystyle \phi(z)= z^2 + z + 1

f:id:equal_l2:20170302060531p:plain

陽的化した台形法 (Heun法)

陽的化した台形法の式は
\displaystyle x_{n+1} = x_n + \frac{h}{2}\left( f(t_n,x_n) + f(t_{n+1},x_n + f(t_n,x_n)h) \right)
安定関数は
\displaystyle\phi(z)=\frac12 z^2 + z + 1
これは陽的中点法と同じ安定関数である。

その他

一般に陰的解法はA-stableである。
また、広義Runge–Kutta法は次数が上がるほど安定領域が広がる傾向があるようだ。

最小二乗法(8) 重み付き多項式フィッティング

前:最小二乗法(7) 項の欠けた線形近似 - equal_l2’s blog
次:なし
\def\fxb{f(x_i,\vec\beta)}\def\csq{\chi^2}\def\ssq{{\sigma_i}^2}
最小二乗法(1)~(7)で行ってきたフィッティングは、各y_iがすべて同じ誤差\sigma_yを持つものとして行ってきた。
しかし、時には各y_iにそれぞれ違う誤差\sigma_iを付けたい気分の日もあるだろう。

このとき、最小にするのはもはや残差の二乗和ではなく、
\displaystyle\csq=\sum_{i=1}^{n} \frac{\left(y_i- \fxb\right)^2}{\ssq}
である。

方法論は同じで、全てのパラメータ\beta_jについて、
{\displaystyle \frac{\partial \csq}{\partial \beta_j}=0 \rm \hspace{10pt} for \hspace{2pt} all \hspace{2pt} \it j}
となるような\beta_jを求めればよい。

実際に偏微分すると、
\begin{eqnarray}
\frac{\partial \csq}{\partial \beta_j}
&=& \frac{\partial}{\partial \beta_j}\sum_{i=1}^{n}\frac{\left(y_i- \fxb\right)^2}{\ssq}\\
&=& -2 \sum_{i=1}^{n}\frac{\partial \fxb}{\partial \beta_j}\frac{y_i- \fxb}{\ssq}= 0
\end{eqnarray}

変形すれば、正規方程式は
{\displaystyle
\sum_{i=1}^{n}\frac{\partial \fxb}{\partial \beta_j} \frac{y_i}{\ssq} =
\sum_{i=1}^{n}\frac{\partial \fxb}{\partial \beta_j} \frac{\fxb}{\ssq}
}
となる。

各パラメータの誤差は、誤差の伝播式より、
\displaystyle
\sigma(\beta_j)=\sqrt{\sum_{i=1}^{n}\left\{ {\sigma_i}^2 \left( \frac{\partial \beta_j}{\partial y_i} \right)^2\right\}}
となる。

ここで、モデル関数として多項式
{\displaystyle \fxb=\sum_{k=0}^{d} \beta_k {x_i}^k}
を用いて、
\displaystyle X_m = \sum_{i=1}^{n}\frac{{x_i}^m}{\ssq},Y_m=\sum_{i=1}^{n}\frac{{x_i}^m y_i}{\ssq}
とおけば、解くべき正規方程式は

    \left(\begin{array}{c}{Y_0 \\ Y_1 \\ \vdots \\ Y_d }\end{array}\right) =
    
    \left(
    \begin{array}{cccc}
      X_0& X_1& \cdots& X_d& \\
      X_1& X_2& \cdots& X_{1+d}& \\
      \vdots& \vdots& \ddots& \vdots& \\
      X_d& X_{d+1}& \cdots& X_{d+d}& \\
    \end{array}
    \right)

    \left(\begin{array}{c}{\beta_0 \\ \beta_1 \\ \vdots \\ \beta_d }\end{array}\right)
となり、(項の中身が変わったことを除けば)重みをすべて同じにした場合と変わらない。

次回は、一応d=1についての結果を述べておこう。

Qiitaに登録してみた

C++でちょっと気付いたこととかをブログに書くのはなんだか気が引けるので、Qiitaに登録してみた。
EqualL2 - Qiita
Tipsみたいなものはこっちに書いていこうと思う。

シリーズ物は引き続きこっちで。

常微分方程式の数値解法(3) 中点法

今回も引き続き次のような式の解法を考える。
\def\disp{\displaystyle}\disp\frac{dx}{dt} = f(t,x)

中点法の準備

Euler法の導出と同様のことを、中央差分で行うと、
x_{n+\frac12} = x_{n-\frac12} + f(t_n,x_n)h
添字からマイナスの項を消せば、
\underline{x_{n+1} = x_{n} + f(t_{n+\frac12},x_{n+\frac12})h}\cdots(1)

この式は代数的に処理してx_{n+1}, x_{n}の式にすることはできない。
x_{n+\frac12}を、近似などでx_{n+1}, x_{n}で表される式にしなければならない。

陽的中点法

x_{n+\frac12}テイラー展開を用いれば、
\begin{align}
x_{n+\frac12} &= x(t_n +\frac12 h) \\
&\approx x(t_n) + \frac12 h \frac{dx}{dt}(t_n) \\
&= x(t_n) + \frac12 h f(t_n,x_n)
\end{align}

この結果を上の(1)式へ代入すれば、
\underline{x_{n+1} = x_{n} + f(t_n+\frac12 h,x_n + \frac12 h f(t_n,x_n))h}

陰的中点法

x_{n+\frac12} \approx \frac12 \left( x_n + x_{n+1} \right)
とする。これは、t_n ~ t_{n+1}の間を1次関数であるものとして近似したと考えられる。

これを(1)式へ代入すると、
\underline{x_{n+1} = x_{n} + f(t_n+\frac12 h,\frac12 \left( x_n + x_{n+1} \right))h}

ちなみに、陽的オイラー法で右辺のx_{n+1}を展開すると、陰的中点法は陽的中点法に一致する。

図的イメージ

高さf(t_{n+\frac12},x_{n+\frac12})、幅hの矩形の面積を足し合わせる。
f:id:equal_l2:20150524192345p:plain
誤差がいい具合に消しあうことが多いらしく、精度はオイラー法よりも良い傾向がある。

常微分方程式の数値解法(2) Euler法

実際に1次の常微分方程式を解く。

基本

\def\disp{\displaystyle}
次のような1次の常微分方程式を考える。
\disp\frac{dx}{dt} = f(t,x)

このような常微分方程式に対して、ある種の数値解法の公式は漸化式の形で表される。
いま、次のような記法を導入する。
t_n = t_0 + nh
x_n = x(t_n)
ここで、nは一般には整数であるが、便宜上は非整数を使ってもよい。

これを用いて、各差分による微分近似を書き直すと、
\displaystyle\begin{eqnarray}
\frac{dx}{dt} 
&=& \frac{x_{n+1} - x_n}{h} \\
&=& \frac{x_n - x_{n-1}}{h} \\
&=& \frac{x_{n+\frac12}-x_{n-\frac12}}{h}
\end{eqnarray}
このとき、常微分方程式
\disp\frac{dx}{dt} = f(t_n,x_n)
となる。

このときのhを、積分の刻み幅とも呼ぶ。
hは小さいほうが微分近似が正確になるから、積分の時にも小さいほうがいい。

陽的Euler法

前進差分による微分の近似は、
\disp\frac{dx}{dt} = \frac{x_{n+1} - x_n}{h}
これと、
\disp\frac{dx}{dt} = f(t_n,x_n)
より、
\disp f(t_n,x_n) = \frac{x_{n+1} - x_n}{h}
ゆえに、
\underline{x_{n+1} = x_n + f(t_n,x_n)h}

陰的Euler法

後退差分で、上と同様のことを行うと、
x_{n} = x_{n-1} + f(t_n,x_n)h

マイナスがあるのが嫌なので、n \to n+1とすれば、
\underline{x_{n+1} = x_n + f(t_{n+1},x_{n+1})h}

この式は、x_{n+1}を求めるのに、fの引数としてx_{n+1}が使われている。
fxに陽に依存する場合は、実際にfの式を代入して、x_{n+1}についての式に整理しなければならない。

このように、実際に個々のfを代入しないとx_{n+1}に関する式が得られないものを陰解法 (implicit method)という。
対して、どのようなfについても同じようにx_{n+1}が得られるものを、陽解法(explicit method)という。

面倒な陰解法が何の役に立つのかというと、陽解法に比べて安定性が高いのである。
安定性が高いとは、刻み幅が大きいときでも解析解との乖離が小さいという意味である。
安定性については別の機会に詳しく説明する。

台形法

上の2つの方法の式を辺々足す。つまり、
\begin{array}{rl}
x_{n+1} & = x_n + f(t_n,x_n)h \\
\ +) x_{n+1} & = x_n + f(t_{n+1},x_{n+1})h \\ \hline
2 x_{n+1} & = 2 x_n + h \left( f(t_n,x_n) + f(t_{n+1},x_{n+1}) \right) 
\end{array}

辺々を2で割れば
\disp\underline{x_{n+1} = x_n + \frac{h}{2}\left( f(t_n,x_n) + f(t_{n+1},x_{n+1}) \right)}

これも陰解法である。他の2法に比べて精度が良い。

各解法の図的イメージ

陽的オイラー法は、高さf(t_n,x_n)、幅hの矩形の面積を足し合わせていくことで積分を近似する。
f:id:equal_l2:20150524011810p:plain

陰的オイラー法は、高さf(t_{n+1},x_{n+1})、幅hの矩形の面積を足し合わせていくことで積分を近似する。
f:id:equal_l2:20150524011825p:plain

台形法は、上底f(t_n,x_n)、下底f(t_{n+1},x_{n+1})、高さhの台形の面積を足し合わせていくことで積分を近似する。
f:id:equal_l2:20150524011818p:plain

台形法の精度の良さがイメージできただろうか。

陰解法を陽解法にする

陰解法は、陽的オイラー法を使って陽解法にすることができる。
右辺のx_{n+1}を陽的オイラー法による推定値\tilde{x}_{n+1} = x_n + f(t_n,x_n)hに置き換えれば良い。
この方法で陽解法にした台形法は、特にホイン法(Heun's method)と呼ばれる。

常微分方程式の数値解法(1) 差分

常微分方程式を数値的に解くためには、積分を離散的に扱う必要がある。

まず、その準備として関数の差分という概念を導入する。

差分は汎関数(のようなもの?)で、3種類ある。

前進差分 \Delta_h[x](t) =  x(t + h) - x(t)

後退差分 \nabla_h[x](t) =  x(t) - x(t-h)

中央差分 \delta_h[x](t) =  x(t+\frac12 h)-x(t-\frac12 h)

十分に小さいhについて、上の差分をhで割れば微分\frac{dx}{dt}の近似になる。
つまりは、
\displaystyle\begin{eqnarray}
\frac{dx}{dt} 
&=& \frac{x(t + h) - x(t)}{h} \\
&=& \frac{x(t) - x(t-h)}{h} \\
&=& \frac{x(t+\tfrac12h)-x(t-\tfrac12h)}{h}
\end{eqnarray}

前進差分、後退差分から得られる微分近似はhの1次の精度であるが、中央差分から得られるそれは2次の精度である。

ただし、振動する関数に対しては中央差分による微分近似が具合が悪くなることもあるようだ。