はじめに
この記事では，常微分方程式の精度保証付き数値計算のさわりを解説する．
精度保証付き数値計算 (validated numerics) は「数値計算によって得られた結果に対して，数学的に厳密な誤差限界を与える手法」である（荻田・柏木・劉，2018）．こう聞くと，精度保証付き数値計算は「従来の数値解法に誤差評価を付け加える理論」であるかのように思われるかもしれない．しかし，それは精度保証付き数値計算のほんの一面にすぎない．
誤差を評価するには，まず解の存在を証明しなければならない．計算の目的によっては，解が一意であることも示す必要があるだろう．こうした要件から，従来の数値解法をそのまま精度保証付き数値解法に修正できるとは限らず，しばしば精度保証ならではのアプローチが必要になる．この記事によって，この記事によって，精度保証付き数値計算が持つそうした特徴を少しでも伝えられたら嬉しい．
Googleドライブから，この記事のPDF版 とMarkdown版 をダウンロードできる．印刷して読みたい方にはPDF版をおすすめする．
浮動小数点演算
計算機において，実数の演算は浮動小数点演算 (floating-point arithmetic; FP) という，Q \mathbb{Q} Q IEEE 754-2008 という標準規格に則った浮動小数点演算を実装している．
この記事を読むには，浮動小数点演算について次のことを抑えていれば十分である．詳細を知りたければBoldo (2023) が参考になるだろう．
浮動小数点数の全体集合F ˉ \mathbb{\bar{F}} F ˉ Q \mathbb{Q} Q F \mathbb{F} F { ± ∞ } \lbrace\pm\infty\rbrace { ± ∞ }
F \mathbb{F} F ある自然数p p p F \mathbb{F} F p p p 2 ± 1 2^{\pm 1} 2 ± 1 p p p 精度 (precision) という．
浮動小数点数でない実数x x x F ˉ \mathbb{\bar{F}} F ˉ x x x 丸め (rounding) という．代表的な丸め関数は以下の3通りである．
x x x F ˉ \mathbb{\bar{F}} F ˉ RD ( x ) \operatorname{RD}(x) RD ( x ) 負の無限大への丸め (round towards negative infinity) という．x x x F ˉ \mathbb{\bar{F}} F ˉ RU ( x ) \operatorname{RU}(x) RU ( x ) 正の無限大への丸め (round towards positive infinity) という．x x x F ˉ \mathbb{\bar{F}} F ˉ RN ( x ) \operatorname{RN}(x) RN ( x ) 最近接丸め (round to nearest) という．
最近接丸め関数は一意ではない．IEEE 754-2008にしたがうシステムは偶数丸め (round ties to even) という最近接丸めを実装しており，指定しない限り偶数丸めが使われる．
区間演算
F \mathbb{F} F R \mathbb{R} R F ˉ \mathbb{\bar{F}} F ˉ
a ≤ x ≤ b ( a , b ∈ F ˉ ) a \leq x \leq b\quad(a,b\in\mathbb{\bar{F}}) a ≤ x ≤ b ( a , b ∈ F ˉ ) の形で表現する．x x x [ a , b ] \lbrack a,b\rbrack [ a , b ] 区間演算 (interval arithmetic; IA) に置き換えて実行する．
はじめから浮動小数点演算の誤差まで考慮して区間演算を議論するのは面倒なので，この記事ではR \mathbb{R} R F \mathbb{F} F
S S S R \mathbb{R} R x 1 , x 2 ∈ S x_{1},x_{2}\in S x 1 , x 2 ∈ S x 1 ≤ x 2 x_{1}\leq x_{2} x 1 ≤ x 2
x 1 ≤ u ≤ x 2 ⟹ u ∈ S x_{1}\leq u\leq x_{2}
\implies u\in S x 1 ≤ u ≤ x 2 ⟹ u ∈ S が成立するとき，S S S 閉区間 (closed interval) という．この定義のもとでは∅ \emptyset ∅ R \mathbb{R} R { a } \lbrace a\rbrace { a } a ∈ R a\in\mathbb{R} a ∈ R
閉区間S S S a a a b b b S S S [ a , b ] \lbrack a,b\rbrack [ a , b ]
[ a , b ] ≔ { x ∈ R ∣ a ≤ x ≤ b } \lbrack a,b\rbrack \coloneqq \lbrace x\in\mathbb{R}\mid a\leq x\leq b\rbrace [ a , b ] : = { x ∈ R ∣ a ≤ x ≤ b } とする．また，しばしば実数a a a { a } = [ a , a ] \lbrace a\rbrace=\lbrack a,a\rbrack { a } = [ a , a ] [ a , + ∞ ] \lbrack a,+\infty\rbrack [ a , + ∞ ] [ − ∞ , b ] \lbrack-\infty,b\rbrack [ − ∞ , b ] ± ∞ \pm\infty ± ∞
m × n m\times n m × n S S S S S S m × n m\times n m × n 区間行列 (interval matrix) という．m = 1 m=1 m = 1 n = 1 n=1 n = 1 区間ベクトル (interval vector) ともいう．
以下では閉区間を[ x ] \lbrack x\rbrack [ x ] x x x
[ x ] = [ x ▽ , x △ ] ( x ▽ = inf [ x ] and x △ = sup [ x ] ) \lbrack x\rbrack = \lbrack x^{\triangledown},x^{\vartriangle}\rbrack\quad(x^{\triangledown}=\inf\lbrack x\rbrack\;\mathrel{\textrm{and}}\;x^{\vartriangle}=\sup\lbrack x\rbrack) [ x ] = [ x ▽ , x △ ] ( x ▽ = inf [ x ] and x △ = sup [ x ]) とする．また，1次元閉区間の全体集合をI R \mathbb{IR} IR
I R n , I R m × n \mathbb{IR}^{n},
\quad\mathbb{IR}^{m\times n} IR n , IR m × n と書く．
区間行列は「各成分が閉区間である行列」とも「行列を端点とする閉区間」ともみなせる．そこで，次のように記号を定める．
m × n m\times n m × n [ X ] = { ( u i j ) ∣ a i j ≤ u i j ≤ b i j } \lbrack\bm{X}\rbrack=\lbrace(u_{i\mskip2mu\relax j})\mid a_{i\mskip2mu\relax j}\leq u_{i\mskip2mu\relax j}\leq b_{i\mskip2mu\relax j}\rbrace [ X ] = {( u i j ) ∣ a i j ≤ u i j ≤ b i j } m × n m\times n m × n X ▽ \bm{X}^{\triangledown} X ▽ X △ \bm{X}^{\vartriangle} X △
X ▽ ≔ ( a i j ) , X △ ≔ ( b i j ) \bm{X}^{\triangledown} \coloneqq (a_{i\mskip2mu\relax j}),
\quad\bm{X}^{\vartriangle} \coloneqq (b_{i\mskip2mu\relax j}) X ▽ : = ( a i j ) , X △ : = ( b i j ) で定義する．また，[ X ] \lbrack\bm{X}\rbrack [ X ]
[ X ] = [ X ▽ , X △ ] = ( [ x i j ] ) \lbrack\bm{X}\rbrack = \lbrack\bm{X}^{\triangledown},\bm{X}^{\vartriangle}\rbrack
= (\lbrack x_{i\mskip2mu\relax j}\rbrack) [ X ] = [ X ▽ , X △ ] = ([ x i j ]) とも表記する．ただし[ x i j ] = [ x i j ▽ , x i j △ ] = [ a i j , b i j ] \lbrack x_{i\mskip2mu\relax j}\rbrack=\lbrack x_{i\mskip2mu\relax j}^{\triangledown},x_{i\mskip2mu\relax j}^{\vartriangle}\rbrack=\lbrack a_{i\mskip2mu\relax j},b_{i\mskip2mu\relax j}\rbrack [ x i j ] = [ x i j ▽ , x i j △ ] = [ a i j , b i j ]
精度保証付き数値計算の文脈では，区間の下限をx ‾ \underline{x} x x ‾ \overline{x} x
S S S R \mathbb{R} R S S S I ∈ I R I\in\mathbb{IR} I ∈ IR
J ⊇ S ⟹ J ⊇ I ( J ∈ I R ) J \supseteq S
\implies J \supseteq I\quad(J\in\mathbb{IR}) J ⊇ S ⟹ J ⊇ I ( J ∈ IR ) を満たすものをS S S 区間包 (interval hull) といい，hull S \operatorname*{hull}S hull S
任意のS ⊆ R S\subseteq\mathbb{R} S ⊆ R S S S S S S [ inf S , sup S ] \lbrack\inf S,\sup S\rbrack [ inf S , sup S ] S S S θ \theta θ { f ( θ ) ∣ θ ∈ P } \lbrace f(\theta)\mid\theta\in P\rbrace { f ( θ ) ∣ θ ∈ P } hull S \operatorname*{hull}S hull S
hull θ ∈ P f ( θ ) \operatorname*{hull}_{\theta\in P}f(\theta) θ ∈ P hull f ( θ ) とも表す．
⋆ \star ⋆ R \mathbb{R} R [ x ] , [ y ] ∈ I R \lbrack x\rbrack,\lbrack y\rbrack\in\mathbb{IR} [ x ] , [ y ] ∈ IR [ x ] ⋆ [ y ] \lbrack x\rbrack\star\lbrack y\rbrack [ x ] ⋆ [ y ]
[ x ] ⋆ [ y ] ≔ hull { u ⋆ v ∣ ( u , v ) ∈ D ∩ ( [ x ] , [ y ] ) } \lbrack x\rbrack\star\lbrack y\rbrack \coloneqq \operatorname*{hull}\lbrace u\star v\mid(u,v)\in D\cap(\lbrack x\rbrack,\lbrack y\rbrack)\rbrace [ x ] ⋆ [ y ] : = hull { u ⋆ v ∣ ( u , v ) ∈ D ∩ ([ x ] , [ y ])} で定義する．ただし，D D D u ⋆ v u\star v u ⋆ v ( u , v ) ∈ R 2 (u,v)\in\mathbb{R}^{2} ( u , v ) ∈ R 2
⋆ \star ⋆ [ x ] ⋆ [ y ] \lbrack x\rbrack\star\lbrack y\rbrack [ x ] ⋆ [ y ]
[ x ] + [ y ] = [ x ▽ + y ▽ , x △ + y △ ] , [ x ] − [ y ] = [ x ▽ − y △ , x △ − y ▽ ] \begin{gathered}
\lbrack x\rbrack+\lbrack y\rbrack = \lbrack x^{\triangledown}+y^{\triangledown},x^{\vartriangle}+y^{\vartriangle}\rbrack,\\
\lbrack x\rbrack-\lbrack y\rbrack = \lbrack x^{\triangledown}-y^{\vartriangle},x^{\vartriangle}-y^{\triangledown}\rbrack
\end{gathered} [ x ] + [ y ] = [ x ▽ + y ▽ , x △ + y △ ] , [ x ] − [ y ] = [ x ▽ − y △ , x △ − y ▽ ] と書ける．⋆ \star ⋆
[ x ] ⋅ [ y ] = hull { x ▽ y ▽ , x △ y ▽ , x ▽ y △ , x △ y △ } , [ x ] [ y ] = hull { x ▽ y ▽ , x △ y ▽ , x ▽ y △ , x △ y △ } if 0 ∉ [ y ] \begin{gathered}
\lbrack x\rbrack\cdot\lbrack y\rbrack = \operatorname*{hull}\lbrace x^{\triangledown}y^{\triangledown},x^{\vartriangle}y^{\triangledown},x^{\triangledown}y^{\vartriangle},x^{\vartriangle}y^{\vartriangle}\rbrace,\\
\frac{\lbrack x\rbrack}{\lbrack y\rbrack} = \operatorname*{hull}\biggl\lbrace\frac{x^{\triangledown}}{y^{\triangledown}},\frac{x^{\vartriangle}}{y^{\triangledown}},\frac{x^{\triangledown}}{y^{\vartriangle}},\frac{x^{\vartriangle}}{y^{\vartriangle}}\biggr\rbrace\mskip18mu\relax\textrm{if}\mskip6mu\relax 0\notin\lbrack y\rbrack
\end{gathered} [ x ] ⋅ [ y ] = hull { x ▽ y ▽ , x △ y ▽ , x ▽ y △ , x △ y △ } , [ y ] [ x ] = hull { y ▽ x ▽ , y ▽ x △ , y △ x ▽ , y △ x △ } if 0 ∈ / [ y ] となる（0 ∈ [ y ] 0\in\lbrack y\rbrack 0 ∈ [ y ]
たとえば，[ x ] = [ 2 , 3 ] \lbrack x\rbrack=\lbrack 2,3\rbrack [ x ] = [ 2 , 3 ] [ y ] = [ − 5 , − 2 ] \lbrack y\rbrack=\lbrack-5,-2\rbrack [ y ] = [ − 5 , − 2 ]
[ x ] + [ y ] = [ − 3 , 1 ] , [ x ] − [ y ] = [ 4 , 8 ] , [ x ] ⋅ [ y ] = [ − 15 , − 4 ] , [ x ] [ y ] = [ − 1.5 , − 0.4 ] \begin{gathered}
\lbrack x\rbrack+\lbrack y\rbrack = \lbrack-3,1\rbrack,
\quad\lbrack x\rbrack-\lbrack y\rbrack = \lbrack 4,8\rbrack,\\
\lbrack x\rbrack\cdot\lbrack y\rbrack = \lbrack-15,-4\rbrack,
\quad\frac{\lbrack x\rbrack}{\lbrack y\rbrack} = \lbrack-1.5,-0.4\rbrack
\end{gathered} [ x ] + [ y ] = [ − 3 , 1 ] , [ x ] − [ y ] = [ 4 , 8 ] , [ x ] ⋅ [ y ] = [ − 15 , − 4 ] , [ y ] [ x ] = [ − 1.5 , − 0.4 ] である．
区間演算における減算は加算の逆演算ではなく，除算は乗算の逆演算ではない．実際，[ x ] = [ 0 , 1 ] \lbrack x\rbrack=\lbrack 0,1\rbrack [ x ] = [ 0 , 1 ] [ y ] = [ 2 , 5 ] \lbrack y\rbrack=\lbrack 2,5\rbrack [ y ] = [ 2 , 5 ]
[ x ] − [ x ] = [ − 1 , 1 ] ≠ 0 , [ y ] [ y ] = [ 0.4 , 2.5 ] ≠ 1 \lbrack x\rbrack-\lbrack x\rbrack = \lbrack-1,1\rbrack
\neq 0,
\quad\frac{\lbrack y\rbrack}{\lbrack y\rbrack} = \lbrack 0.4,2.5\rbrack
\neq 1 [ x ] − [ x ] = [ − 1 , 1 ] = 0 , [ y ] [ y ] = [ 0.4 , 2.5 ] = 1 である．また，[ z ] = [ − 1 , 3 ] \lbrack z\rbrack=\lbrack -1,3\rbrack [ z ] = [ − 1 , 3 ]
[ x ] ⋅ ( [ y ] + [ z ] ) = [ 0 , 1 ] ⋅ [ 1 , 8 ] = [ 0 , 8 ] , [ x ] ⋅ [ y ] + [ x ] ⋅ [ z ] = [ 0 , 5 ] + [ − 1 , 3 ] = [ − 1 , 8 ] \begin{gathered}
\lbrack x\rbrack\cdot(\lbrack y\rbrack+\lbrack z\rbrack) = \lbrack 0,1\rbrack\cdot\lbrack 1,8\rbrack
= \lbrack 0,8\rbrack,\\
\lbrack x\rbrack\cdot\lbrack y\rbrack+\lbrack x\rbrack\cdot\lbrack z\rbrack = \lbrack 0,5\rbrack+\lbrack -1,3\rbrack
= \lbrack -1,8\rbrack
\end{gathered} [ x ] ⋅ ([ y ] + [ z ]) = [ 0 , 1 ] ⋅ [ 1 , 8 ] = [ 0 , 8 ] , [ x ] ⋅ [ y ] + [ x ] ⋅ [ z ] = [ 0 , 5 ] + [ − 1 , 3 ] = [ − 1 , 8 ] となるので，分配法則も成り立たない．その代わり，劣分配法則 (sub-distributive law)
[ x ] ⋅ ( [ y ] + [ z ] ) ⊆ [ x ] ⋅ [ y ] + [ x ] ⋅ [ z ] \lbrack x\rbrack\cdot(\lbrack y\rbrack+\lbrack z\rbrack) \subseteq \lbrack x\rbrack\cdot\lbrack y\rbrack+\lbrack x\rbrack\cdot\lbrack z\rbrack [ x ] ⋅ ([ y ] + [ z ]) ⊆ [ x ] ⋅ [ y ] + [ x ] ⋅ [ z ] が成立する．
区間ベクトルと区間行列に対しても，和と差は成分ごとに定義される．たとえば[ x ] = ( [ x 1 ] [ x 2 ] ) T \lbrack\bm{x}\rbrack=(\begin{matrix}\lbrack x_{1}\rbrack & \lbrack x_{2}\rbrack\end{matrix})^{\mathsf{T}} [ x ] = ( [ x 1 ] [ x 2 ] ) T [ y ] = ( [ y 1 ] [ y 2 ] ) T \lbrack\bm{y}\rbrack=(\begin{matrix}\lbrack y_{1}\rbrack & \lbrack y_{2}\rbrack\end{matrix})^{\mathsf{T}} [ y ] = ( [ y 1 ] [ y 2 ] ) T
[ x ] ± [ y ] = ( [ x 1 ] ± [ y 1 ] [ x 2 ] ± [ y 2 ] ) T \lbrack\bm{x}\rbrack\pm\lbrack\bm{y}\rbrack = (\begin{matrix}\lbrack x_{1}\rbrack\pm\lbrack y_{1}\rbrack & \lbrack x_{2}\rbrack\pm\lbrack y_{2}\rbrack\end{matrix})^{\mathsf{T}} [ x ] ± [ y ] = ( [ x 1 ] ± [ y 1 ] [ x 2 ] ± [ y 2 ] ) T である．
機械区間演算
区間演算を計算機に実装しようとすると，扱える数が浮動小数点数に制限される．上限と下限がともにF ˉ \mathbb{\bar{F}} F ˉ 機械区間 (machine interval) といい，全体集合をI F \mathbb{IF} IF
S S S R \mathbb{R} R S S S I ∈ I F I\in\mathbb{IF} I ∈ IF
J ⊇ S ⟹ J ⊇ I ( J ∈ I F ) J \supseteq S
\implies J \supseteq I\quad(J\in\mathbb{IF}) J ⊇ S ⟹ J ⊇ I ( J ∈ IF ) を満たすものをS S S F \mathbb{F} F hull F S \operatorname{hull}_{\mathbb{F}}S hull F S
⋆ \star ⋆ R \mathbb{R} R [ x ] , [ y ] ∈ I F \lbrack x\rbrack,\lbrack y\rbrack\in\mathbb{IF} [ x ] , [ y ] ∈ IF [ x ] ⋆ [ y ] \lbrack x\rbrack\star\lbrack y\rbrack [ x ] ⋆ [ y ]
[ x ] ⋆ [ y ] ≔ hull F { u ⋆ v ∣ ( u , v ) ∈ D ∩ ( [ x ] , [ y ] ) } \lbrack x\rbrack\star\lbrack y\rbrack \coloneqq \operatorname{hull}_{\mathbb{F}}\lbrace u\star v\mid(u,v)\in D\cap(\lbrack x\rbrack,\lbrack y\rbrack)\rbrace [ x ] ⋆ [ y ] : = hull F { u ⋆ v ∣ ( u , v ) ∈ D ∩ ([ x ] , [ y ])} で定義する．ただし，D D D u ⋆ v u\star v u ⋆ v ( u , v ) ∈ R 2 (u,v)\in\mathbb{R}^{2} ( u , v ) ∈ R 2
機械区間演算について，⋆ \star ⋆
[ x ] + [ y ] = [ RD ( x ▽ + y ▽ ) , RU ( x △ + y △ ) ] , [ x ] − [ y ] = [ RD ( x ▽ − y △ ) , RU ( x △ − y ▽ ) ] \begin{gathered}
\lbrack x\rbrack+\lbrack y\rbrack = \lbrack\operatorname{RD}(x^{\triangledown}+y^{\triangledown}),\operatorname{RU}(x^{\vartriangle}+y^{\vartriangle})\rbrack,\\
\lbrack x\rbrack-\lbrack y\rbrack = \lbrack\operatorname{RD}(x^{\triangledown}-y^{\vartriangle}),\operatorname{RU}(x^{\vartriangle}-y^{\triangledown})\rbrack
\end{gathered} [ x ] + [ y ] = [ RD ( x ▽ + y ▽ ) , RU ( x △ + y △ )] , [ x ] − [ y ] = [ RD ( x ▽ − y △ ) , RU ( x △ − y ▽ )] となる．⋆ \star ⋆ Yamanaka (2015) を参照のこと．
べき級数演算
区間演算は，実数の演算を浮動小数点数の演算で不等式評価する技法といえる．同様に，関数の演算を浮動小数点数係数多項式の演算で不等式評価する技法が，Kashiwagi (1994) のべき級数演算である．
簡単のため，この記事ではI R \mathbb{IR} IR I R \mathbb{IR} IR I F \mathbb{IF} IF
文字T T T
f ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n ( [ a i ] ∈ I R ) f(T) = \lbrack a_{0}\rbrack+\lbrack a_{1}\rbrack T+\dotsb+\lbrack a_{n}\rbrack T^{n}\quad(\lbrack a_{i}\rbrack\in\mathbb{IR}) f ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n ([ a i ] ∈ IR ) を区間多項式 (interval polynomial) という．
実数t t t
f ( t ) ≔ ∑ i = 0 n [ a i ] t i = { ∑ i = 0 n u i t i ∣ u i ∈ [ a i ] } f(t) \coloneqq \sum_{i=0}^{n}\lbrack a_{i}\rbrack t^{i}
= \Biggl\lbrace\sum_{i=0}^{n}u_{i}t^{i}\Biggm\vert u_{i}\in\lbrack a_{i}\rbrack\Biggr\rbrace f ( t ) : = i = 0 ∑ n [ a i ] t i = { i = 0 ∑ n u i t i u i ∈ [ a i ] } と定める．f ( t ) f(t) f ( t ) f ( T ) ∣ T = t f(T)\rvert_{T=t} f ( T ) ∣ T = t
f ( T ) f(T) f ( T ) D D D x = inf f ( t ) x=\inf f(t) x = inf f ( t ) x = sup f ( t ) x=\sup f(t) x = sup f ( t ) x = ϕ ( t ) x=\phi(t) x = ϕ ( t ) t ∈ D t\in D t ∈ D band D f ( T ) \operatorname*{band}_{D}f(T) band D f ( T )
始域X X X Y Y Y C ( X , Y ) C(X,Y) C ( X , Y ) D D D band D f ( T ) \operatorname*{band}_{D}f(T) band D f ( T )
band D f ( T ) ≔ { ϕ ∈ C ( D , R ) ∣ ϕ ( t ) ∈ f ( t ) for all t } \operatorname*{band}_{D}f(T) \coloneqq \lbrace\phi\in C(D,\mathbb{R})\mid\phi(t)\in f(t)\mskip18mu\relax\textrm{for all}\mskip6mu\relax t\rbrace D band f ( T ) : = { ϕ ∈ C ( D , R ) ∣ ϕ ( t ) ∈ f ( t ) for all t } で定義する．D D D band f ( T ) \operatorname*{band}f(T) band f ( T )
たいていの文献でband f ( T ) \operatorname*{band}f(T) band f ( T ) f ( T ) f(T) f ( T )
連続関数ϕ ( t ) \phi(t) ϕ ( t ) band f ( T ) \operatorname*{band}f(T) band f ( T )
inf f ( t ) ≤ ϕ ( t ) ≤ sup f ( t ) for all t \inf f(t) \leq \phi(t)
\leq \sup f(t)\mskip18mu\relax\textrm{for all}\mskip6mu\relax t inf f ( t ) ≤ ϕ ( t ) ≤ sup f ( t ) for all t が成り立つことだから，区間多項式は関数に対する不等式評価に相当する．
この節では，区間多項式に対して，以下の基本的な演算を定める．どの演算も，不等式評価と矛盾しないように定義されることが重要である．
加減乗除
関数合成
積分
減次
区間多項式f ( T ) = ∑ i = 0 m [ a i ] T i f(T)=\sum_{i=0}^{m}\lbrack a_{i}\rbrack T^{i} f ( T ) = ∑ i = 0 m [ a i ] T i g ( T ) = ∑ i = 0 n [ b i ] T i g(T)=\sum_{i=0}^{n}\lbrack b_{i}\rbrack T^{i} g ( T ) = ∑ i = 0 n [ b i ] T i
f ( T ) ± g ( T ) ≔ ∑ i = 0 max { m , n } ( [ a i ] ± [ b i ] ) T i , f ( T ) g ( T ) ≔ ∑ i = 0 m + n ( ∑ j = 0 i [ a j ] [ b i − j ] ) T i \begin{gathered}
f(T)\pm g(T) \coloneqq \sum_{i=0}^{\max\lbrace m,n\rbrace}(\lbrack a_{i}\rbrack\pm\lbrack b_{i}\rbrack) T^{i},\\
f(T)g(T) \coloneqq \sum_{i=0}^{m+n}\Biggl(\sum_{j=0}^{i}\lbrack a_{j}\rbrack\lbrack b_{i-j}\rbrack\Biggr)T^{i}
\end{gathered} f ( T ) ± g ( T ) : = i = 0 ∑ m a x { m , n } ([ a i ] ± [ b i ]) T i , f ( T ) g ( T ) : = i = 0 ∑ m + n ( j = 0 ∑ i [ a j ] [ b i − j ] ) T i で定義される．ただし，[ a m + 1 ] = [ a m + 2 ] = ⋯ = 0 \lbrack a_{m+1}\rbrack=\lbrack a_{m+2}\rbrack=\dotsb=0 [ a m + 1 ] = [ a m + 2 ] = ⋯ = 0 [ b n + 1 ] = [ b n + 2 ] = ⋯ = 0 \lbrack b_{n+1}\rbrack=\lbrack b_{n+2}\rbrack=\dotsb=0 [ b n + 1 ] = [ b n + 2 ] = ⋯ = 0
演算⋆ \mathord{\star} ⋆ ϕ 1 ∈ band D f ( T ) \phi_{1}\in\operatorname*{band}_{D}f(T) ϕ 1 ∈ band D f ( T ) ϕ 2 ∈ band D g ( T ) \phi_{2}\in\operatorname*{band}_{D}g(T) ϕ 2 ∈ band D g ( T ) t t t ϕ 1 ( t ) ⋆ ϕ 2 ( t ) \phi_{1}(t)\star\phi_{2}(t) ϕ 1 ( t ) ⋆ ϕ 2 ( t ) band D ( f ( T ) ⋆ g ( T ) ) \operatorname*{band}_{D}(f(T)\star g(T)) band D ( f ( T ) ⋆ g ( T ))
⋆ \star ⋆ ϕ 1 ∈ band f ( T ) \phi_{1}\in\operatorname*{band}f(T) ϕ 1 ∈ band f ( T ) ϕ 2 ∈ band g ( T ) \phi_{2}\in\operatorname*{band}g(T) ϕ 2 ∈ band g ( T ) t ∈ D t\in D t ∈ D u i ∈ [ a i ] u_{i}\in\lbrack a_{i}\rbrack u i ∈ [ a i ] v j ∈ [ b j ] v_{j}\in\lbrack b_{j}\rbrack v j ∈ [ b j ]
ϕ 1 ( t ) = ∑ i = 0 m u i t i , ϕ 2 ( t ) = ∑ j = 0 n v j t j \phi_{1}(t) = \sum_{i=0}^{m}u_{i}t^{i},
\quad\phi_{2}(t) = \sum_{j=0}^{n}v_{j}t^{j} ϕ 1 ( t ) = i = 0 ∑ m u i t i , ϕ 2 ( t ) = j = 0 ∑ n v j t j を満たす．k = i + j k=i+j k = i + j
ϕ 1 ( t ) ϕ 2 ( t ) = ∑ i = 0 m ∑ j = 0 n u i v j t i + j = ∑ k = 0 m + n ∑ i = 0 k u i v k − i t k , ∑ k = 0 m + n ∑ i = 0 k u i v k − i t k ∈ ∑ k = 0 m + n ( ∑ i = 0 k [ a i ] [ b k − i ] ) t k \begin{gathered}
\phi_{1}(t)\phi_{2}(t) = \sum_{i=0}^{m}\sum_{j=0}^{n}u_{i}v_{j}t^{i+j}
= \sum_{k=0}^{m+n}\sum_{i=0}^{k}u_{i}v_{k-i}t^{k},\\
\sum_{k=0}^{m+n}\sum_{i=0}^{k}u_{i}v_{k-i}t^{k} \in \sum_{k=0}^{m+n}\Biggl(\sum_{i=0}^{k}\lbrack a_{i}\rbrack\lbrack b_{k-i}\rbrack\Biggr)t^{k}
\end{gathered} ϕ 1 ( t ) ϕ 2 ( t ) = i = 0 ∑ m j = 0 ∑ n u i v j t i + j = k = 0 ∑ m + n i = 0 ∑ k u i v k − i t k , k = 0 ∑ m + n i = 0 ∑ k u i v k − i t k ∈ k = 0 ∑ m + n ( i = 0 ∑ k [ a i ] [ b k − i ] ) t k なのでϕ 1 ( t ) ϕ 2 ( t ) ∈ ( f ( T ) g ( T ) ) ∣ T = t \phi_{1}(t)\phi_{2}(t)\in(f(T)g(T))\rvert_{T=t} ϕ 1 ( t ) ϕ 2 ( t ) ∈ ( f ( T ) g ( T )) ∣ T = t
D D D 0 0 0 y = ψ ( x ) y=\psi(x) y = ψ ( x ) [ R ] = ⋃ t ∈ D f ( t ) \lbrack R\rbrack=\bigcup_{t\in D}f(t) [ R ] = ⋃ t ∈ D f ( t ) ψ ( x ) \psi(x) ψ ( x )
ψ ( f ( T ) ) ≔ ∑ i = 0 k − 1 1 i ! hull x ∈ [ a 0 ] ( d i y d x i ) ( [ a 1 ] T + [ a 2 ] T 2 + ⋯ ) i ≔ + 1 k ! hull x ∈ [ R ] ( d k y d x k ) ( [ a 1 ] T + [ a 2 ] T 2 + ⋯ ) k \begin{aligned}
\psi(f(T)) &\coloneqq \sum_{i=0}^{k-1}\frac{1}{i!}\operatorname*{hull}_{x\in\lbrack a_{0}\rbrack}\biggl(\frac{\mathrm{d}^{i}y}{\mathrm{d}x^{i}}\biggr)(\lbrack a_{1}\rbrack T+\lbrack a_{2}\rbrack T^{2}+\dotsb)^{i}\\
&\hphantom{{}\coloneqq{}}+\frac{1}{k!}\operatorname*{hull}_{x\in\lbrack R\rbrack}\biggl(\frac{\mathrm{d}^{k}y}{\mathrm{d}x^{k}}\biggr)(\lbrack a_{1}\rbrack T+\lbrack a_{2}\rbrack T^{2}+\dotsb)^{k}
\end{aligned} ψ ( f ( T )) : = i = 0 ∑ k − 1 i ! 1 x ∈ [ a 0 ] hull ( d x i d i y ) ([ a 1 ] T + [ a 2 ] T 2 + ⋯ ) i : = + k ! 1 x ∈ [ R ] hull ( d x k d k y ) ([ a 1 ] T + [ a 2 ] T 2 + ⋯ ) k とおく．ただし，k k k f ( T ) 0 = 1 f(T)^{0}=1 f ( T ) 0 = 1
すべてのϕ ∈ band D f ( T ) \phi\in\operatorname*{band}_{D}f(T) ϕ ∈ band D f ( T ) ψ ( ϕ ( t ) ) \psi(\phi(t)) ψ ( ϕ ( t )) band D ψ ( f ( T ) ) \operatorname*{band}_{D}\psi(f(T)) band D ψ ( f ( T ))
t ∈ D t\in D t ∈ D u i ∈ [ a i ] u_{i}\in\lbrack a_{i}\rbrack u i ∈ [ a i ] ϕ ( t ) = ∑ u i t i \phi(t)=\sum u_{i}t^{i} ϕ ( t ) = ∑ u i t i u 0 ∈ f ( 0 ) ⊆ [ R ] u_{0}\in f(0)\subseteq\lbrack R\rbrack u 0 ∈ f ( 0 ) ⊆ [ R ] x 0 = ϕ ( t ) x_{0}=\phi(t) x 0 = ϕ ( t ) ξ ∈ hull { u 0 , x 0 } \xi\in\operatorname*{hull}\lbrace u_{0},x_{0}\rbrace ξ ∈ hull { u 0 , x 0 }
ψ ( x 0 ) = ∑ i = 0 k − 1 1 i ! d i y d x i ∣ x = u 0 ( x 0 − u 0 ) i = + 1 k ! d k y d x k ∣ x = ξ ( x 0 − u 0 ) k \begin{aligned}
\psi(x_{0}) &= \sum_{i=0}^{k-1}\frac{1}{i!}\frac{\mathrm{d}^{i}y}{\mathrm{d}x^{i}}\biggr\rvert_{x=u_{0}}(x_{0}-u_{0})^{i}\\
&\hphantom{{}={}}+\frac{1}{k!}\frac{\mathrm{d}^{k}y}{\mathrm{d}x^{k}}\biggr\rvert_{x=\xi}(x_{0}-u_{0})^{k}
\end{aligned} ψ ( x 0 ) = i = 0 ∑ k − 1 i ! 1 d x i d i y x = u 0 ( x 0 − u 0 ) i = + k ! 1 d x k d k y x = ξ ( x 0 − u 0 ) k を満たす（テイラーの定理）．よって
ψ ( ϕ ( t ) ) = ∑ i = 0 k − 1 1 i ! d i y d x i ∣ x = u 0 ( u 1 t + u 2 t 2 + ⋯ ) i = + 1 k ! d k y d x k ∣ x = ξ ( u 1 t + u 2 t 2 + ⋯ ) k \begin{aligned}
\psi(\phi(t)) &= \sum_{i=0}^{k-1}\frac{1}{i!}\frac{\mathrm{d}^{i}y}{\mathrm{d}x^{i}}\biggr\rvert_{x=u_{0}}(u_{1}t+u_{2}t^{2}+\dotsb)^{i}\\
&\hphantom{{}={}}+\frac{1}{k!}\frac{\mathrm{d}^{k}y}{\mathrm{d}x^{k}}\biggr\rvert_{x=\xi}(u_{1}t+u_{2}t^{2}+\dotsb)^{k}
\end{aligned} ψ ( ϕ ( t )) = i = 0 ∑ k − 1 i ! 1 d x i d i y x = u 0 ( u 1 t + u 2 t 2 + ⋯ ) i = + k ! 1 d x k d k y x = ξ ( u 1 t + u 2 t 2 + ⋯ ) k なのでψ ( ϕ ( t ) ) ∈ ψ ( f ( T ) ) ∣ T = t \psi(\phi(t))\in\psi(f(T))\rvert_{T=t} ψ ( ϕ ( t )) ∈ ψ ( f ( T )) ∣ T = t
この命題を使うと，区間多項式の除算を定義できる．y = 1 / x y=1/x y = 1/ x
1 i ! d i y d x i = ( − 1 ) i x i + 1 \frac{1}{i!}\frac{\mathrm{d}^{i}y}{\mathrm{d}x^{i}} = \frac{(-1)^{i}}{x^{i+1}} i ! 1 d x i d i y = x i + 1 ( − 1 ) i なので
g ( T ) f ( T ) ≔ g ( T ) ⋅ 1 f ( T ) , 1 f ( T ) ≔ ∑ i = 0 k − 1 hull x ∈ [ a 0 ] ( ( − 1 ) i x i + 1 ) ( [ a 1 ] T + [ a 2 ] T 2 + ⋯ ) i ≔ + hull x ∈ [ R ] ( ( − 1 ) k x k + 1 ) ( [ a 1 ] T + [ a 2 ] T 2 + ⋯ ) k \begin{gathered}
\frac{g(T)}{f(T)} \coloneqq g(T)\cdot\frac{1}{f(T)},\\
\begin{aligned}
\frac{1}{f(T)} &\coloneqq \sum_{i=0}^{k-1}\operatorname*{hull}_{x\in\lbrack a_{0}\rbrack}\biggl(\frac{(-1)^{i}}{x^{i+1}}\biggr)(\lbrack a_{1}\rbrack T+\lbrack a_{2}\rbrack T^{2}+\dotsb)^{i}\\
&\hphantom{{}\coloneqq{}}+\operatorname*{hull}_{x\in\lbrack R\rbrack}\biggl(\frac{(-1)^{k}}{x^{k+1}}\biggr)(\lbrack a_{1}\rbrack T+\lbrack a_{2}\rbrack T^{2}+\dotsb)^{k}
\end{aligned}
\end{gathered} f ( T ) g ( T ) : = g ( T ) ⋅ f ( T ) 1 , f ( T ) 1 : = i = 0 ∑ k − 1 x ∈ [ a 0 ] hull ( x i + 1 ( − 1 ) i ) ([ a 1 ] T + [ a 2 ] T 2 + ⋯ ) i : = + x ∈ [ R ] hull ( x k + 1 ( − 1 ) k ) ([ a 1 ] T + [ a 2 ] T 2 + ⋯ ) k とすればよい．ただし，関数y = 1 / x y=1/x y = 1/ x [ R ] \lbrack R\rbrack [ R ] 0 ∉ [ R ] 0\notin\lbrack R\rbrack 0 ∈ / [ R ] k = 1 k=1 k = 1
1 f ( T ) = 1 [ a 0 ] − ∑ i = 1 n [ a i ] hull x ∈ [ R ] ( x 2 ) T i \frac{1}{f(T)} = \frac{1}{\lbrack a_{0}\rbrack}-\sum_{i=1}^{n}\frac{\lbrack a_{i}\rbrack}{\operatorname*{hull}_{x\in\lbrack R\rbrack}(x^{2})}T^{i} f ( T ) 1 = [ a 0 ] 1 − i = 1 ∑ n hull x ∈ [ R ] ( x 2 ) [ a i ] T i であり，0 ∉ [ R ] 0\notin\lbrack R\rbrack 0 ∈ / [ R ]
1 f ( T ) = 1 [ a 0 ] − ∑ i = 1 n [ a i ] [ R ] 2 T i ( [ R ] = ⋃ t ∈ D f ( t ) ) \frac{1}{f(T)} = \frac{1}{\lbrack a_{0}\rbrack}-\sum_{i=1}^{n}\frac{\lbrack a_{i}\rbrack}{\lbrack R\rbrack^{2}}T^{i}\quad({\textstyle\lbrack R\rbrack=\bigcup_{t\in D}f(t)}) f ( T ) 1 = [ a 0 ] 1 − i = 1 ∑ n [ R ] 2 [ a i ] T i ( [ R ] = ⋃ t ∈ D f ( t ) ) となる．
区間多項式f ( T ) = ∑ i = 0 n [ a i ] T i f(T)=\sum_{i=0}^{n}\lbrack a_{i}\rbrack T^{i} f ( T ) = ∑ i = 0 n [ a i ] T i
∫ f ( T ) d T ≔ ∑ i = 0 n [ a i ] i + 1 T i + 1 \int f(T)\,\mathrm{d}T \coloneqq \sum_{i=0}^{n}\frac{\lbrack a_{i}\rbrack}{i+1}T^{i+1} ∫ f ( T ) d T : = i = 0 ∑ n i + 1 [ a i ] T i + 1 で定義する．
0 ∈ D 0\in D 0 ∈ D ϕ ∈ band D f ( T ) \phi\in\operatorname*{band}_{D}f(T) ϕ ∈ band D f ( T ) t t t ∫ 0 t ϕ ( t ′ ) d t ′ \int_{0}^{t}\phi(t')\,\mathrm{d}t' ∫ 0 t ϕ ( t ′ ) d t ′ band D ( ∫ f ( T ) d T ) \operatorname*{band}_{D}(\int f(T)\,\mathrm{d}T) band D ( ∫ f ( T ) d T )
区間多項式のm m m
∑ i [ a i ] T i ⟶ ∑ i < m [ a i ] T i + ( ⋃ t ∈ D ∑ i ≥ m [ a i ] t i − m ) T m \sum_{i}\lbrack a_{i}\rbrack T^{i} \longrightarrow \sum_{i\lt m}\lbrack a_{i}\rbrack T^{i}+\Biggl(\bigcup_{t\in D}\sum_{i\geq m}\lbrack a_{i}\rbrack t^{i-m}\Biggr)T^{m} i ∑ [ a i ] T i ⟶ i < m ∑ [ a i ] T i + ( t ∈ D ⋃ i ≥ m ∑ [ a i ] t i − m ) T m で定義する．
定義から，減次f ( T ) → g ( T ) f(T)\to g(T) f ( T ) → g ( T )
band D f ( T ) ⊆ band D g ( T ) \operatorname*{band}_{D}f(T) \subseteq \operatorname*{band}_{D}g(T) D band f ( T ) ⊆ D band g ( T ) が成立する．
以上で，I R \mathbb{IR} IR I R d \mathbb{IR}^{d} IR d I R \mathbb{IR} IR
文字T T T
f ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n ( [ a i ] ∈ I R d ) \bm{f}(T) = \lbrack\bm{a}_{0}\rbrack+\lbrack\bm{a}_{1}\rbrack T+\dotsb+\lbrack\bm{a}_{n}\rbrack T^{n}\quad(\lbrack\bm{a}_{i}\rbrack\in\mathbb{IR}^{d}) f ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n ([ a i ] ∈ IR d ) を区間多項式という．
f ( T ) \bm{f}(T) f ( T ) [ a i j ] \lbrack a_{i\mskip2mu\relax j}\rbrack [ a i j ] f i ( T ) f_{i}(T) f i ( T )
( [ a 1 j ] [ a 2 j ] ⋮ [ a d j ] ) = [ a j ] , f i ( T ) = ∑ j = 1 n [ a i j ] T j \mathopen{}\mathclose{\left(\begin{matrix}\lbrack a_{1\mskip2mu\relax j}\rbrack \\ \lbrack a_{2\mskip2mu\relax j}\rbrack \\ \vdots \\ \lbrack a_{d\mskip2mu\relax j}\rbrack\end{matrix}\right)} = \lbrack\bm{a}_{j}\rbrack,
\quad f_{i}(T) = \sum_{j=1}^{n}\lbrack a_{i\mskip2mu\relax j}\rbrack T^{j} [ a 1 j ] [ a 2 j ] ⋮ [ a d j ] = [ a j ] , f i ( T ) = j = 1 ∑ n [ a i j ] T j で定めると，f ( T ) \bm{f}(T) f ( T ) f 1 ( T ) , f 2 ( T ) , … , f d ( T ) f_{1}(T),f_{2}(T),\dotsc,f_{d}(T) f 1 ( T ) , f 2 ( T ) , … , f d ( T )
( f 1 ( T ) f 2 ( T ) ⋮ f d ( T ) ) \mathopen{}\mathclose{\left(\begin{matrix}f_{1}(T) \\ f_{2}(T) \\ \vdots \\ f_{d}(T)\end{matrix}\right)} f 1 ( T ) f 2 ( T ) ⋮ f d ( T ) と同一視できる．f ( T ) \bm{f}(T) f ( T )
区間多項式と基本演算の組をべき級数演算 (power series arithmetic; PSA) という．区間多項式の次数をm m m
∑ i [ a i ] T i ⟶ ∑ i ≤ m [ a i ] T i \sum_{i}\lbrack a_{i}\rbrack T^{i} \longrightarrow \sum_{i\leq m}\lbrack a_{i}\rbrack T^{i} i ∑ [ a i ] T i ⟶ i ≤ m ∑ [ a i ] T i で次数下げする算法をType-I PSA といい，減次
∑ i [ a i ] T i ⟶ ∑ i < m [ a i ] T i + ( ⋃ t ∈ D ∑ i ≥ m [ a i ] t i − m ) T m \sum_{i}\lbrack a_{i}\rbrack T^{i} \longrightarrow \sum_{i\lt m}\lbrack a_{i}\rbrack T^{i}+\Biggl(\bigcup_{t\in D}\sum_{i\geq m}\lbrack a_{i}\rbrack t^{i-m}\Biggr)T^{m} i ∑ [ a i ] T i ⟶ i < m ∑ [ a i ] T i + ( t ∈ D ⋃ i ≥ m ∑ [ a i ] t i − m ) T m で次数下げする算法をType-II PSA という．
べき級数演算のこの定義は原形から少しずれている．本来の定義は柏木雅英（2018），または柏木啓一郎・柏木雅英（2011） を参照のこと．
次数下げをしなければ，べき級数演算によって区間多項式の次数は限りなく大きくなりうる．計算資源には限りがあるので，現実的にはどこかで次数下げが必要になる．Type-I PSAは，求めたいものが多項式の係数であり，band f ( T ) \operatorname*{band}f(T) band f ( T )
ピカール反復
べき級数演算は，初期値問題の精度保証付き数値計算に応用される．この節では，初期値問題の数学について簡単に説明する．
微分方程式と初期条件の組
{ d x / d t = f ( x , t ) at t > t 0 , x = x 0 at t = t 0 \begin{cases}
\mathrm{d}\bm{x}/\mathrm{d}t = \bm{f}(\bm{x},t)\mskip12mu\relax\textrm{at}\mskip6mu\relax t\gt t_{0},\\
\bm{x} = \bm{x}_{0}\mskip12mu\relax\textrm{at}\mskip6mu\relax t=t_{0}
\end{cases} { d x / d t = f ( x , t ) at t > t 0 , x = x 0 at t = t 0 を初期値問題 (initial value problem; IVP) という．d x / d t = f ( x , t ) \mathrm{d}\bm{x}/\mathrm{d}t=\bm{f}(\bm{x},t) d x / d t = f ( x , t ) t t t
x = ∫ d x d t d t = ∫ f ( x , t ) d t \bm{x} = \int\frac{\mathrm{d}\bm{x}}{\mathrm{d}t}\,\mathrm{d}t
= \int\bm{f}(\bm{x},t)\,\mathrm{d}t x = ∫ d t d x d t = ∫ f ( x , t ) d t となるから，解x = ϕ ( t ) \bm{x}=\bm{\phi}(t) x = ϕ ( t )
ϕ ( t ) = x 0 + ∫ t 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ \bm{\phi}(t) = \bm{x}_{0}+\int_{t_{0}}^{t}\bm{f}(\bm{\phi}(t'),t')\,\mathrm{d}t' ϕ ( t ) = x 0 + ∫ t 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ を満たす．両辺をt t t ϕ \bm{\phi} ϕ ϕ \bm{\phi} ϕ P f ( ϕ ) P_{\bm{f}}(\bm{\phi}) P f ( ϕ )
P f ( ϕ ) : t ↦ x 0 + ∫ t 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ P_{\bm{f}}(\bm{\phi})\colon t \mapsto \bm{x}_{0}+\int_{t_{0}}^{t}\bm{f}(\bm{\phi}(t'),t')\,\mathrm{d}t' P f ( ϕ ) : t ↦ x 0 + ∫ t 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ で定義すると，初期値問題を解くことは，連続関数ϕ \bm{\phi} ϕ
ϕ = P f ( ϕ ) \bm{\phi} = P_{\bm{f}}(\bm{\phi}) ϕ = P f ( ϕ ) を解くことと同値である．ϕ = P f ( ϕ ) \bm{\phi}=P_{\bm{f}}(\bm{\phi}) ϕ = P f ( ϕ ) ϕ \bm{\phi} ϕ P f P_{\bm{f}} P f 不動点 (fixed point) という．
以下ではt 0 = 0 t_{0}=0 t 0 = 0 ϕ 0 , ϕ 1 , … \bm{\phi}_{0},\bm{\phi}_{1},\dotsc ϕ 0 , ϕ 1 , …
ϕ 0 ( t ) = x 0 , ϕ n + 1 ( t ) = x 0 + ∫ 0 t f ( ϕ n ( t ′ ) , t ′ ) d t ′ \bm{\phi}_{0}(t) = \bm{x}_{0},
\quad\bm{\phi}_{n+1}(t) = \bm{x}_{0}+\int_{0}^{t}\bm{f}(\bm{\phi}_{n}(t'),t')\,\mathrm{d}t' ϕ 0 ( t ) = x 0 , ϕ n + 1 ( t ) = x 0 + ∫ 0 t f ( ϕ n ( t ′ ) , t ′ ) d t ′ で定める．この関数列はϕ n + 1 = P f ( ϕ n ) \bm{\phi}_{n+1}=P_{\bm{f}}(\bm{\phi}_{n}) ϕ n + 1 = P f ( ϕ n ) n → ∞ n\to\infty n → ∞
ϕ ∞ = P f ( ϕ ∞ ) \bm{\phi}_{\infty} = P_{\bm{f}}(\bm{\phi}_{\infty}) ϕ ∞ = P f ( ϕ ∞ ) となることが予想される．つまり，極限関数ϕ ∞ \bm{\phi}_{\infty} ϕ ∞ ピカール反復 (Picard iteration) という．
ピカール反復について，次の定理が成り立つ．証明は付録に回す．
U U U R d \mathbb{R}^{d} R d x 0 \bm{x}_{0} x 0 U U U D = [ 0 , h ] ⊆ R D=\lbrack 0,h\rbrack\subseteq\mathbb{R} D = [ 0 , h ] ⊆ R
f : { ( x , t ) ∣ x ∈ U and t ∈ D } → R d \bm{f}\colon\lbrace(\bm{x},t)\mid\bm{x}\in U\;\mathrel{\textrm{and}}\;t\in D\rbrace \to \mathbb{R}^{d} f : {( x , t ) ∣ x ∈ U and t ∈ D } → R d はリプシッツ条件を満たすとする．すなわち，定数L > 0 L\gt 0 L > 0 x , x ′ ∈ U \bm{x},\bm{x}'\in U x , x ′ ∈ U
sup t ∈ D ∣ f ( x , t ) − f ( x ′ , t ) ∣ ≤ L ∣ x − x ′ ∣ \sup_{t\in D}\lvert\bm{f}(\bm{x},t)-\bm{f}(\bm{x}',t)\rvert \leq L\lvert\bm{x}-\bm{x}'\rvert t ∈ D sup ∣ f ( x , t ) − f ( x ′ , t )∣ ≤ L ∣ x − x ′ ∣ となると仮定する．一様ノルムに関する閉部分集合K ⊆ C ( D , R d ) K\subseteq C(D,\mathbb{R}^{d}) K ⊆ C ( D , R d ) K K K P f P_{\bm{f}} P f
K ≠ ∅ K\neq\emptyset K = ∅ すべてのϕ ∈ K \bm{\phi}\in K ϕ ∈ K t ∈ D t\in D t ∈ D ϕ ( t ) ∈ U \bm{\phi}(t)\in U ϕ ( t ) ∈ U
すべてのϕ ∈ K \bm{\phi}\in K ϕ ∈ K P f ( ϕ ) ∈ K P_{\bm{f}}(\bm{\phi})\in K P f ( ϕ ) ∈ K
K K K K K K K K K
柏木の方法
この節では，Kashiwagi (1994) により提案された，初期値問題の精度保証付き数値解法の概要を示す．
重要なのは，p ( T ) \bm{p}(T) p ( T ) band D p ( T ) \operatorname*{band}_{D}\bm{p}(T) band D p ( T ) C ( D , R d ) C(D,\mathbb{R}^{d}) C ( D , R d )
p ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n \bm{p}(T) = \lbrack\bm{a}_{0}\rbrack+\lbrack\bm{a}_{1}\rbrack T+\dotsb+\lbrack\bm{a}_{n}\rbrack T^{n} p ( T ) = [ a 0 ] + [ a 1 ] T + ⋯ + [ a n ] T n を得たとする．K = band D p ( T ) K=\operatorname{band}_{D}\bm{p}(T) K = band D p ( T ) K K K
関数ϕ ∈ band D p ( T ) \bm{\phi}\in\operatorname{band}_{D}\bm{p}(T) ϕ ∈ band D p ( T ) P f ( ϕ ) P_{\bm{f}}(\bm{\phi}) P f ( ϕ ) band D p ( T ) \operatorname{band}_{D}\bm{p}(T) band D p ( T ) t ∈ D t\in D t ∈ D
x 0 + ∫ 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ ∈ p ( t ) \bm{x}_{0}+\int_{0}^{t}\bm{f}(\bm{\phi}(t'),t')\,\mathrm{d}t' \in \bm{p}(t) x 0 + ∫ 0 t f ( ϕ ( t ′ ) , t ′ ) d t ′ ∈ p ( t ) が成り立つことである．式
x 0 + ∫ f ( p ( T ) , T ) d T \bm{x}_{0}+\int\bm{f}(\bm{p}(T),T)\,\mathrm{d}T x 0 + ∫ f ( p ( T ) , T ) d T をType-II PSAで計算して，得られる区間多項式をp ( T ) ′ = ∑ [ a i ′ ] T i \bm{p}(T)'=\sum\lbrack\bm{a}_{i}'\rbrack T^{i} p ( T ) ′ = ∑ [ a i ′ ] T i i i i [ a i ′ ] ⊆ [ a i ] \lbrack\bm{a}_{i}'\rbrack\subseteq\lbrack\bm{a}_{i}\rbrack [ a i ′ ] ⊆ [ a i ]
P f ( ϕ ) ∈ band D p ( T ) ′ ⊆ band D p ( T ) P_{\bm{f}}(\bm{\phi}) \in \operatorname*{band}_{D}\bm{p}(T)' \subseteq \operatorname*{band}_{D}\bm{p}(T) P f ( ϕ ) ∈ D band p ( T ) ′ ⊆ D band p ( T ) より，K = band D p ( T ) K=\operatorname*{band}_{D}\bm{p}(T) K = band D p ( T )
p ( T ) \bm{p}(T) p ( T )
p 0 ( T ) = x 0 , p n + 1 ( T ) = x 0 + ∫ f ( p n ( T ) , T ) d T \bm{p}_{0}(T) = \bm{x}_{0},
\quad\bm{p}_{n+1}(T) = \bm{x}_{0}+\int\bm{f}(\bm{p}_{n}(T),T)\,\mathrm{d}T p 0 ( T ) = x 0 , p n + 1 ( T ) = x 0 + ∫ f ( p n ( T ) , T ) d T をType-I PSAで計算してp 1 ( T ) , p 2 ( T ) , … \bm{p}_{1}(T),\bm{p}_{2}(T),\dotsc p 1 ( T ) , p 2 ( T ) , … p n ( T ) \bm{p}_{n}(T) p n ( T ) n n n p n ( T ) \bm{p}_{n}(T) p n ( T ) x = ϕ ( t ) \bm{x}=\bm{\phi}(t) x = ϕ ( t )
x 0 + ϕ ′ ( 0 ) t 1 ! + ϕ ′ ′ ( 0 ) t 2 2 ! + ⋯ + ϕ ( n ) ( 0 ) t n n ! \bm{x}_{0}+\frac{\bm{\phi}'(0)t}{1!}+\frac{\bm{\phi}''(0)t^{2}}{2!}+\dotsb+\frac{\bm{\phi}^{(n)}(0)t^{n}}{n!} x 0 + 1 ! ϕ ′ ( 0 ) t + 2 ! ϕ ′′ ( 0 ) t 2 + ⋯ + n ! ϕ ( n ) ( 0 ) t n を近似することがわかっている．そのため，ピカール反復を適当な回数n n n [ v ] T n \lbrack\bm{v}\rbrack T^{n} [ v ] T n
p ( T ) = p n ( T ) + [ v ] T n \bm{p}(T) = \bm{p}_{n}(T)+\lbrack\bm{v}\rbrack T^{n} p ( T ) = p n ( T ) + [ v ] T n と定めれば，先述の方法がうまく機能すると期待できる．[ v ] \lbrack\bm{v}\rbrack [ v ]
おわりに
最後に，意欲的な方のために，本文では書けなかったキーワードをいくつか示しておく．
閉区間を中心と半径の組⟨ x ⟩ = ⟨ c , r ⟩ = [ c − r , c + r ] \langle x\rangle=\langle c,r\rangle=\lbrack c-r,c+r\rbrack ⟨ x ⟩ = ⟨ c , r ⟩ = [ c − r , c + r ] Yamanaka (2015) を見よ．
関数の不等式評価を得る方法はType-II PSAだけではない．平均値形式や，Makino & Berz (1999) のテイラーモデルなどがよく知られている．
常微分方程式の精度保証付き数値解法では，Moore法とLohner法が有名である．Moore法の解説はAlefeld & Mayer (2000) に，Lohner法の解説は柏木雅英（2018）にある．
この記事で説明した解法だけだと，長い区間0 ≤ t ≤ h 0\leq t\leq h 0 ≤ t ≤ h 柏木啓一郎・柏木雅英（2011） はアフィン演算という技法を応用して，長い区間にわたる解の不等式評価を計算している．
付録A（0を元に持つ区間の除算）
0 ∈ [ y ] 0\in\lbrack y\rbrack 0 ∈ [ y ] [ x ] / [ y ] \lbrack x\rbrack/\lbrack y\rbrack [ x ] / [ y ]
[ y ] = 0 \lbrack y\rbrack=0 [ y ] = 0 [ x ] / [ y ] = ∅ \lbrack x\rbrack/\lbrack y\rbrack=\emptyset [ x ] / [ y ] = ∅ [ y ] ≠ 0 \lbrack y\rbrack\neq 0 [ y ] = 0 [ x ] = 0 \lbrack x\rbrack=0 [ x ] = 0 [ x ] / [ y ] = 0 \lbrack x\rbrack/\lbrack y\rbrack=0 [ x ] / [ y ] = 0 y ▽ < 0 < y △ y^{\triangledown}\lt 0\lt y^{\vartriangle} y ▽ < 0 < y △ [ x ] ≠ 0 \lbrack x\rbrack\neq 0 [ x ] = 0 [ x ] / [ y ] = [ − ∞ , + ∞ ] \lbrack x\rbrack/\lbrack y\rbrack=\lbrack-\infty,+\infty\rbrack [ x ] / [ y ] = [ − ∞ , + ∞ ]
いずれも成り立たないとき，[ x ] / [ y ] \lbrack x\rbrack/\lbrack y\rbrack [ x ] / [ y ]
y ▽ < y △ = 0 y^{\triangledown}\lt y^{\vartriangle}=0 y ▽ < y △ = 0 0 = y ▽ < y △ 0=y^{\triangledown}\lt y^{\vartriangle} 0 = y ▽ < y △ x △ < 0 x^{\vartriangle}\lt 0 x △ < 0 [ x △ / y ▽ , + ∞ ] \lbrack x^{\vartriangle}/y^{\triangledown},+\infty\rbrack [ x △ / y ▽ , + ∞ ] [ − ∞ , x △ / y △ ] \lbrack-\infty,x^{\vartriangle}/y^{\vartriangle}\rbrack [ − ∞ , x △ / y △ ] x ▽ < x △ = 0 x^{\triangledown}\lt x^{\vartriangle}=0 x ▽ < x △ = 0 [ 0 , + ∞ ] \lbrack 0,+\infty\rbrack [ 0 , + ∞ ] [ − ∞ , 0 ] \lbrack-\infty,0\rbrack [ − ∞ , 0 ] x ▽ < 0 < x △ x^{\triangledown}\lt 0\lt x^{\vartriangle} x ▽ < 0 < x △ [ − ∞ , + ∞ ] \lbrack-\infty,+\infty\rbrack [ − ∞ , + ∞ ] [ − ∞ , + ∞ ] \lbrack-\infty,+\infty\rbrack [ − ∞ , + ∞ ] 0 = x ▽ < x △ 0=x^{\triangledown}\lt x^{\vartriangle} 0 = x ▽ < x △ [ − ∞ , 0 ] \lbrack-\infty,0\rbrack [ − ∞ , 0 ] [ 0 , + ∞ ] \lbrack 0,+\infty\rbrack [ 0 , + ∞ ] 0 < x ▽ 0\lt x^{\triangledown} 0 < x ▽ [ x ▽ / y ▽ , + ∞ ] \lbrack x^{\triangledown}/y^{\triangledown},+\infty\rbrack [ x ▽ / y ▽ , + ∞ ] [ − ∞ , x ▽ / y △ ] \lbrack-\infty,x^{\triangledown}/y^{\vartriangle}\rbrack [ − ∞ , x ▽ / y △ ]
y ▽ < 0 < y △ y^{\triangledown}\lt 0\lt y^{\vartriangle} y ▽ < 0 < y △ x △ < 0 or 0 < x ▽ x^{\vartriangle}\lt 0\;\mathrel{\textrm{or}}\;0\lt x^{\triangledown} x △ < 0 or 0 < x ▽
R = { u v ∣ u ∈ [ x ] and v ∈ [ y ] ∖ { 0 } } R = \biggl\lbrace\frac{u}{v}\biggm\vert u\in\lbrack x\rbrack\;\mathrel{\textrm{and}}\;v\in\lbrack y\rbrack\setminus\lbrace 0\rbrace\biggr\rbrace R = { v u u ∈ [ x ] and v ∈ [ y ] ∖ { 0 } } は閉区間の和集合で
R = [ − ∞ , a ] ∪ [ b , + ∞ ] = { [ − ∞ , x △ / y △ ] ∪ [ x △ / y ▽ , + ∞ ] ( x △ < 0 ) , [ − ∞ , x ▽ / y ▽ ] ∪ [ x ▽ / y △ , + ∞ ] ( x ▽ > 0 ) \begin{aligned}
R &= \lbrack-\infty,a\rbrack\cup\lbrack b,+\infty\rbrack\\
&= \begin{cases}\lbrack-\infty,x^{\vartriangle}/y^{\vartriangle}\rbrack\cup\lbrack x^{\vartriangle}/y^{\triangledown},+\infty\rbrack & (x^{\vartriangle}\lt 0),\\ \lbrack-\infty,x^{\triangledown}/y^{\triangledown}\rbrack\cup\lbrack x^{\triangledown}/y^{\vartriangle},+\infty\rbrack & (x^{\triangledown}\gt 0)\end{cases}
\end{aligned} R = [ − ∞ , a ] ∪ [ b , + ∞ ] = { [ − ∞ , x △ / y △ ] ∪ [ x △ / y ▽ , + ∞ ] [ − ∞ , x ▽ / y ▽ ] ∪ [ x ▽ / y △ , + ∞ ] ( x △ < 0 ) , ( x ▽ > 0 ) と書ける．この場合，区間演算の定義通りに計算すると
[ x ] [ y ] = hull R = [ − ∞ , + ∞ ] \frac{\lbrack x\rbrack}{\lbrack y\rbrack} = \operatorname*{hull}R
= \lbrack-\infty,+\infty\rbrack [ y ] [ x ] = hull R = [ − ∞ , + ∞ ] となる．しかし，この計算結果からは「値域が開区間( a , b ) (a,b) ( a , b ) [ − ∞ , a ] \lbrack-\infty,a\rbrack [ − ∞ , a ] [ b , + ∞ ] \lbrack b,+\infty\rbrack [ b , + ∞ ]
CAPD (version 5.2.0) のように，0 ∈ [ y ] 0\in\lbrack y\rbrack 0 ∈ [ y ] [ x ] / [ y ] \lbrack x\rbrack/\lbrack y\rbrack [ x ] / [ y ]
付録B（解の存在と一意性）
本文では証明を省略した，ピカール反復に関する定理を示す．
C ( D , R d ) C(D,\mathbb{R}^{d}) C ( D , R d )
∥ ϕ ∥ = sup t ∈ D ∣ e − 2 L t ϕ ( t ) ∣ \lVert\bm{\phi}\rVert = \sup_{t\in D}\lvert\mathrm{e}^{-2Lt}\bm{\phi}(t)\rvert ∥ ϕ ∥ = t ∈ D sup ∣ e − 2 L t ϕ ( t )∣ で定める．一様ノルム∥ ϕ ∥ ∞ = sup { ∣ ϕ ( t ) ∣ ∣ t ∈ D } \lVert\bm{\phi}\rVert_{\infty}=\sup\lbrace\lvert\bm{\phi}(t)\rvert\mid t\in D\rbrace ∥ ϕ ∥ ∞ = sup {∣ ϕ ( t )∣ ∣ t ∈ D }
e − 2 L h ∥ ϕ ∥ ∞ ≤ ∥ ϕ ∥ ≤ ∥ ϕ ∥ ∞ \mathrm{e}^{-2Lh}\lVert\bm{\phi}\rVert_{\infty} \leq \lVert\bm{\phi}\rVert \leq \lVert\bm{\phi}\rVert_{\infty} e − 2 L h ∥ ϕ ∥ ∞ ≤ ∥ ϕ ∥ ≤ ∥ ϕ ∥ ∞ が成り立つので，C ( D , R d ) C(D,\mathbb{R}^{d}) C ( D , R d ) ∥ ∙ ∥ \lVert\mathord{\bullet}\rVert ∥ ∙ ∥ K K K P f P_{\bm{f}} P f ( K , ∥ ∙ ∥ ) (K,\lVert\mathord{\bullet}\rVert) ( K , ∥ ∙ ∥) ϕ 0 , ϕ 1 ∈ K \bm{\phi}_{0},\bm{\phi}_{1}\in K ϕ 0 , ϕ 1 ∈ K
∥ P f ( ϕ 1 ) − P f ( ϕ 0 ) ∥ ≤ sup t ∈ D ( e − 2 L t ∫ 0 t ∣ f ( ϕ 1 ( t ′ ) , t ′ ) − f ( ϕ 0 ( t ′ ) , t ′ ) ∣ d t ′ ) ≤ sup t ∈ D ( e − 2 L t ∫ 0 t L ∣ ϕ 1 ( t ′ ) − ϕ 0 ( t ′ ) ∣ d t ′ ) , ∫ 0 t L ∣ ϕ 1 ( t ′ ) − ϕ 0 ( t ′ ) ∣ d t ′ ≤ ∫ 0 t L e 2 L t ′ ∥ ϕ 1 − ϕ 0 ∥ d t ′ = e 2 L t − 1 2 ∥ ϕ 1 − ϕ 0 ∥ \begin{gathered}
\begin{aligned}
&\lVert P_{\bm{f}}(\bm{\phi}_{1})-P_{\bm{f}}(\bm{\phi}_{0})\rVert\\
&\leq \sup_{t\in D}\biggl(\mathrm{e}^{-2Lt}\int_{0}^{t}\lvert\bm{f}(\bm{\phi}_{1}(t'),t')-\bm{f}(\bm{\phi}_{0}(t'),t')\rvert\,\mathrm{d}t'\biggr)\\
&\leq \sup_{t\in D}\biggl(\mathrm{e}^{-2Lt}\int_{0}^{t}L\lvert\bm{\phi}_{1}(t')-\bm{\phi}_{0}(t')\rvert\,\mathrm{d}t'\biggr),
\end{aligned}\\
\begin{aligned}
\int_{0}^{t}L\lvert\bm{\phi}_{1}(t')-\bm{\phi}_{0}(t')\rvert\,\mathrm{d}t' &\leq \int_{0}^{t}L\mathrm{e}^{2Lt'}\lVert\bm{\phi}_{1}-\bm{\phi}_{0}\rVert\,\mathrm{d}t'\\
&= \frac{\mathrm{e}^{2Lt}-1}{2}\lVert\bm{\phi}_{1}-\bm{\phi}_{0}\rVert
\end{aligned}
\end{gathered} ∥ P f ( ϕ 1 ) − P f ( ϕ 0 )∥ ≤ t ∈ D sup ( e − 2 L t ∫ 0 t ∣ f ( ϕ 1 ( t ′ ) , t ′ ) − f ( ϕ 0 ( t ′ ) , t ′ )∣ d t ′ ) ≤ t ∈ D sup ( e − 2 L t ∫ 0 t L ∣ ϕ 1 ( t ′ ) − ϕ 0 ( t ′ )∣ d t ′ ) , ∫ 0 t L ∣ ϕ 1 ( t ′ ) − ϕ 0 ( t ′ )∣ d t ′ ≤ ∫ 0 t L e 2 L t ′ ∥ ϕ 1 − ϕ 0 ∥ d t ′ = 2 e 2 L t − 1 ∥ ϕ 1 − ϕ 0 ∥ より
∥ P f ( ϕ 1 ) − P f ( ϕ 0 ) ∥ ≤ sup t ∈ D ( 1 − e − 2 L t 2 ) ∥ ϕ 1 − ϕ 0 ∥ ≤ 1 2 ∥ ϕ 1 − ϕ 0 ∥ \begin{aligned}
\lVert P_{\bm{f}}(\bm{\phi}_{1})-P_{\bm{f}}(\bm{\phi}_{0})\rVert &\leq \sup_{t\in D}\biggl(\frac{1-\mathrm{e}^{-2Lt}}{2}\biggr)\lVert\bm{\phi}_{1}-\bm{\phi}_{0}\rVert\\
&\leq \frac{1}{2}\lVert\bm{\phi}_{1}-\bm{\phi}_{0}\rVert
\end{aligned} ∥ P f ( ϕ 1 ) − P f ( ϕ 0 )∥ ≤ t ∈ D sup ( 2 1 − e − 2 L t ) ∥ ϕ 1 − ϕ 0 ∥ ≤ 2 1 ∥ ϕ 1 − ϕ 0 ∥ だから，P f P_{\bm{f}} P f
参考文献
Alefeld, Götz; Mayer, Günter. Interval analysis: theory and applications. J. Comput. Appl. Math . 2000, vol. 121, p. 421-464. https://doi.org/10.1016/S0377-0427(00)00342-3 , (accessed 2023-12-11).
Boldo, Sylvie et al. Floating-point arithmetic. Acta Numerica . 2023, vol. 32, p. 203-290. https://doi.org/10.1017/S0962492922000101 , (accessed 2023-12-05).
Kapela, Tomasz et al. CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems. Commun. Nonlinear Sci. Numer. Simul . 2021, vol. 101, 105578.
柏木啓一郎, 柏木雅英. 平均値形式とアフィン演算を用いた常微分方程式の精度保証法. 日本応用数理学会論文誌. 2011, vol. 21, no. 1, p. 37-58, (online), 入手先, J-STAGE , (参照 2023-12-15).
Kashiwagi, Masahide. “Numerical Validation for Ordinary Differential Equations Using Power Series Arithmetic”. Proc. 1994 Symposium on Nonlinear Theory and its Applications . Kagoshima, Japan, 1994-10-06/08, IEICE, 1994, p. 243-246.
柏木雅英. “常微分方程式の精度保証付き数値解法”. 精度保証付き数値計算の基礎. 大石進一編. コロナ社, 2018, p. 165-196.
柏木雅英. kv. version 0.4.55, 2022-09-16. http://verifiedby.me/kv/ , (参照 2023-12-06).
Makino, Kyoko; Berz, Martin. Efficient Control of the Dependency Problem Based on Taylor Model Methods. Reliable Computing . 1999, vol. 5, p. 3-12, (online), available from SpringerLink , (accessed 2023-12-15).
荻田武史, 柏木雅英, 劉雪峰. “序論”. 精度保証付き数値計算の基礎. 大石進一編. コロナ社, 2018, p. 11-32.
Yamanaka, Naoya; Oishi, Shin’ichi. Interval Arithmetic and Its Implementations. RIMS Kôkyûroku Bessatsu . 2015, B54, p. 71-98. http://hdl.handle.net/2433/241301 , (accessed 2023-12-05).
IEEE Std 754: 2008. Standard for Floating-Point Arithmetic .
