クラフチックの方法による非線形方程式の解の精度保証と高速化

(1)

平成 16 年度

卒業論文

クラフチックの方法による非線形方程式の解の精度保証と高速化

平成 17 年 2 月 2 日

指導教授 : 大石進一 ^教授

早稲田大学理工学部情報学科

1G01P087-9 平野晃

(2)

第 1 _章

序論

(6)

1.1 _背景

今日,いろいろな現象をモデル化し,その現象を微分方程式などの基礎方程式で解析する研究が進められてきた.したがって,微分方程式を解いたり,積分を計算したりするなどの連続数学の問題を解くことの需要が非常に高い.ところが,連続数学の問題を解くことは困難である場合が多い.そこで,計算機を用いて連続数学の問題を解こうとする研究が進められてきた.数値計算の理論もそういった研究から作られてきた.このような数値計算による連続数学の問題へのアプローチは,現代科学技術の大幅な進歩を促してきた.

従来,近似計算は近似解を求めることのみに使用されることが多く,近似解の近傍に真の解が存在するかといったことや,近似解と真の解の間の誤差評価などは,数値計算をする人の経験による判断や常識によって判断されてきた部分があった．これは,数値計算によって数学的に正しいことが保証された結論を得るためには計算のコストがあまりにもかかり過ぎて現実的ではなかったためである．しかし,こういった経験による判断や常識が通用しない,例外的な数値計算がある場合に関しては何ら保証されていない.

そこで近年,計算機の能力が大幅に向上したので,数値計算を近似計算のみにとどめておくのではなく,数値計算によって連続数学の問題に対しても数学的に正しいことが保証された結論を導くための理論と技術の開発が活発に進められている.これを精度保証数値計算という.計算機上での浮動小数点数演算における丸め誤差を含む演算結果の保証が,理論的にも実用的にも高い精度で効率よく実現できることが明らかにされ,計算の信頼性という問題が世界的にも注目を浴びている.

また,数値計算アルゴリズムにも影響を与え,今後の科学技術計算のあるべきひとつの方法として考えられるようになっている.

1.2 _{本論文の目的}

丸め誤差の厳密な評価を行う上で現在主流となっている手法に区間演算と呼ばれるものがある.この区間演算を用いて非線形方程式の解を求める場合,反復法による近似解が広く

(7)

用いられている.この際,反復法の過程における計算を区間演算に置き換えても,近似解の含まれる区間が得られるだけであり,真の解の精度保証には何も役立たない.そこで,与えられた方程式をこれと同値な不動点問題に置き換え,不動点定理の成立条件を区間演算を用いて数値的に評価することにより,解の存在保証と誤差評価を行うのである. 本論文では,非線形方程式に対する精度保証付き数値計算法の中で,非線形方程式に対する解の存在証明の代表的なものとしてクラフチック法を利用している.この方法により非線形方程式の解の精度保証付き数値計算を,従来の方法から計算量を減らし,高速かつ容易に行うことが本論文の目的である.

1.3 _{本論文の構成}

本論文の構成は以下の通りである.

第 2 章では,精度保証付き数値計算の準備として,浮動小数点数システム,区間演算,精度保証の基本的概念について述べる.

第 3 章では,ニュートン法,クラフチック法を用いた非線形方程式の解の精度保証付き数値計算について述べる.

第 4 章では,クラフチック法の高速化について述べる.

第 5 章では,数値例とMATLABでの実行結果,データについて述べる.

第 6 章では,統括として結果を元に考察について述べる.

(8)

第 2 _章

準備

(9)

2.1 _{浮動小数点数}

現代の数値計算では,浮動小数点数が標準として用いられている. 浮動小数点数は現在では標準化が進められ,ほとんどのパソコンやワークステーション, ベクトル計算機,並列計算機などで同じ形式の浮動小数点体系が用いられている. この体系に基づいて,高速に浮動小数点演算が実行される回路がコンピュータに実装されているため, 数値計算は浮動小数点数体系上で行われることがほとんどである. そこで,このような浮動小数点数の体系の上での計算の精度を保証することが重要である.

そのために,まず浮動小数点数について説明する. 浮動小数点数については IEEE 標準 754 ( IEEE standard 754 ,以下 IEEE 754 と略記する. ) がパソコンやワークステーションなどをはじめとして,多くのコンピュータで標準的に用いられている. 本論文では,以下

IEEE 754 に基づく 2 進数浮動小数点数システムを考えることにする.

2.1.1 浮動小数点数

IEEE 754では4つのタイプの数が用意されている. それは規格化2進浮動小数点数,零,

非規格化 2 進浮動小数点数, NaN ( Not a Number , 非数) である.

2 進規格化浮動小数点数

2 進規格化浮動小数点数とは a=±

Ã1 2 +d2

2² + d3

2³ +· · ·+dN

2^N

!

×2^e, (d_i = 0または1). (2.1) と書ける数をいう. emin を負の整数,emax を正の整数として,e は emin ≤e≤emax となる整数である.

m = 1 2 +d₂

2² +d₃

2³ +· · ·+d_N

2^N. (2.2)

