平成 16 年度
卒業論文
コンピュータ代数システムによる
ガレルキン近似の数式処理及び数値計算の検証
平成１７年２月２日
指導教授 : 大石 進一 教授
早稲田大学 理工学部 情報学科
1G01P017-7 岩田 揚平
目 次
1 序論 3
1.1 背景 . . . 4
1.2 本論文の目的 . . . 5
1.3 本論文の構成 . . . 5
2 準備 6 2.1 はじめに . . . 7
2.2 浮動小数点数 . . . 7
2.2.1 浮動小数点数 . . . 7
2.2.2 浮動小数点数と丸め . . . 9
2.2.3 浮動小数点数の性質 . . . 10
2.3 区間解析 . . . 11
2.4 区間演算 . . . 11
2.5 機械区間演算 . . . 15
2.6 むすび . . . 17
3 関数解析 18 3.1 はじめに . . . 19
3.2 関数解析 . . . 19
3.3 バナッハ空間と作用素 . . . 19
3.4 縮小写像原理 . . . 23
3.5 ニュートン法 . . . 25
3.6 むすび . . . 26
4 ガレルキン法 27 4.1 はじめに . . . 28
4.2 ガレルキン法の概要 . . . 28
4.3 ガレルキン近似解の例 . . . 31
4.4 むすび . . . 31
5 実行結果 32 5.1 はじめに . . . 33
5.2 実行環境 . . . 33
5.2.1 実行環境 . . . 33
5.2.2 MuPAD . . . 33
5.3 実行内容 . . . 34
5.4 実行結果 . . . 34
5.5 検証 . . . 35
5.6 近似解 . . . 37
5.7 むすび . . . 37
6 統括 38 6.1 はじめに . . . 39
6.2 統括 . . . 39
6.3 考察 . . . 39
6.4 むすび . . . 40
謝辞 41
参考文献 43
第 1 章
序論
1.1 背景
今日,いろいろな現象をモデル化し,その現象を微分方程式などの基礎方程式で解析する 研究が進められてきた. したがって,微分方程式を解いたり,積分を計算したりするなどの 連続数学の問題を解くことの需要が非常に高い. ところが,連続数学の問題を解くことは困 難である場合が多い.そこで,計算機を用いて連続数学の問題を解こうとする研究が進めら れてきた.数値計算の理論はそういった研究から作られてきた.このような数値計算による 連続数学の問題へのアプローチは,現代科学技術の大幅な進歩を促してきた.
従来,近似計算は近似解を求めることのみに使用されることが多く,近似解の近傍に真の 解が存在するかといったことや, 近似解と真の解の間の誤差評価などは,数値計算をする人 の経験による判断や常識によって判断されてきた部分があった. これは,数値計算によって 数学的に正しいことが保証された結論を得るためには計算のコストがあまりにもかかり過 ぎて現実的ではなかったためである. しかし,こういった経験による判断や常識が通用しな い,例外的な数値計算がある場合に関しては何ら保証されていない.
そこで近年,計算機の能力が大幅に向上したので,数値計算を近似計算のみにとどめてお くのではなく, 数値計算によって連続数学の問題に対しても数学的に正しいことが保証さ れた結論を導くための理論と技術の開発が活発に進められている.これを精度保証付き数値 計算という. 計算機上での浮動小数点数演算における丸め誤差を含む演算結果の保証が,理 論的にも実用的にも高い精度で効率よく実現できることが明らかにされ,計算の信頼性とい う問題が世界的にも注目を浴びている.
また,この数値計算の他にコンピュータで数学的問題を解くための方法として,数式処理 とも呼ばれるコンピュータ代数がある. コンピュータ代数は,数を任意精度で表し,かつ計 算をコンピュータで記号的に扱うものであり,記号的あるいは代数的な解法が知られていな い問題においては数値計算による近似的な解法を用いる必要があるが, 記号的な計算が行 われるために計算結果は厳密となる.
これらは今後の科学技術計算においてあるべき方法として考えられている.
1.2 本論文の目的
ダフィング方程式の周期解をガレルキン法を用いて近似する際,非常に複雑な数式処理及 び数値計算が必要となる. 本論文ではこれを代数処理システムを用いて行い,その性能を検 証する.
1.3 本論文の構成
本論文の構成は以下の通りである.
第２章では,数値計算の基礎となる,浮動小数点数システム,区間演算,機械区間演算につ いて述べる.
第３章では,非線形現象の解析に役立つ関数解析の基礎及びニュートン法について述べる.
第４章では,非線形微分方程式の周期解を求めるガレルキン法について述べる.
第５章では,コンピュータ代数システムを用いてガレルキン法の決定方程式を導出する際 の扱うデータ量と実行時間の分布を検証し,また数値計算により近似解を求める.
第６章では,本論文の統括を述べる.
第 2 章
準備
2.1 はじめに
本章では,数値計算の基礎となる,浮動小数点システム,区間演算,機械区間演算について 述べる.
2.2 浮動小数点数
現代の数値計算では,浮動小数点数が標準として用いられている.浮動小数点数は現在で は標準化が進められ,ほとんどのパソコンやワークステーション,ベクトル計算機,並列計算 機などで同じ形式の浮動小数点体系が用いられている.この体系に基づいて,高速に浮動小 数点演算が実行される回路がコンピュータに実装されているため,数値計算は浮動小数点数 体系上で行われることがほとんどである.
そのため,まず浮動小数点数について説明する.浮動小数点数については IEEE 標準754
( IEEE standard 754 ,以下 IEEE 754 と略記する ) がパソコンやワークステーションな
どをはじめとして,多くのコンピュータで標準的に用いられている.本論文では,以下 IEEE 754 に基づく 2 進数浮動小数点数システムを考えることにする.
2.2.1 浮動小数点数
IEEE 754では 4 つのタイプの数が用意されている.それは規格化 2進浮動小数点数,零,
非規格化 2 進浮動小数点数, NaN ( Not a Number , 非数) である.
2 進規格化浮動小数点数
2 進規格化浮動小数点数とは a=±
Ã1 2 +d2
22 + d3
23 +· · ·+dN
2N
!
×2e, (di = 0または1). (2.1) と書ける数をいう. emin を負の整数,emax を正の整数として,e は emin ≤e≤emax と なる整数である.
m = 1 2 +d2
22 +d3
23 +· · ·+dN
2N. (2.2)
を符合付き仮数(signed mantissa)といい,eを指数(exponent) という.指数eも2進 数で表される.通常,単精度,倍精度(8 byte = 64 bit) ,拡張倍精度(10 byte = 80 bit) の浮動小数点システムがあるが,それぞれつぎのような浮動小数点システムである.
N = 24, (−126 ≤e≤127), N = 53, (−1022≤e≤1023), N = 64, (−16382≤e≤16383).
規格化 2進浮動小数点システムにおいて表される数の絶対値の最大値は xmax=
µ1 2 + 1
22 + 1
23 +· · ·+ 1 2N
¶
2emax. (2.3)
であり,その最小値は
xmin = 1
22emin. (2.4)
である.倍精度数では(2.3),(2.4)はそれぞれ
(2−252)×21023 = 1.7976931348623157· · · ×10308, (2.5) 2−1023 = 2.22507358507201· · · ×10−308. (2.6) である. |x|> xmax のときにオーバーフロー (overflow) が生じたという.
倍精度浮動小数点数においては,仮数部が 53bit である.
2−53= 1.1102230246· · · ×10−16. (2.7) より,倍精度浮動小数点数は 10進数で約 16桁の精度がある.
零
零は規格化されて
+
µ0 2+ 0
22 + 0
23 +· · ·+ 0 2N
¶
2emin. (2.8)
と表される.
非規格化 2 進浮動小数点数
IEEE 754 では浮動小数点数は指数部がemin となったとき,仮数部の最初の桁が 1よ
り小さい数を表すために,デフォルトで最初の桁を1とすることをやめ,ここが 0とな る数を置くことを許す規格となっている.これを非規格化数 (denormalized number) という.非規格化数の範囲に数が入ることを,漸近アンダーフロー(gradual underflow) という.
このような非規格化浮動小数点数の正で最も小さな数は
2−1074 = 4.94065645841246544· · · ×10−324. (2.9) である.これ以下の数になると,アンダーフロー (underflow) が生じたという.
NaN
このほかに, IEEE 754 ではつぎのような特別な数が用意されている.
(a) NaN(Not a Number)は √
−5,∞
∞,+∞+ (−∞) など不当な演算の結果として得
られる.
(b) ±∞ はオーバーフローの結果や零で割った結果として得られる.
(c) ±0はアンダーフローか ±∞での割り算の結果として得られる.
2.2.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 と表す.アルゴリズムでの表 記は NEAR とする.もし,このような点が 2 点ある場合には,仮数部の最後のビット が偶数である浮動小数点数に丸める.これを偶数丸め方式 (round to even) という.
切り捨て (round toward 0)
絶対値が c 以下の浮動小数点数の中で, c に最も近いものに丸める.アルゴリズムで の表記は ZERO とする.
2.2.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) に一致する ように計算することを表している.
また,平方根も
(√
x)f p=°(√
x), (任意のx∈Fについて). (2.13) と,浮動小数点数演算によって計算された平方根 (√
x)f p は,正確な実数演算で計算された 平方根 √
xを指定された丸めの方向へ丸めた数となる.注意すべきことは,指数関数や三角 関数などはこのような規格を満たしていない.つまり,初等関数の値を精度保証付きで求め るためには工夫が必要である.
2.3 区間解析
前述の精度保証付き数値計算は,実数値で与えられる真の解の上限と下限を浮動小数点演 算により計算しようとするものである. このような考え方は区間解析と呼ばれる数学的理 論と密接な関係がある.この区間解析についてまず簡単に述べる.
区間解析では区間を数の拡張と考える.そして,その四則演算が定義される.これを区間演 算という. 区間演算を浮動小数点システムの上で展開するものを機械区間演算という. ま た,区間演算および機械区間演算はそれを基本演算として直接的に考えると実用上大きな問 題を生じることが多い. このような問題点を解決するために前述の精度保証付き数値計算 の原理が考えられた.
以下,区間演算について述べていく.
2.4 区間演算
機械区間演算を定義するのが目標であるが，まず，実数上での演算を仮定する区間演算 の概念を説明するここでの議論において四則演算は実数の四則演算とする．
区間解析においては，区間とは
[x, x] = {x∈R| x≤x≤x} (2.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)
二つの区間[x],[y]が与えられたときその二つの区間の四則演算を次で定義する．
[x]◦[y] = {x◦y | x∈[x], y ∈[y]} (2.18) ただし，◦ ∈ {+,−,×, /}．これを区間演算という．この定義では区間演算は無限回の実数 演算をしないと実行できないような印象を与えるが，実は次が成立する．
[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.19)
これから区間演算が有限回の実数演算で実行できることがわかる．かけ算(Table2.1)と
割算(Table2.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]
区間演算については，包含関係における単調性
[x]⊆[x0], [y]⊆[y0]⇒[x]◦[y]⊆[x0]◦[y0], ◦ ∈ {+,−,×, /} (2.20) が成立する．
また，加法と乗法に関し交換則と結合則が成立する．
[x]◦[y] = [y]◦[x], ◦ ∈ {+,×}
[x]◦([y]◦[z]) = ([x]◦[y])◦[z], ◦ ∈ {+,×} (2.21)
しかし，加法と乗法の逆元は存在しない．すなわち，−[x] = [−x,−x]であるが 0 = [0]⊆[x]−[x] = [x−x, x−x]
1 = [1]⊆[x]/[x] (2.22)
となることがわかる．上式で等号は[x]が点区間のときのみ成立する．
また，分配則も区間演算に対しては成立しない．そのかわり次の劣分配則が成立する．
[x]×([y] + [z])⊆[x]×[y] + [x]×[z] (2.23) この式で等号は例えば区間[y]と[z]が同じ符合をもつときに成立する．
関数f :D ⊂R→Rを領域Dの任意の閉区間の上で連続な初等関数とする．関数fを 次のようにしてIR上の関数に拡張することができる．
f([x]) ={f(x) | x∈[x]} (2.24)
IRを実数上の有界閉区間の集合とする．関数f : D ⊂ R → Rが初等関数の合成や四 則演算に基づく演算規則で定義されているとき，f をIR上の関数に拡張したい．しかし，
f([x])を厳密に計算することは非常な計算時間を要することがある．そこで，fの演算規則
を単純に区間演算で置き換えた演算規則によって定義される関数をfの区間拡張といいf□ で表す．
f([x])⊆f□([x]) (2.25)
が成立する．区間拡張は一意的に定まるわけではなく，f を定義する数式を同等な他の表 現形式に書き直したとき，それらの表現に基づく区間拡張が一般に異なる評価を与えるこ とに注意する．
上式で等号が成立するとき，厳密な包み込みができたという．厳密な包み込みは一般に 難しいが，fの数学的な表現式に区間値パラメータが現れず，また，[x]が一度しかその表 現に現れなければ，その表現に基づく区間拡張が厳密な包み込みを与えることが知られて
いる．一般に式の表現に[x]が少なく現れれば現れるほど，区間包囲は厳密な包み込みに近 くなることが多いことが経験されている．区間[x]の幅d([x])が小さいとき，
q(f([x]), f□([x]))≤αd([x]) (2.26) となることが知られている．ただし，αは正の[x]によらない定数である．
なお，一般に関数の値域f([a, b])をf([a, b]) ⊂ [c, d]と区間[c, d]で評価するとき，区間
[c, d]のことをf([a, b])の区間包囲という．区間拡張は区間包囲の一種である．
2.5 機械区間演算
前節で説明した区間演算においては，区間の両端の数は実数であるとし，厳密な実数演 算に基づいて区間演算が定義された．しかし，浮動小数点数システム上で数値計算すると きには，厳密な区間演算は実行できない．そこで，数学的な正解が含まれる区間が演算の 結果得られることを念頭に区間演算で与えられた浮動小数点数システムR上に展開する方 法を説明する．
IR={[x]∈IR | x, x∈R} (2.27)
とする．IRの要素を機械区間という．機械区間は区間の定義で上端と下端の数が浮動小 数点数で与えられるものをいう．区間[a, b] ∈ IRを区間[c, d] ∈ IRに丸める区間の丸め
♦:IR→IRを
[x]⊂ ♦[x] (2.28)
および，次を満たすように定義する．
♦[x] = [x] (すべての[x]∈IRについて)
[x]⊂[y]⇒ ♦[x]⊂ ♦[y] (すべての[x],[y]∈IRについて)
♦(−[x]) = −♦[x] (すべての[x]∈IRについて) (2.29)
IEEE754上の浮動小数点数システムにおいては
♦[x] = [5x,4x] (2.30)
と定義すれば上記の条件は満たされる．機械区間演算の四則演算を
[x]♦[y] = ♦([x]◦[y]) (2.31)
と定義する．ただし，◦ ∈ {+,−,×, /}である．
具体的には機械区間演算は次のように定義される．
[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.32) 解を有限回の四則演算の組合せで求める方法(直接解法という)において，浮動小数点数 を機械区間に置き換え，四則演算を機械区間演算に置き換えて計算をすれば，真の解を含 む区間が得られる．こうして,数値計算により数学的に正しい結論を導き出すために，直接 解法が知られている問題については, このような置き換えを行なうことが盛んに行なわれ た．ガウスの消去法において，途中の計算における数を機械区間に，四則演算を機械区間 演算に置き換えたアルゴリズムを区間ガウス消去法という．区間ガウス消去法はこのよう なアルゴリズムの代表的なものである．
このようなアプローチはある程度の成功を収めたが，反面，機械区間演算を演算の単位 と考えると，いろいろな問題が生じることもわかった．その一つは，このような計算では 丸めの切替が頻繁におき，計算速度の極端な低下が見られることや，プログラムが複雑に なってしまうことである．また，例えば，区間ガウス消去法では，実行の途中で現れる区 間の幅の爆発が生じてしまうことが多い．したがって，区間ガウス消去法は，高々，数十 次元の線形方程式の解法としてのみ有効であると考えられている．
本論文では，このような点の反省に基づき，区間演算を基本単位とは考えず，内積計算 程度のもっと大きな単位を計算の基礎とすることで精度保証をすることを論じていく．
2.6 むすび
本章では,浮動小数点システム,区間演算,機械区間演算の概説を行った. 次章では,非線 形現象を解析するのに役立つ,関数解析について説明する.
第 3 章
関数解析
3.1 はじめに
本章では,非線形現象を解析するのに役立つ,関数解析および,非線形方程式の解を求める 問題に対して有力なニュートン法について説明する.
3.2 関数解析
内積計算や行列の積の数値計算結果などは，区間演算によって精度保証できる．しかし，
線形代数方程式の解や，非線形方程式の解の存在を示し，その近似解を精度保証するため には工夫が必要である．そのために考え出されたのが関数解析における不動点定理(縮小写 像原理)を利用するというアイディアである．
以下,非線形現象の解析のために重要となる関数解析の基礎について述べる．
3.3 バナッハ空間と作用素
まず，Rを実数の集合とし，XをR上のベクトル空間とする．ベクトル空間のことを線 形空間ともいう．
有界閉区間[a, b]上の実数値連続関数の全体をC[a, b]と書く．C[a, b]の任意の元に対し て和を(x+y)(t) = x(t) +y(t), t ∈ [a, b]で定義する．また，実数αによるスカラー倍を (αx)(t) = αx(t), t ∈[a, b]で定義する．このとき，C[a, b]はベクトル空間となる．C[a, b]の ように関数のつくるベクトル空間を関数空間という．
ベクトル空間Xの要素となる任意のベクトルxに対して，非負の実数値kxkが対応し，
これが次の(N1)から (N3)までの条件を満たすとき，これをノルムと呼ぶ．ノルムはベク トルの長さを拡張した概念である．
(N1) kxk ≥0, kxk= 0とx= 0は同値 (N2) kαxk=|α|kxk
(N3) kx+yk ≤ kxk+kyk ノルムは次の不等式を満たす．
|kxk − kyk| ≤ kx±yk ≤ kxk+kyk (3.1) 同じベクトル空間の上でもノルムは幾通りも定義できることがある．例えば，n次元ベ クトル空間の上ではx= (x1, x2,· · ·, xn)に対して，ユークリッドノルムkxk =
vu utXn
i=1
x2i が 普通よく用いられるが
kxk= max
1≤i≤n|xi| (3.2)
もまたノルムとなる．これを最大値ノルムという．以下では，ノルム空間Xのノルムであ ることを明示することが必要であるときには，x ∈ Xに対しkxkX と，どのノルム空間の ノルムであるかを下つきで表すことにする．
ノルム空間X内の点列{xn}∞n=1が点x∈Xに収束すること lim
n→∞xn =xを
n→∞lim =kxn−xk= 0 (3.3)
によって定義する．式 (3.3) を，xn →x, (n→ ∞)とも書く．
ノルム空間Xの部分集合Mが有界であるとは，ある非負定数cが存在して，任意のx∈M に対して
kxk ≤c (3.4)
が成立することである．
連続関数の空間C[a, b]の任意の元xに対してノルムを kxk= max
a≤t≤b|x(t)| (3.5)
で定義する．これは，ノルムの条件を満たす．このノルムを最大値ノルムといい．他のノ ルムと区別するときにはkxk∞と書く．
ノルム空間X の無限点列{xn}∞n=1が基本列またはコーシー列であるとは，任意の実数 ε >0に対して，ある正の整数Nが存在して，任意の整数m, n≥N に対して
kxm−xnk< ε (3.6)
となることである．これは実数の基本列の定義で，絶対値をノルムで置き換えたものである．
ノルム空間Xがバナッハ空間であるとは，その任意の基本列{xn}∞n=1に対して，ある x∗ ∈Xが唯一つ定まり
kxn−x∗k →0, (n → ∞) (3.7)
が成立することである．任意の基本列が収束列になるとき，そのノルム空間は完備である ともいう．完備ノルム空間をバナッハ空間という．連続関数の最大値ノルムによるノルム 空間C[a, b]はバナッハ空間となる．
MとY を集合とする．Aが集合Mから集合Y への作用素であるとは，Mの各点にただ 一つのY の点yを対応させることである．この対応関係を
A:M →Y, y =Ax (3.8)
と書く．
作用素A:M →Y が 全単射のとき，逆作用素A−1 :Y →Mが
A−1y=x⇔Ax=y (3.9)
によって定義される．これは，Aが全単射であることから，任意のy∈Y に対して，ある x ∈M が存在してAx =yとなり，かつ単射であることからそのようなxはただ一つに定 まるからである．
Xの線形部分空間DからY への作用素L:D→Y が線形であるとは，次の重ねの理が 成立することをいう：
1. L(x+y) =Lx+Ly (任意のx, y ∈Dに対して) 2. L(αx) = αLx (任意のx∈Dと任意の実数αに対して)
X, Y をバナッハ空間とし，M をXの部分集合とする．作用素A:M →Y がM 上で連 続であるとは，任意のx∈M とxに収束するM の任意の点列{xn}∞n=1 ⊂Mに対し
kAx −Axk →0, (n→ ∞) (3.10)
が成立することである．これをAxn→Ax, (n→ ∞)と書く．
線形作用素L: X → Y に対し，どんなx ∈Xに対しても，xによらない非負定数Kが 存在して
kLxk ≤Kkxk (3.11)
が成立するとき，作用素Lを有界作用素という．線形作用素L:X →Y が有界作用素であ るとき
sup
x∈X,x6=0
kLxk
kxk (3.12)
が有限値となる．これをkLkと書き，線形作用素Lの作用素ノルムという．
任意のx∈Xにおいて
kLxk ≤ kLkkxk (3.13)
が成り立つ．
kLk= sup
x∈X,kxk=1
kLxk (3.14)
となる．
X, Y をバナッハ空間とする．線形作用素L:X →Y が有界作用素であることとXで連 続作用素であることは同値である．
X, Y をバナッハ空間とする．X からY への有界線形作用素の集合をL(X, Y)と書く．
L1, L2 ∈L(X, Y)とa, b∈Rに対し
L=aL1+bL2 (3.15)
も有界線形作用素となる．これはL(X, Y)がベクトル空間となることを意味している．さ らに，この空間に作用素ノルムによってノルムを導入するとL(X, Y)は作用素ノルムによ りバナッハ空間となる．
3.4 縮小写像原理
数値計算の原理を考える際の基本的方法は，問題を不動点定理に帰着させることである．
今，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)が 存在して，任意のx1, x2 ∈D(g)に対して
kg(x1)−g(x2)k ≤kkx1−x2k (3.17) が成立することである．
このとき，次の定理が成立する：
定理 3.1 (バナッハの縮小写像の原理)
Xをバナッハ空間とし，Sをその空でない閉部分集合とする．g :S → Sを縮小写像と すると，Sの中にgの浮動点が存在する．しかもこの浮動点は唯一である．
さらに，Sの中の任意の点x0をとり，逐次反復(iteration)
xn+1 =g(xn), n = 0,1,2,· · · (3.18) を考えると,{xn}∞n=1はSの中の点列となる．Sの中のgの浮動点をx∗とするとき，xnは x∗へ収束し，かつ事前誤差評価(a priori error estimate)
kxn−x∗k ≤ kn
1−kkx0−g(x0)k (3.19)
が成立する．また，事後誤差評価(a posteriori error estimate) kxn+1−x∗k ≤ k
1−kkxn+1−xnk (3.20)
も成立する． □
数値計算の精度を保証する有効な原理の一つは，問題を適当なバナッハ空間X上の写像g の不動点を求める問題に帰着させることである．すなわち，適切な閉集合Sを選び，その 上でgが縮小写像となること，および，
g(S)⊂S (3.21)
が満たされることを数値計算により示そうとするものである．Sが十分小さな閉集合であ れば，解の存在のみならず，シャープな誤差評価も得られたことになる．
縮小写像の原理の応用として次の定理が成立する．
定理 3.2 (バナッハの摂動定理)
X, Y をバナッハ空間，A:X →Y を有界線形作用素で有界な逆作用素A−1 :Y →Xを もつとする．このとき，有界線形作用素B :X →Y が条件
kA−1BkL(X,Y) ≤1 (3.22)
を満たすとき，有界線形作用素A+B :X → Y は有界な逆作用素(A+B)−1 :Y → Xを もち
k(A+B)−1kL(X,Y) ≤ kA−1kL(X,Y) 1− kA−1BkL(X,Y)
(3.23)
となる． □
この定理は，例えば，A, Bをn×n行列とすると，与えられた行列T = A+Bの近似 逆行列をA−1とすると，T−1の存在が逆行列の計算(3.22)に帰着させられることがわかる．
これから，前節までの議論により，精度保証付きでこの条件の成立がチェックできること になる．
3.5 ニュートン法
一般に非線形方程式は直接的に解を求めることは不可能であり,その近似解は反復法に より求められる.ニュートン法はその中でももっとも基本的な非線形方程式の反復解法であ る.そこで,ニュートン法について説明する.
f :RntoRnを1回連続微分可能な関数とする．このとき，方程式
f(x) = 0 (3.24)
を解くことを考える．x0 ∈Rnをf(x) = 0の近似解とするが，f(x0) = 0とする．このとき x1 =x0+hがf(x) = 0のよりよい近似解となるようにhを定めることを考える．||h|| →0 のとき
f(x1) =f(x0 +h) = f(x0) +f0(x0)h+o(||h||) (3.25) であるから，f(x0+h)をf(x0) +f0(x0)hと近似し
f(x0) +Df(x0)h= 0 (3.26)
によってhを定める．ただし，|h| →0のとき，o(h)/h→0とする．f(x0) = 0という仮定 より，
h=−f(x0)
f0(x0) (3.27)
となる．よって
x1 =x0− f(x0)
f0(x0) (3.28)
となる．ここで
g(x) =x− f(x)
f0(x) (3.29)
とおくと上式は
x1 =g(x0) (3.30)
と書け，これを繰り返して
f(xn+1) =g(xn) =xn−f0x(xn)−1f(xn)(n = 0,1,2,· · ·) (3.31) によって解を改良していく方法をニュートン法という．ただし，f0(x)はfの点xにおける ヤコビ行列である．
ニュートン法の変形として，次のようなものも考えられる．G∈Rn×nとして
g(x) = x−Gf(x) (3.32)
と置き，近似解法として
xn+1 =g(xn) (n = 0,1,2,· · ·) (3.33) を用いる方法である．これを簡易ニュートン法という．
G=f0(x0)−1 (3.34)
とする場合や，Gをf0(x0)−1の何らかの近似とする場合が多い．
3.6 むすび
本章では，関数解析の基礎およびニュートン法について説明した．次章では,非線形微分 方程式の近似解を求める方法として,ガレルキン法について述べる.
第 4 章
ガレルキン法
4.1 はじめに
本章では非線形微分方程式dxdt = f(x, t)の周期解を求める方法であるガレルキン法につ いて述べていく.
4.2 ガレルキン法の概要
非線形微分方程式
dx
dt =f(x, t) (4.1)
を考える.ただし，各tについてx(t), f(x(t), t)∈Rd で
f(x, t+ 2π) = f(x, t) (4.2)
を満たすものとする．このような周期的な微分方程式に外力に同調して現れる周期2π の 周期解を求めるガレルキン法の説明をする．
ガレルキン法は，周期解x(t) のフーリエ級数展開 x(t) = 1
2a0 +
X∞ n=1
(ancosnt+bnsinnt) (4.3) を利用する方法である．ただし，an, bn はそれぞれがd 次元のベクトルフーリエ係数であ る．周期解を求めるためには，フーリエ係数{a0, a1, b1, a2, b2,· · ·} が決まればよい．この フーリエ級数をもとの微分方程式から決定しようというのが，ガレルキン法の基本的アイ デアである．解が滑らかであればn が大きくなるとan, bn は0 に近くなる．そこで近似解 としてフーリエ級数をm 位までで打ち切った
xm(t) = 1 2a0+
Xm n=1
(ancosnt+bnsinnt) (4.4) を考える．こうすると決定すべきベクトルフーリエ係数の数が2m+ 1 個と有限個になる．
ここで、以下のために記号を用意する．
Pm :x→xm (4.5)
Pm は
x(t) = 1 2a0+
X∞ n=1
(ancosnt+bnsinnt) (4.6)
とフーリエ級数で表される関数をm 位で打ち切って xm(t) = 1
2a0+
Xm n=1
(ancosnt+bnsinnt) (4.7) と三角多項式で表される関数に写す写像である．
xm をもとの微分方程式 dxdt =f(x, t) に代入する．このときf(xm(t), t)はf が非線形 関数であれば，一般に，m 位より高位のフーリエ級数となる．そこでこれをm 位で打ち 切ってPmf(xm(t), t)とする．こうして
dxm(t)
dt =Pmf(xm(t), t) (4.8)
を満たすように係数 {a0, a1, b1, a2, b2,· · ·am, bm} を決めるのがガレルキン法である．
このプロセスを具体的に書き表す.
F0(y) = 1 π
Z 2π
0 f(xm(t), t)dt (4.9)
F2n−1(y) = 1 π
Z 2π
0 f(xm(t), t) cosntdt−nbn (4.10) F2n(y) = 1
π
Z 2π
0 f(xm(t), t) sinntdt−nan(n= 1,2,· · ·, m) (4.11) と定義する.ただし,y= (a0, a1, b1, a2, b2,· · ·am, bm)T ∈R2n+1 とする.
また,F(y) = (F0(y), F1(y),· · ·, F2m(y))T とするとF :R(2m+1)d →R(2m+1)d となる.この とき
dxm(t)
dt =Pmf(xm(t), t) (4.12)
はyを未知数とする有限次元連立非線形方程式
F(y) = 0 (4.13)
と同値となる.F(y) = 0 をガレルキン近似の決定方程式という.
周期解を求めるためのガレルキン法の数学的構造を明らかにするために,内積の記号を用 意する.
(a(t), b(t)) =
Z 2π
0 a(t)b(t)dt (4.14)
また
L(t) = dx
dt −f(x, t) (4.15)
と書く.このときガレルキン法の決定方程式は
(L(x), φn) = 0, n = 0,1,2,· · ·, M (4.16) と書かれる.ただし, φ0 = 1/2, φ1 = cost, φ2 = sint,· · ·, φ2m−1 = cosmt, φ2m = sinmtであ る.また,M = 2mである.さらに,c0 = a0, c1 = a1, c2 = b1,· · ·, c2m−1 = am, c2m = bmとす る.このとき
x=
XM n=0
cnφn (4.17)
と書ける.
一般に,解きたい微分方程式を
L(x) = 0 (4.18)
とする.余弦関数cosnt,正弦関数sinntが周期的境界条件をみたすように,所望の境界条件 をみたす一次独立な関数系φn, n= 0,1,2,· · ·によって,近似解を
x=
XM n=0
cnφn (4.19)
と仮定し,未定係数cn, n = 0,1,2,· · ·, M を
(L(x), φn) = 0, i= 0,1,2,· · ·, M (4.20) によって決定するのが一般のガレルキン法である.ただし,内積は与えられた境界値問題に 対応して適当なものを選ぶものとする.
4.3 ガレルキン近似解の例
例としてダフィングの方程式 d2x
dt2 +kdx
dt +x3 =Bcost (4.21)
を考える.
ダフィング方程式においてはx(t) が周期解となるとき, y(t) = −x(t+π) も周期解とな るという性質がある.これを対称性という.
このことから,対称周期解は x(t) =
X∞ n=1
(ancos (2n−1)t+bnsin (2n−1)t) (4.22) の形のフーリエ展開をもつことが分かる.
x(t) =acost+bsint (4.23)
を式4.21に代入するとガレルキン近似は costの係数= 0:
−a+B +kb+ 3
4(a3 +ab2) = 0 (4.24)
sintの係数= 0:
−b+ka+ 3
4(a2b+b3) = 0 (4.25)
を解いて求められる.この連立方程式の近似解は第３章で述べたニュートン法などで容易に 求めることが出来る.
4.4 むすび
本章では非線形微分方程式の周期解を求める方法であるガレルキン法について述べた. 次 章では,コンピュータ代数システムを用いて,本章で述べたガレルキン法の決定方程式を導 出する際の扱うデータ量と実行時間の分布を検証し,さらに実際に数値計算により近似解を 求める.
第 5 章
実行結果
5.1 はじめに
本章では,次節で示すMuPADを用いて, 第４章で示したガレルキン法の決定方程式を導 出する際の扱うデータ量と実行時間の分布を検証する.
さらに,実際に数値計算を行い,ガレルキン近似解を導出する.
5.2 実行環境
5.2.1 実行環境
実行環境は以下の通りである.また今回使用したコンピュータ代数システム,MuPADに ついて次の小節で述べる.
• OS：Windows XP
• CPU：Pentium4 1.70GHz
• メモリ：512MB RAM
• MuPAD：MuPAD Light 2.5.3
5.2.2 MuPAD
MuPADは,Multi-Processing Algebra Data Toolの略である.MuPADは総合的なコンピ ュータ代数システムとしてドイツのパデルボン（Paderborn）大学のフックススタイナー
（Fuchssteiner）教授のグループにより1992年以降以降開発されたものである.
今回使用したMuPAD Lightは正式版であるMuPAD Proのサブセットであり,フリーで あるため完全なGUIは提供されていないが, 数学的機能についてはMuPAD Proと同等で ある.
また,MuPADはオブジェクト指向プログラミングが可能であり,さらに数値計算システ
ムであるScilabの機能を取り込んでおり, この数値計算ツールもフリーで利用できる.
5.3 実行内容
前節で示したコンピュータ代数システム,MuPADを用いて,以下の作業を行う.
ダフィングの方程式
d2x
dt2 +kdx
dt +x3 =Bcost (5.1)
を取り上げる.
これに,対称周期解のフーリエ展開 x(t) =
X∞ n=1
(ancos (2n−1)t+bnsin (2n−1)t) (5.2)
をm位で打ち切った
xm(t) =
Xm n=1
(ancos (2n−1)t+bnsin (2n−1)t) (5.3) を代入する.これに数式処理を行いガレルキン法の決定方程式を求める.
さらに,求まったガレルキン法の決定方程式を,数値計算によって解く.
5.4 実行結果
m= 1,2,· · · 12において,ガレルキン法の決定方程式をテキストファイル形式で保存した ときの容量をTable5.1に示す.また,m = 1,2,· · · 12において,ガレルキン法の決定方程式を 求めるための実行時間をTable5.2に示す.
Table 5.1: ガレルキン法の決定方程式のファイル容量
m 1 2 3 4 5 6
ファイル容量(B) 144 609 1867 4436 8776 15389
m 7 8 9 10 11 12
ファイル容量(B) 24726 37293 53562 74718 100967 132857
Table 5.2: ガレルキン法の決定方程式導出の実行時間
m 1 2 3 4 5 6
実行時間(s) 0.321 0.551 4.867 25.62 92.47 268.3
m 7 8 9 10 11 12
実行時間(s) 661.1 1457 2895 5415 9792 15665
5.5 検証
ガレルキン法の決定方程式において,m= 1,2,· · · 12としたときの合計項数をTable5.3に
示す.また,これとTable5.1に示したガレルキン法の決定方程式のファイル容量の分布との
比較を図5.1に示す.さらに,実行時間とファイル容量の分布の比較を図5.2に示す.
Table 5.3: ガレルキン法の決定方程式の合計項数
m 1 2 3 4 5 6
合計項数 9 35 99 223 429 737
m 7 8 9 10 11 12
合計項数 1169 1747 2491 3423 4565 5937
図 5.1: ファイル容量と合計項数の比較
図 5.2: 実行時間とファイル容量の比較
図5.1より,m = 1,2,· · · 12のとき,ガレルキン法の決定方程式を格納したファイルの容 量と,ガレルキン法の決定方程式における合計項数は同程度の変化を示していることがわか る.また,図5.2より,mを大きくしたとき,ガレルキン法の決定方程式の導出時間は,求まる 決定方程式のファイル容量に比べて大きな変化を示すことがわかる.
5.6 近似解
m= 12としたときのガレルキン法の決定方程式を,k = 0.1, B = 0.15としてMuPADで 数値計算によって解き.次の近似解を得た.
x(t) =−0.15108 cost+ 0.015374 sint−0.000092085 cos 3t+ 0.000032334 sin 3t
−0.000000057291 cos 5t+ 0.000000035963 sin 5t−0.00000000003262 cos 7t + 0.000000000032538 sin 7t−0.000000000000016701 cos 9t
+ 0.000000000000026475 sin 9t−0.0000000000000000070649 cos 11t + 0.000000000000000020047 sin 11t−1.6538∗10−21cos 13t
+ 0.000000000000000000014337 sin 13t+ 1.0673∗10−24cos 15t
+ 9.7402∗10−24sin 15t+ 2.1754∗10−27cos 17t+ 6.2867∗10−27sin 17t + 2.3881∗10−30cos 19t+ 3.8335∗10−30sin 19t+ 2.1589∗10−33cos 21t + 2.1775∗10−33sin 21t+ 1.7575∗10−36cos 23t+ 1.1159∗10−36sin 23t
5.7 むすび
本章では,MuPADを用いてガレルキン法の決定方程式を導出する際の扱うデータ量と実 行時間の分布を検証し, また実際に数値計算により近似解を求めた.次章では,この結果を 元に考察を述べていく.
第 6 章
統括
6.1 はじめに
本章では,本論文の結びとして,本論文の統括を述べる.
6.2 統括
本論文では,コンピュータ代数システムを用いて,対称周期解をもつダフィング方程式に おける, ガレルキン法の決定方程式の導出の処理を行い,その性能を検証した.
第２章では,数値計算の基礎となる,浮動小数点数システム,区間演算,機械区間演算につ いて述べた.
第３章では,非線形現象の解析に役立つ関数解析の基礎及びニュートン法について述べた.
第４章では,非線形微分方程式の周期解を求めるガレルキン法について述べた.
第５章では,MuPADを用いてガレルキン法の決定方程式を導出する際の扱うデータ量と 実行時間の分布を検証し, また数値計算により,近似解を導出した.
本論文の構成は以上のようであった．
6.3 考察
第５章で検証した通り,数式処理システムにおいて複雑な数式処理を行う際,処理は記号 的に行われるため計算結果は厳密になるが, 扱うデータ量が大きくなると,計算量は非常に 大きくなる.
また,近年の計算機の能力の大幅な向上に伴い,数値計算の技術は非常に高精度で演算結 果の保証を実現するまでになっている.
これらの技術の発展は今後の理工学の分野において非常に有益なものとなる. 理工学の 分野においては,本論文で扱った,強制振動を表すモデル方程式であるダフィング方程式の ような, 様々な現象のモデル化によって得られた微分方程式を解析することが基本となり, これにおいては数式処理技術と数値計算技術双方の応用が望まれるためである.
6.4 むすび
本論文ではコンピュータ代数システムによるガレルキン近似の数式処理及び数値計算に ついて論じた. 考察を踏まえ,今後の数式処理及び数値計算技術の発展を期待したい.
謝辞
本論文の製作にあたり,たくさんの方々から多くの助言や指導を頂きました.
まず,本研究を進めるに当たり,終始丁寧な御指導及び御激励を賜り,最後まで面倒を見て 下さいました大石 進一教授に深く感謝いたします.
また,終始丁寧な御指導,御教示を下さいました, 理工学研究科客員講師,荻田 武史氏,丸 山 晃佐氏,中谷 祐介氏に深く感謝いたします.
また,研究室内や合宿など,また日常的な場面でも大変お世話になりました,大石研究室修 士課程2年,大山 博毅氏,森山 敦史氏, 並びに大石研究室修士課程1年,坂内 太郎氏,島本 誠氏,関本 竜平氏,徳永 克久氏,福地 健氏,細田 幸裕氏,横山 大輔氏に深く感謝いたします.
さらに,研究活動及びその他あらゆる場面で意見の交換や協力などをして下さいました, 大石研究室4年,有滝 貴広氏,大成 高顕氏,川崎 文裕氏,栗原 崇氏,佐藤 友規氏,平野 晃氏, ト 詣各氏,水島 裕氏,吉岡 毅氏,吉田 昂央氏に深く感謝いたします.
参考文献
[1] 大石進一：”非線形解析入門”，コロナ社，1998年，
[2] 大石進一：”精度保障付き数値計算”，コロナ社，2001年，
[3] 赤間世紀，山口喜博：”ＭｕＰＡＤで学ぶ基礎数学”，丸善株式会社，2004年