を符合付き仮数 (signed mantissa) といい, e を指数 (exponent) という. 指数 e も 2 進数で表される.通常,単精度,倍精度 (8 byte = 64 bit) ,拡張倍精度 (10 byte =

80 bit) の浮動小数点システムがあるが,それぞれつぎのような浮動小数点システムで

(10)

ある.

N = 24, (−126 ≤e≤127), N = 53, (−1022≤e≤1023), N = 64, (−16382≤e≤16383).

規格化 2進浮動小数点システムにおいて表される数の絶対値の最大値は x_max=

µ1 2 + 1

2² + 1

2³ +· · ·+ 1 2^N

¶

2^e^max. (2.3)

であり,その最小値は

xmin = 1

22^e^min. (2.4)

である.倍精度数では(2.3),(2.4)はそれぞれ

(2−2⁵²)×2¹⁰²³ = 1.7976931348623157· · · ×10³⁰⁸, (2.5)

2⁻¹⁰²³ = 2.22507358507201· · · ×10⁻³⁰⁸. (2.6)

である. |x|> x_max のときにオーバーフロー (overflow) が生じたという.

倍精度浮動小数点数においては,仮数部が 53bit である.

2⁻⁵³= 1.1102230246· · · ×10⁻¹⁶. (2.7) より,倍精度浮動小数点数は 10進数で約 16桁の精度がある.

零

零は規格化されて

+

µ0 2+ 0

2² + 0

2³ +· · ·+ 0 2^N

¶

2^e^min. (2.8)

と表される.

非規格化 2 進浮動小数点数

IEEE 754 では浮動小数点数は指数部がe_min となったとき,仮数部の最初の桁が 1よ

(11)

り小さい数を表すために,デフォルトで最初の桁を1とすることをやめ,ここが0となる数を置くことを許す規格となっている. これを非規格化数(denormalized number) という. 非規格化数の範囲に数が入ることを,漸近アンダーフロー(gradual underflow) という. このような非規格化浮動小数点数の正で最も小さな数は

2⁻¹⁰⁷⁴ = 4.94065645841246544· · · ×10⁻³²⁴. (2.9)

である.これ以下の数になると,アンダーフロー (underflow) が生じたという.

NaN

このほかに, IEEE 754 ではつぎのような特別な数が用意されている.

(a) NaN(Not a Number)は √

−5,∞

∞,+∞+ (−∞) など不当な演算の結果として得

られる.

(b) ±∞ はオーバーフローの結果や零で割った結果として得られる.

(c) ±0はアンダーフローか ±∞での割り算の結果として得られる.

2.1.2 浮動小数点数と丸め

IEEE 754 では,つぎの4 つの丸めのモードが指定できる. c を実数 (c∈R)とする.

上向きの丸め (round upward)

c 以上の浮動小数点数の中で最も小さい数に丸める. これを 4 :R → F と表す.アルゴリズムでの表記は UP とする.

下向きの丸め (round downward)

c 以下の浮動小数点数の中で最も大きい数に丸める.これを 5 : R →F と表す. アルゴリズムでの表記は DOWNとする.

最近点への丸め (round to nearest)

c に最も近い浮動小数点数に丸める.これを □: R →F と表す. アルゴリズムでの

(12)

表記は NEAR とする.もし,このような点が 2 点ある場合には,仮数部の最後のビットが偶数である浮動小数点数に丸める. これを偶数丸め方式(round to even)という.

切り捨て (round toward 0)

絶対値が c 以下の浮動小数点数の中で, c に最も近いものに丸める.アルゴリズムでの表記は ZERO とする.

2.1.3 浮動小数点数の性質

丸めの演算を写像として°:R→F と書く. すなわち,° は 4,5,□ のいずれかと考える. IEEE 754 では,丸めの演算はつぎの条件を満たす.

°x=x, (任意のx∈Fについて),

x <=y→ °x <=°y, (任意のx, y ∈Rについて). (2.10) また, x∈F のとき,符号を変えることにより, −x や |x| が得られるので,これらは正確に計算される.

IEEE 754 では,つぎの性質が成立する.

□(−x) = −□x, (任意のx∈Rについて), 4(−x) = − 5x, (任意のx∈Rについて),

5(−x) = − 4x, (任意のx∈Rについて). (2.11)

IEEE 754では,浮動小数点数演算(F 上での四則演算)は丸めとの関係により,つぎのよ

うに定義されている. · ∈ {+,−,×, /} , ° ∈ {4,5,□} のとき

x°· y=°(x·y), (任意のx, y ∈Rについて). (2.12) この式は,左辺の浮動小数点数の四則演算の結果 x°· y は,右辺の数学的に正しい (実数としての) 四則演算の結果 x·y を指定された丸めを行って得られた数 °(x·y) に一致するように計算することを表している.

(13)

また,平方根も

(√

x)f p=°(√

x), (任意のx∈Fについて). (2.13) と,浮動小数点数演算によって計算された平方根 (√

x)_{f p} は, 正確な実数演算で計算された平方根 √

x を指定された丸めの方向へ丸めた数となる. 注意すべきことは,指数関数や三角関数などはこのような規格を満たしていない. つまり,初等関数の値を精度保証付きで求めるためには工夫が必要である.

2.2 _区間解析

精度保証付き数値計算の原理は, 実数値で与えられる真の解の上限と下限を浮動小数点演算により計算しようとするものであった．このような考え方は区間解析と呼ばれる数学的理論と密接な関係がある．この区間解析についてまず簡単に述べる．

区間解析では区間を数の拡張と考える．そして, その四則演算が定義される．これを区間演算という．区間演算を浮動小数点システムの上で展開するものを機械区間演算という．

以下,区間演算について述べていく．

2.3 _区間演算

機械区間演算を定義するのが目標であるが,まず, 実数上での演算を仮定する区間演算の概念を説明するここでの議論において四則演算は実数の四則演算とする．

2.3.1 四則演算の定義

区間解析においては,区間とは

[x, x] = {x∈R| x≤x≤x} (2.14)

(14)

と表される閉区間であるとする．ただし,x ≤ x ∈ Rで,それぞれの区間の下端,上端とする．以下,記号の節約のため,

[x] = [x, x] (2.15)

と区間を表すことにし,区間[x]と単に書けば,それは[x] = [x, x]を表すものとする．x=x となる区間[x]は点区間という．点区間は実数であるので,点区間[x]の表す実数をxと書くことにする．

区間[x]について以下を定義する．

d([x]) = x−x, r([x]) = x−x

2 , m([x]) = x+x

2 (2.16)

d([x]), r([x]), m([x])をそれぞれ区間[x]の直径,半径,中心という．

区間[x]の最小絶対値と最大絶対値をそれぞれ次で定義する．

h[x]i= min{|x| |x∈[x]}

|[x]|= max{|x| | x∈[x]}= max{|x|,|x|} (2.17) 2つの区間[x],[y]の距離を

q([x]) = max{|x−y|,|x−y|} (2.18) で定義する.

2つの区間[x],[y]が与えられたときその二つの区間の四則演算を次で定義する．

[x]◦[y] = {x◦y | x∈[x], y ∈[y]} (2.19) ただし,◦ ∈ {+,−,×, /}．これを区間演算という．この定義では区間演算は無限回の実数演 算をしないと実行できないような印象を与えるが,実は次が成立する．

[x] + [y] = [x+y, x+y]

[x]−[y] = [x−y, x−y]

[x]×[y] = [min{xy, xy, xy, xy},max{xy, xy, xy, xy}]

[x]/[y] = [x]×

"

1 y,1

y

#

, (0∈/ [y]) (2.20)

(15)

これから区間演算が有限回の実数演算で実行できることがわかる．かけ算(Table 2.1)と

割算(Table 2.2)においては,場合分けにより,より少ない手間で計算する次のような計算規

則も簡単に導ける．

Table 2.1: 区間[x],[y]の乗算[x][y]

y≥0 y30 y≤0

x≥0 [xy, xy] [xy, xy] [xy, xy]

x30 [xy, xy] [min{xy, xy},max{xy, xy}] [xy, xy]

x≤0 [xy, xy] [xy, xy] [xy, xy]

Table 2.2: 区間[x],[y]の除算[x]/[y]

y≥0 y≤0 x≥0 [x/y, x/y] [x/y, x/y]

x30 [x/y, x/y] [x/y, x/y]

x≤0 [x/y, x/y] [x/y, x/y]

2.3.2 区間演算の性質

区間演算については,包含関係における単調性

[x]⊆[x⁰], [y]⊆[y⁰]⇒[x]◦[y]⊆[x⁰]◦[y⁰], ◦ ∈ {+,−,×, /} (2.21) が成立する．

また,加法と乗法に関し交換則と結合則が成立する．

[x]◦[y] = [y]◦[x], ◦ ∈ {+,×}

[x]◦([y]◦[z]) = ([x]◦[y])◦[z], ◦ ∈ {+,×} (2.22)

(16)

しかし,加法と乗法の逆元は存在しない．すなわち,−[x] = [−x,−x]であるが 0 = [0]⊆[x]−[x] = [x−x, x−x]

1 = [1]⊆[x]/[x] (2.23)

となることがわかる．上式で等号は[x]が点区間のときのみ成立する．

また,分配則も区間演算に対しては成立しない．そのかわり次の劣分配則が成立する．

[x]×([y] + [z])⊆[x]×[y] + [x]×[z] (2.24) この式で等号は例えば区間[y]と[z]が同じ符合をもつときに成立する．

2.3.3 区間拡張

関数f :D ⊂R→Rを領域Dの任意の閉区間の上で連続な初等関数とする．関数fを次のようにしてIR上の関数に拡張することができる．

f([x]) ={f(x) | x∈[x]} (2.25)

IRを実数上の有界閉区間の集合とする．関数f :D⊂R→Rが初等関数の合成や四則演算に基づく演算規則で定義されているとき,fをIR上の関数に拡張したい．しかし,f([x]) を厳密に計算することは非常な計算時間を要することがある．そこで,fの演算規則を単純に区間演算で置き換えた演算規則によって定義される関数をfの区間拡張といいf_□で表す．

f([x])⊆f_□([x]) (2.26)

が成立する．区間拡張は一意的に定まるわけではなく,fを定義する数式を同等な他の表現形式に書き直したとき, それらの表現に基づく区間拡張が一般に異なる評価を与えることに注意する．

上式で等号が成立するとき,厳密な包み込みができたという．厳密な包み込みは一般に難しいが,f の数学的な表現式に区間値パラメータが現れず,また, [x]が一度しかその表現に

(17)

現れなければ,その表現に基づく区間拡張が厳密な包み込みを与えることが知られている．

一般に式の表現に[x]が少なく現れれば現れるほど,区間包囲は厳密な包み込みに近くなることが多いことが経験されている．区間[x]の幅d([x])が小さいとき,

q(f([x]), f_□([x]))≤αd([x]) (2.27) となることが知られている．ただし, αは正の[x]によらない定数である．

なお,一般に関数の値域f([a, b])をf([a, b])⊂[c, d]と区間[c, d]で評価するとき,区間[c, d]

のことをf([a, b])の区間包囲という．区間拡張は区間包囲の一種である．

一般の整式の区間包囲は区間を変数として与え計算したもの,すなわち区間拡張したものを用いればよい. 初等関数のうちe^xや√

xなどの単調関数に関しても,区間の両端の計算によって区間包囲が求められる.

2.4 _{機械区間演算}

前節で説明した区間演算においては,区間の両端の数は実数であるとし,厳密な実数演算に基づいて区間演算が定義された．しかし,浮動小数点数システム上で数値計算するときには,厳密な区間演算は実行できない．そこで,数学的な正解が含まれる区間が演算の結果得られることを念頭に区間演算で与えられた浮動小数点数システムR上に展開する方法を説明する．

IR={[x]∈IR | x, x∈R} (2.28)

とする．IRの要素を機械区間という．機械区間は区間の定義で上端と下端の数が浮動小数点数で与えられるものをいう．区間[a, b] ∈ IRを区間[c, d] ∈ IRに丸める区間の丸め

♦:IR→IRを

[x]⊂ ♦[x] (2.29)

および,次を満たすように定義する．

♦[x] = [x] (すべての[x]∈IRについて)

(18)

[x]⊂[y]⇒ ♦[x]⊂ ♦[y] (すべての[x],[y]∈IRについて)

♦(−[x]) = −♦[x] (すべての[x]∈IRについて) (2.30)

IEEE754上の浮動小数点数システムにおいては

♦[x] = [5x,4x] (2.31)

と定義すれば上記の条件は満たされる．機械区間演算の四則演算を

[x]♦[y] = ♦([x]◦[y]) (2.32)

と定義する．ただし,◦ ∈ {+,−,×, /}である．

具体的には機械区間演算は次のように定義される．

[x] + [y] = [5(x+y),4(x+y)]

[x]−[y] = [5(x−y),4(x−y)]

[x]∗[y] = [5(min{xy, xy, xy, xy),4(max{xy, xy, xy, xy)]

[x]/[y] = [5(min{x/y, x/y, x/y, x/y),4(max{x/y, x/y, x/y, x/y)] (2.33) 解を有限回の四則演算の組合せで求める方法(直接解法という)において, 浮動小数点数を機械区間に置き換え,四則演算を機械区間演算に置き換えて計算をすれば,真の解を含む区間が得られる．こうして,数値計算により数学的に正しい結論を導き出すために, 直接解法が知られている問題については, このような置き換えを行なうことが盛んに行なわれた．

ガウスの消去法において,途中の計算における数を機械区間に,四則演算を機械区間演算に置き換えたアルゴリズムを区間ガウス消去法という．区間ガウス消去法はこのようなアルゴリズムの代表的なものである．

区間ガウス消去法では, 実行の途中で現れる区間の幅の爆発が生じてしまうことが多い．

したがって,区間ガウス消去法は,高々,数十次元の線形方程式の解法としてのみ有効であると考えられている．

(19)

第 3 _章

非線形方程式の解の精度保証付き数値計算

(20)

3.1 _{ニュートン法}

3.1.1 はじめに

一般に非線形方程式は直接的に解を求めることは不可能であり，その近似解は反復法により求められる．ニュートン法はその中でももっとも基本的な非線形方程式の反復解法である．そこで，ニュートン法について説明する．

3.1.2 ニュートン法

f :Rⁿ →Rⁿを1回連続微分可能な関数とする．このとき，方程式

f(x) = 0 (3.1)

を解くことを考える．x₀ ∈Rⁿをf(x) = 0の近似解とするが,f(x₀) = 0とする．このとき x₁ =x₀+hがf(x) = 0のよりよい近似解となるようにhを定めることを考える．||h|| →0 のとき

f(x₁) =f(x₀ +h) = f(x₀) +f⁰(x₀)h+o(||h||) (3.2) であるから,f(x₀+h)をf(x₀) +f⁰(x₀)hと近似し

f(x₀) +Df(x₀)h= 0 (3.3)

によってhを定める．ただし,|h| → 0のとき,o(h)/h → 0とする．f(x₀) = 0 という仮定より,

h=−f(x0)

f⁰(x₀) (3.4)

となる．よって

x₁ =x₀− f(x₀)

f⁰(x₀) (3.5)

となる．ここで

g(x) =x− f(x)

f⁰(x) (3.6)

(21)

とおくと上式は

x1 =g(x0) (3.7)

と書け,これを繰り返して

f(xn+1) =g(xn) =xn−f⁰x(xn)⁻¹f(xn)(n = 0,1,2,· · ·) (3.8) によって解を改良していく方法をニュートン法という．ただし,f⁰(x)はf の点x におけるヤコビ行列である．

3.1.3 簡易ニュートン法

ニュートン法の変形として,次のようなものも考えられる．G∈R^n×nとして

g(x) = x−Gf(x) (3.9)

と置き,近似解法として

xn+1 =g(xn) (n = 0,1,2,· · ·) (3.10) を用いる方法である．これを簡易ニュートン法という．

G=f⁰(x₀)⁻¹ (3.11)

とする場合や,Gをf⁰(x₀)⁻¹の何らかの近似とする場合が多い．

3.1.4 簡易ニュートン法収束定理

fを一回連続微分可能な関数とする．方程式f(x) = 0を考える．このとき,その近似解 x0に対し,あるδ >0 が存在して,次の条件を満たすとする．

[x0−δ, x0+δ]∈I (3.12)

maxx∈[x0−δ, x0]|1−rf⁰(x)|=k <1 (3.13)

|rf(x₀)|+kδ ≤δ (3.14)

(22)

このとき,f(x) = 0の真の解x∗は[x₀−δ, x₀+δ]内にただ１つ存在し,誤差評価

|x₀−x∗ | ≤ |rf(x₀)|

1−k (3.15)

が成り立つ．

3.2 不動点の存在を保証する理論

非線形方程式の解の精度保証付き数値計算は,その解を不動点とする不動点問題に帰着される.そこで不動点の存在を保証するために縮小写像と縮小写像原理,そしてBrouwerの定理を説明する.

3.2.1 縮小写像原理

数値計算の原理を考える際の基本的方法は, 問題を不動点定理に帰着させることである．

今,Xをバナッハ空間,gをXからその中への写像とする．

g(x) = x (3.16)

を満たすXの点を不動点という．写像gによって動かしても動かない点という意味である．

写像gの不動点を求める問題を不動点問題という．

写像gに不動点が存在することを保証する定理を不動点定理という．以下で述べる縮小写像原理は不動点定理の中でも, もっとも基本的かつ有用なものである．また,不動点の逐次反復による求め方も示されており, 具体的な解法に直接結び付く場合も多い．

写像gの定義域をD(g),値域をR(g)とする．それぞれバナッハ空間Xの部分集合とする．このとき,写像g :D(g) →R(g)が縮小写像であるとは, ある定数k, (0≤ k <1)が存在して, 任意のx₁, x₂ ∈D(g)に対して

kg(x₁)−g(x₂)k ≤kkx₁−x₂k (3.17) が成立することである．

このとき,次の定理が成立する：

(23)

定理 3.1 (バナッハの縮小写像の原理)

Xをバナッハ空間とし,Sをその空でない閉部分集合とする．g :S →Sを縮小写像とすると, Sの中にgの不動点が存在する．しかもこの不動点は唯一である．

さらに,Sの中の任意の点x₀をとり,逐次反復(iteration)

x_n+1 =g(x_n), n = 0,1,2,· · · (3.18) を考えると,{x_n}^∞_n=1はSの中の点列となる. Sの中のgの不動点をx^∗とするとき,x_nはx^∗ へ収束し, かつ事前誤差評価(a priori error estimate)

kxn−x^∗k ≤ kⁿ

1−kkx0−g(x0)k (3.19)

が成立する.また,事後誤差評価(a posteriori error estimate) kx_n+1−x^∗k ≤ k

1−kkx_n+1−x_nk (3.20)

も成立する. □

数値計算の精度を保証する有効な原理の一つは, 問題を適当なバナッハ空間X上の写像g の不動点を求める問題に帰着させることである.すなわち, 適切な閉集合Sを選び,その上でgが縮小写像となること,および,

g(S)⊂S (3.21)

が満たされることを数値計算により示そうとするものである. Sが十分小さな閉集合であれば,解の存在のみならず,シャープな誤差評価も得られたことになる.

縮小写像の原理の応用として次の定理が成立する.

定理 3.2 (バナッハの摂動定理)

X, Y をバナッハ空間, A: X → Y を有界線形作用素で有界な逆作用素A⁻¹ :Y →Xをもつとする. このとき,有界線形作用素B :X →Y が条件

kA⁻¹Bk_L(X,Y₎ ≤1 (3.22)

(24)

を満たすとき,有界線形作用素A+B : X → Y は有界な逆作用素(A+B)⁻¹ : Y → Xをもち

k(A+B)⁻¹kL(X,Y) ≤ kA⁻¹k_L(X,Y₎

1− kA⁻¹Bk_L(X,Y₎ (3.23)

となる. □

この定理は,例えば,A, Bをn×n行列とすると, 与えられた行列T =A+Bの近似逆行列をA⁻¹とすると,T⁻¹の存在が逆行列の計算(3.22)に帰着させられることがわかる. これから,前節までの議論により, 精度保証付きでこの条件の成立がチェックできることになる.

定理 3.3 (Brouwerの定理)

Dが有解閉凸集合で、写像g :D →Dを連続写像とする.このとき

g(D)⊂D (3.24)

なら,D内に少なくとも一つは不動点が存在する. □ これらの定理から不動点の存在を示すことが出来るので,精度保証付き数値計算が可能となる.

(25)

3.3 _{クラフチック法}

簡易ニュートン法の収束定理の条件をより容易に,かつ自然な形で評価する方法として, クラフチック(Krawczyk)法がある.これは,ニュートン法の区間版であり,縮小写像の原理と一体化されて真の解を含む区間を生成する方法である.

3.3.1 クラフチック法

簡易ニュートン法の収束定理の条件を容易かつ,自然な形で評価する手法.クラフチック法はニュートン法を区間解析的に行ったものである.

3.3.2 平均値形式

平均値形式を用いると,区間演算における区間幅の増大を抑えることができる.X =Rⁿ, Y = R^m, D ⊂Xとする.以下,I(D)からI(Y)への写像を区間写像と言う.ここで,m :I(D)→Y をm(I)∈Iを満たす任意の関数とする.f :D→Y に対して,

F(I) = f(c) +F⁰(I)(I−c), c =M(I) (3.25) で定義された区間写像F :I(D)→I(Y)を中心形式と言う.とくにm(I) = Mid(I)のとき, 平均値形式と言う.

3.3.3 精度保証

f :Rⁿ →Rⁿを1回連続微分可能な非線形写像として

f(x) = 0 (3.26)

を考える.Lをf⁰(x) の近似行列,R をLの逆行列とする.Rと cは通常の浮動小数点演算

(Newton法)により近似的に求めたとしよう.このとき,f(x) = 0の近似解cの近くに真の解

x∗が存在するか否かや,存在する場合には||x^∗−c||を評価することを考える.

T =c+X,(X = [−e, e], e= 2|Rf(c)|) (3.27)

(26)

とする.クラフチック作用素K(X+c)−c=Hは

H =−Rf(c) + (E−RF⁰(X+c))X (3.28)

で与えられる.H ⊂ X ならばc + X 内にf(x) = 0 の真の解が存在する.F⁰(X + c) ⊂ (F_c⁰, F_w⁰), f(x)∈(f_c, f_w)とする.ただし,f_c, F_c⁰は中心,f_w, F_w⁰ は半径とする.

setround(down);

M =E+ (−R)F⁰_c+ (−|R|)F⁰_w; r = (−R)f_c+ (−|R|)f_w;

setround(up);

M =E+ (−R)F⁰c+|R|F⁰w; r = (−R)f_c+|R|f_w;

M_max =max|M|,|M|;

q =Mmaxe;

K =r+q;

setround(down);

K =r−q;

このとき,H = [K, K]で与えられる.

3.3.4 計算量

この計算方法によるクラフチックの方法の計算には逆行列を求めるのにn³+O(n²), M , M の計算にそれぞれ2n³ +O(n²)の計算量を必要とする.その他の計算にはO(n²)の計算量しか必要としないので,合計してこの精度保証付き数値計算には5n³+O(n²)の計算量がかかる.

(27)

第 4 _章

高速化

(28)

4.1 _はじめに

ここでクラフチック作用素の計算量に注目する．非線形項の計算量は与えられた問題によって異なる.また,現在のところ,非線形項の計算は区間演算による方法が,細かい差異はあっても,ほぼ唯一の方法であるので,これに要する計算量はどんなアルゴリズムであってもほぼ同一であると考える。

前述のHの式を直接用いるアルゴリズムでは,（行列）×（行列）の計算が行われているため,計算量がO(n³)となる．そこでその計算量を減らすためにHに対して,以下のような作業をしてやる．

4.2 高速化理論

クラフチック作用素Hに対して,次のような作用素LH を考える．

LH =L(−L⁻¹f(c) + (E−L⁻¹F⁰(c+X))X) (4.1)

=−f(c) +L(E−L⁻¹F⁰(c+X))X (4.2)

=−f(c) + (L−F⁰(c+X))X (4.3)

これより,区間ベクトルV を

V =−f(c) + (L−F⁰(X+c))X (4.4)

とおくと,Hは線形連立方程式

LH =V (4.5)

を解くことにより求められる．上記のV =−f(c) + (L−F⁰(X+c))Xより,区間ベクトル V を評価するには,まずL−F⁰(c+X)の評価が必要になる.この部分は,次式の右辺のように計算する事により,精度保証付きで評価できる.

L−F⁰(c+X)⊂[L−F_c⁰−F_r⁰, L−F_c⁰ +F_r⁰] (4.6)

(29)

これより,M を

M =maxⁿ|L−F_c⁰−F_r⁰|,|L−F_c⁰+F_r⁰|^o (4.7) とおくことにより,

(L−F⁰(c+X))X ⊂[−Me, Me] (4.8)

と評価できる.よって,区間ベクトルV は,

V ⊂[−f_c−f_r−Me,−f_c+f_r+Me] (4.9) によって精度保証付きで評価される.

区間ベクトルV を求めるアルゴリズムは以下のようになる.

setround(down);

M =L−F⁰c−F⁰r; setround(up);

M =L−F⁰_c+F⁰_w; Mmax =max{|M|,|M|};

p=M_maxe;

q =f_w+p;

V =−fc+q;

setround(down);

V =−f_c−q;

そして、

V ⊂[V , V] (4.10)

が成り立つ．そして,右辺が区間ベクトルである線形連立方程式

LH = [V , V] (4.11)

(30)

の解をすべて含む区間を[H, H]とすると

H ⊂[H, H] (4.12)

と評価される．

4.3 計算量

この高速化されたクラフチック法による精度保証付き数値計算では,LH = [V , V]の計算に ²₃n³+O(n²)の浮動小数点数乗算でできる.区間V を求めるアルゴリズムではO(n²)回の浮動小数点数乗算で実行できることから,²₃n³+O(n²)の計算量でクラフチック作用素H を評価できることが出来る.

よって,クラフチックの方法を普通に計算するときと比較すると、高速に数値計算が可能になっている.

次章では,具体的な数値例を解いた結果を算出し,高速化理論と普通に計算する場合との比較を行う.

(31)

第 5 _章

数値シミュレーション

(32)

5.1 _はじめに

本章では具体的な数値例を用いて,実際に非線型方程式の解の存在検証をMATLABを用いて行い,その有効性を検証する.そして,クラフチック法と高速化したプログラムでの実行時間と解の精度の比較を行う.

5.1.1 MATLAB

MATLABは高品質で用意にプログラムが組める数値計算用語であり,高速に数値計算を

実行する命令が組み込まれており,安定性,高速性に評判がある.

また,INTLABというMATLABで書かれたツールを組み込んでい用いれば,浮動小数点の丸めモードを制御しながら区間演算や微分計算なども高速に,そして容易に実現してくれる.

非線形方程式の解の精度保証付き数値計算をするときには,非線型方程式の評価が必要となるので,区間演算は不可避である.MATLABにINTLABを組み込めば,区間演算が容易になると同時に高速な行列計算も行う.したがって,非線型方程式の解の数値計算には適したツールであるといえる.

MATLABの計算環境は次の通りである.

Table 5.1: MATLABの実行環境

OS Windows 2000

CPU x86 Family 6 Model 7 Stepping2 Memory 256MB

Version MATLAB version 6

(33)

5.2 _数値例

本節では、非線型方程式f(x) = 0に具体的に定め,その解の精度保証付き数値計算を実行する. ここで,xをnベクトルとしx= (x₁, x₂,· · ·, x_n)で表し,f(x) = 0を

f(x) = 0 ⇐⇒











f1(x) = 0 f₂(x) = 0

...

f_n(x) = 0

(5.1)

で表される非線型方程式とする.

今回は数値例として,エサキダイオードからなる非線形抵抗回路に関する非線形方程式を考える.この回路方程式は,

f₁(x) = g(x₁) +x₁+x₂+· · ·+x_n−1 (5.2) f₂(x) = g(x₂) +x₁+x₂+· · ·+x_n−n (5.3)

... (5.4)

f_n(x) = g(x_n) +x₁+x₂+· · ·+x_n−n (5.5) で示される.また,エサキダイオードのI−V 特性を

g(x) = 2.5x³−10.5x²+ 11.8x (5.6) とする.このf(x)のヤコビ行列f⁰(x)を計算すると,

f⁰(x) =







g⁰(x₁) + 1 1 · · · · 1 1 g⁰(x₂) + 1 1 · · · 1

...

1 1 · · · · g⁰(x_n) + 1







(5.7)

となる.

(34)

5.3 _解の検証

まず,解の存在検証に成功する例を具体的に示す.n = 10の場合について考える. ニュートン法により,この非線型方程式の近似解cを１つ求めると,

c=







−0.33545308322972

−0.28285123145438

−0.22613660635450

−0.16436462311202

−0.09617922680485

−0.01951290678153 0.06906744069606 2.20298819542304 2.47107520698319 2.61563562075326







(5.8)

となる.この近似解cをもとに,解の存在性を検証する区間X = [−e, e]を計算すると,

X = 10⁻¹³∗







·

−0.00454601242579, 0.00454601242579

¸

·

−0.00488028934477, 0.00488028934477

¸

·

−0.01000685757526, 0.01000685757526

¸

·

−0.00579154161660, 0.00579154161660

¸

·

−0.00644417607327, 0.00644417607327

¸

·

−0.00732881462610, 0.00732881462610

¸

·

−0.02401227947181, 0.02401227947181

¸

·

−0.28409380300438, 0.28409380300438

¸

·

−0.12780338838763, 0.12780338838763

¸

·

−0.32352059743798, 0.32352059743798

¸







(5.9)

(35)

となる.このとき,e= 2|L⁻¹f(c)|として計算している. また,f_c, f_w, F_c⁰, F_w⁰ をeを元に計算しておき,これよりHを計算すると,

H = 10⁻¹⁴∗







·

−0.03937785617337, 0.03841287862725

¸

·

−0.04697745186037, 0.04109864592197

¸

·

−0.04310792343714, 0.04723131363711

¸

·

−0.05006209005331, 0.04883272526193

¸

·

−0.05564246569900, 0.05427456683918

¸

·

−0.06318687745085, 0.06163119724334

¸

·

−0.07378411355387, 0.08050693231758

¸

·

−0.82178034322933, 0.72020895302417

¸

·

−0.46909553122972, 0.38790919297676

¸

·

−0.12659634242161, 0.41732137285859

¸







(5.10)

となり,H ⊂Xが満たされることから区間T 内に解が存在する.

5.4 計算時間の比較

次にnを100から1000まで値を変えながら,クラフチック作用素Hの計算にかかる時間, 解の精度を測定し,以下の表にまとめた.解の精度はeの値を任意で決め,解が包み込まれた際のeの最大値ノルムを計算した.

(36)

Table 5.2: 数値例のMATLABにおける精度保証付き数値計算の結果クラフチックの方法高速化したプログラム

次元実行時間(sec) 解の精度実行時間(sec) 解の精度

100 0.090 1.168∗10⁻¹² 0.134 1.168∗10⁻¹²

200 0.50 4.102∗10⁻¹¹ 0.49 4.102∗10⁻¹¹

300 1.35 1.586∗10⁻¹¹ 1.21 1.586∗10⁻¹¹

400 2.84 6.273∗10⁻¹² 2.34 6.273∗10⁻¹²

500 5.14 1.908∗10⁻¹¹ 4.04 1.908∗10⁻¹¹

600 8.54 6.294∗10⁻¹¹ 6.57 6.294∗10⁻¹¹

700 13.4 2.742∗10⁻¹¹ 10.1 2.742∗10⁻¹¹

800 19.6 5.662∗10⁻¹¹ 14.3 5.662∗10⁻¹¹

900 27.7 2.714∗10⁻¹⁰ 19.9 2.714∗10⁻¹⁰

1000 38.1 4.142∗10⁻⁸ 27.7 4.142∗10⁻⁸

(37)

5.5 _考察

本論文の結論として,今回実施した高速化プログラムはクラフチック法と比較して速度の点で有効な手段であるということが言えた.大幅な速度の改善とまではいかなかったが,解く非線形方程式の解の次元が大きくなればなるほど,理論的に速度の差が大きくなる.

解の精度については両方の場合でほとんど差はない.これは,近似点をニュートン法から求め,その結果を利用して解の存在範囲を調べる区間を計算している.したがって,精度に差が見られないのも妥当な結果といえる.

以上の結果から大規模な計算を行う時には今回の方法を利用する価値があると言える.

(38)

第 6 _章

統括と今後の課題

(39)

6.1 _統括

本論文では,MATLABによる非線形方程式の精度保証付き数値計算の高速化について述べた.

第 2 章では,準備として浮動小数点,区間演算について述べた．

第 3 章では,非線形方程式の解法,ニュートン法とクラフチック法について述べた．

第 4 章では,クラフチック法の高速化について述べた．

第 5 章では,数値例と実行結果,また実行結果を元に考察を述べた．

本論文の構成は以上のようであった．

6.2 今後の課題

今回はある近似点を求めて,その周りに任意の区間を与えてその範囲内に解が存在するのかを検証した.そのため,二つの方法には解の精度に違いは見られなかった. 今後の課題として,近似点を求めず,ある区間のみを与えた場合にどのような差が見られるかを検証できるようにしたい. また,今回の数値例ではヤコビ行列f⁰(x)を求める際に,単純な形になったのでn次元のf⁰(x)を単純なループで作ることが出来た.しかし,複雑な関数が出てきた場合,全てを手計算で作ることは難しいし,時間の無駄である.そこで,関数を計算するプログラムから関数の微分をプログラムを計算するプログラムを自動生成できれば望ましい.これを自動微分という.この自動微分を用いれば,複雑な関数におけるニュートン法が可能であり,たくさんの例を実験することが出来たと思う.

(40)

6.3 _むすび

本論文では非線形方程式の精度保証数値計算の高速化について論じた.高速化についてある程度の結果を出すことに成功した.

(41)

謝辞

(42)

この卒業論文の完成にあたってたくさんの方々から多くの助言や指導を頂きました. まず,本研究を進めるに当たり,終始丁寧な御指導及び御激励を賜り,最後まで面倒を見て下さいました大石進一教授に深く感謝いたします.

大変お忙しい中でも私の質問に対して熱心に答えて下さいました, 理工学研究科客員講師,丸山晃佐氏,中谷祐介氏,荻田武史氏に大変感謝いたします.

研究室内や合宿など,日常的な場面でお世話になりました,大石研究室博士課程1年,尾崎克久氏, 大石研究室修士課程2年,大山博毅氏,森山敦史氏, 並びに大石研究室修士課程1 年,坂内太郎氏,島本誠氏,関本竜平氏,徳永克久氏,福地健氏,細田幸裕氏,横山大輔氏に深く感謝いたします.

研究内容だけでなくあらゆる場面で意見の交換や協力などをして下さいました, 大石研究室4年,有滝貴広氏,岩田揚平氏,大成高顕氏,川崎文裕氏,栗原崇氏,佐藤友規氏,ト詣各氏,水島裕氏,吉岡毅氏,吉田昂央氏には様々な助言をもらい大変感謝しています.

(43)

参考文献

[1] 大石進一,数値計算,裳華房,1999年,

[2] 大石進一,精度保証付き数値計算,コロナ社,2001年, [3] 大石進一,MATLABによる数値計算,培風館,2002年,

[4] 中谷祐介大石進一,非線形方程式に対する精度保証付き数値,日本シミュレーション学会第20回シミュレーションテクノロジーコンファレンス,pp261-264,2001年6月

クラフチックの方法による非線形方程式の 解の精度保証と高速化

平成 16 年度

卒業論文