計算機で微分方程式の厳密解を求める

大石進一(早稲田大学理工学部情報学科)


本文の構成は次のようである.
1. はじめに
2. 構成的数学
2.1 計算可能実数
2.2 区間演算と計算可能実数の四則演算
2.3 区間解析の発展
2.4 自動微分
3. 精度保証付き数値計算と構成的数学
3.1 ニュートン法と区間解析
3.2 構成的数学との関係
4. 精度保証付き数値計算で微分方程式を解く
4.1 ダフィング方程式の周期解の分岐
4.2 ホモクリニック分岐
5. むすび
参考文献

1.はじめに

(目次に戻る)

組み合わせ数学や離散数学などに現れる,整数演算の有限回の操作で正解 が得られる問題は計算機で厳密に解けるが,微分方程式の境界値問題などの 連続数学の問題は計算機では厳密に解くことはできない.連続数学の問題 が解けるとしても,記号処理による代数的な問題に限られる. このような考え方は一般に広く信じられているものではなかろうか.

これに対して,近年,精度保証付き数値計算と呼ばれる分野が発展し, 微分方程式の厳密な解を求める有効な手段を提供すること, 特に,カオス理論に対する数学的な道具として極めて有用なことが 示されつつある.従来の数値計算では,実数を有限固定桁数の浮動小数点数 で近似し,その四則演算も丸め誤差を含み近似的に計算するに過ぎなかった. また,微分方程式を差分化して扱うなど打ち切り誤差も含まれている. それでも,精度を変えて数値実験を繰り返したり、別の実験との比較をしたりして、 数値実験によって得られた近似解の近くに真の解の存在を予見し、 誤差を含む数値計算により行われた数値実験よりソリトン,カオス 等の新概念の発見を含め,非常に多くの科学的発見がなされてきた. ソリトンの研究においては、数値実験によるソリトン概念の発見に引き続いて 逆散乱法などの解析的手法が発明され、手計算により厳密解が求められる ようになった。これはソリトン系が可積分(解ける系)であることを反映 している。一方、カオスの研究においては、数値実験結果を抽象化し、 カオスを生成する力学モデルが作られ、カオス概念が普遍的であることが示された。 しかし、実際の系を記述する方程式にこのようなモデルが埋め込まれている ことを証明することは困難であった。これは、カオス系 が非可積分系であることを反映している。

精度保証付き数値計算では,数値計算に伴う誤差を厳密に評価して,数学的に 正しい結果を数値計算によって導くことを 目標にしている[1],[2].特に、十分よい近似解があったとき、 その近くに真の解が存在する十分条件を精度保証付き数値計算によって 検証することが目標となる。カオスの研究など、数値実験によって、よい近似解が 多数蓄積されている分野が多いからである。 従来,丸め誤差 の見積もりはとてつもなく過大評価となり,また,打ち切り誤差の見積もりは極めて 困難と思われていた.しかし,幾つかのブレークスルーにより,実際的な問題 に対しても,両者ともに,充分よい評価ができ,数値計算の精度が実用的に評価 できるとともに,カオス理論に現れるような微分方程式の解の存在証明も 精度保証付き数値計算によりできるようになりつつある.

本稿では,精度保証付き数値計算に至る連続数学を計算機上で展開する 考え方について応用とともに論じる.

2.構成的数学

(目次に戻る)

計算機の上では,整数や有理数はメモリが十分であれば厳密に扱うことが できる.といっても,プログラム的には多倍長数の高速演算法などかなり 高度な技術を要するが.一方,円周率(以下,\piと表す) などの無理数は,巡回しない無限小数となるので,そのすべての桁を メモリに入れることはできない.こうして,無理数は厳密には計算機に 乗せられないと考えられるようになったと思われる.しかし,そう簡単に あきらめてしまって良いのであろうか.ここでは,発想を転換して 別の方法で厳密に無理数を計算機に乗せることを考えよう.

2.1 計算可能実数

(目次に戻る)

簡単のため,を考えよう.今,計算機と計算機言語を1つ固定して いるとする.例えば,パソコンでC言語を使っているとしよう. \piの値を任意の精度まで計算する プログラムが存在することは良く知られている.すなわち,円周率が含まれる ``いれこ''の閉区間列

[an, bn], n=1,2,...;

を表す有理数an <= bnで, 区間

[an,bn]

が区間

[an+1,bn+1]

を含み,nが無限大となるときanbnが同じ極限に収束する ものを計算 するプログラムが存在する.このプログラム自身は有限長であるので, それは計算機に記憶することができる.この円周率のように その値を任意の精度で計算するプログラムが存在する実数を計算可能実数 あるいは構成的実数という.自然対数の底eなどの無理数も その値を任意精度で計算するプログラムが存在するので計算可能実数である. むしろ,自然に考えられる無理数で計算可能でないものを探すのは難しいとも いえる.理論的には1つのプログラム言語によって書ける異なるプログラムは 可算個(整数の集合と1:1対応が取れること)なので, 計算可能実数の総数も可算個となる.実数は非可算 無限個(整数の集合と1:1対応できない、より多い無限)あるので, 計算可能でない実数の方が圧倒的に多いということもできる. しかし,物理や工学に現れる方程式の係数などにこのような計算可能でない 実数が現れた例は知られていないようである.

こうして,無理数の値そのものでなく,それを計算するプログラムをその 数と思うことにすれば,計算可能な無理数は計算機の上に厳密に乗せることが できると考えることができる.このように,数学的な対象で,その値を 任意精度で計算できるプログラムが存在するとき,そのプログラムをその数学 的対象と考える(同一視する)という考え方によれば,無理数でも計算機の上に 厳密に乗せられることになる.このような考え方は,実は,数学基礎論では 比較的古くから提唱されており,マルコフ流の構成的数学と 呼ばれている[3].

本稿では,精度保証付き数値計算の基礎がこのマルコフ流の構成的数学によって 与えられることを論じて行く. しかし,数学基礎論での構成的数学は,実在する プログラムと対応する数学的対象だけに数学的に扱う対象を制限したとき, どの位のことができるかという, やや否定的な側面を扱う雰囲気が強いのに対し,精度保証付き数値計算は カオスの存在証明など,従来の数学では難しい問題を計算機を利用して 突破しようという積極的な側面が強い.このように,連続数学の問題を計算機を武器 にして解くための道具として構成的数学が利用できるというのが本稿の 新しい論点である.

2.2 区間演算と計算可能実数の四則演算

(目次に戻る)

計算機の上に無理数が小数としては厳密に表現できないというのを 改善するために区間解析[4]が考えだされた.区間解析では,円周率は区間 [3.1415,3.1416]に含まれるというように(無理)数を有理数を両端とする区間で表す. これは円周率を完全には捕らえていないが, 正しい情報は与えている.そこで,今,pは区間[a,b]に含まれること 及びqは区間[c,d]に含まれる,という ように2つの実数p,qの情報がわかっているとき,pqを四則演算 した結果はどのような区間に入るか考えてみよう.p,qに関し 付帯的な情報が なければ#を+,-,*,/のいずれかとすると

p#qは区間[a,b]#[c,d]の中に入る

と評価することになろう.ただし,

[a,b]+[c,d] = [a+c,b+d],

[a,b]-[c,d] = [a-d,b-c],

[a,b]*[c,d] = [min A, max A],

[a,b]/[c,d] = [min B,max B].

ここに,A={ac,ad,bc,bd}, B={a/c,a/d,b/c,b/d}である.上式による 区間の四則演算を区間演算という.この区間演算を利用するとpが区間 [an,bn]及びqが区間 [cn,dn]にn=1,2,...について 含まれるというように2つの計算可能実数が 与えられたとき,その2数の四則演算を

[en,fn] =[an,bn]#[cn,dn], n=1,2,...

によって定義することができる.nが無限大となるとき enfnは同じ極限に収束し,それはp#qとなって, 上式は通常の実数の四則演算の定義に一致し, 計算可能実数の集合はこの四則演算によって体(実数の集合と同じ様に 四則演算ができる集合のこと)をなすことが知られている. 区間演算のように,数値計算によって真の解が含まれている区間を計算 するような,精度が保証された数値計算を精度保証付き数値計算という. 精度保証付き数値計算の基礎となる区間演算はこのように計算可能実数の 四則演算を定義するのに直接に関連する.本稿の目標は,構成的数学の 議論に区間解析という数値計算の分野で既に30年以上の歴史をもち,高度に 発展した技法を持ち込むこむと,どうなるかの議論ということもできる. 計算可能実数の集合は体の構造をもち,通常の理工学のモデリングで現れる 数は計算可能実数であるから,理工学の問題を解くとき,計算可能実数体の 上で議論を展開してもよさそうである.従来,数学基礎論側からの議論では, このような立場から,積極的に新しい応用ができるというというものは 少なかったようであるが,ここではそこそこの議論はできそうであるという 積極論を展開してみたい.

2.3 区間解析の発展

(目次に戻る)

ガウスの消去法など数値計算における数値演算を区間演算に置き換えると, 連立方程式の解を含む区間が得られるなど,結果は数学的に正しいものが 得られる.しかし,区間演算を繰り返すと,得られる区間の幅は,すぐに 増大し,結果として数学的に正しい結論が得られるとしても,答えは -百万から+百万の間にあるといった,意味のない結果しか得られない ことが多い.また,ニュートン法などのような反復解法では,数値演算を 区間演算で置き換えただけでは,数学的に正しい結果は得られない.ニュートン 法の途中で得られる近似解と真の解の間の関係は一般には不明であるから である.

このような困難から,区間解析は一時その研究者の数を減らした.このよう なとき,区間幅を減少させるテクニックが幾つか開発され,現代区間解析が誕生 したのである.そのようなブレークスルーの一つは平均値形式という 考え方である.今,fを実数の上で定義された1回微分可能な 実数値関数とする.平均値形式は、写像fによる 区間[a,b]の像f([a,b]) を過大にならずに評価することを目的に考え出された。 f'(x)f(x)xに関する微分として

F([a,b])=f(c)+f'([a,b])([a,b]-c)

とする.ただし、c=(a+b)/2.これを,fの平均値形式 という.平均値の定理から f([a,b])F([a,b])に含まれることがわかる. この平均値形式を用いると,区間演算の欠点である 区間幅の増大が抑えられることがある.例えば,f(x)=x-xとしよう. 明らかに,f=0であるから,f([a,b])=0である. 区間演算では f([a,b])[a,b]-[a,b]= [a-b,b-a]に含まれる と評価していまい, もとの区間[a,b]の2倍の幅の区間として評価される. 一方,平均値形式を用いると,f'(x)=1-1=0であるから

F([a,b])=0+0*[a-c,b-c]=0

となる.こうして,f([a,b])F([a,b])に含まれる ことよりf([a,b])=0が示され、 区間演算では決して縮小することのなかった区間幅が 平均値形式を用いることによって,縮小する場合があることがわかった. 一般に,f'([a,b])の値が小さいとき, すなわち,何等かの打ち消し合いが発生 し,[a,b]が小さな幅のとき,平均値形式は良い評価を与える.

以上の例は簡単なもので,めのこでも,式変形ができるが,fが 大変な数のステップからなるプログラムで実現されているとき,このような 打ち消し合いが生じたか否かは容易に判定できないので,平均値形式は 極めて有用な技法となる.しかし,一方で,このような場合,f' の計算は 容易ではない.そこで,fを計算するプログラムが与えられた場合に, f'を計算するプログラムを自動合成する手法が研究されるようになった. これを自動微分という.自動微分は数式処理による微分と違って, fを計算するプログラムを変形してf' を計算するプログラムを作るのと, 特定のxに対して,f(x)の値とf'(x) の値のみを計算すればよいので, 数式処理の微分より極めて効率的である.

2.4 自動微分

(目次に戻る)

自動微分と数式処理の記号微分の違いは,一見わかりにくいので,ここで 簡単に説明をしよう.数値計算においては,微分は差分商を用いて近似するのが 伝統的であった.差分近似では,小さなものを,小さなもので割ることにより, 桁落ちで有効数字が大幅に減少するという欠点がある.この差分近似を自動 微分で置き換えるだけでも,数値計算の信頼性が大幅に向上する.

自動微分の説明をするために,その1つの実現法を紹介しよう.そのために プログラミング技法おけるオブジェクト指向を簡単に説明する.筆者 はオブジェクト指向言語はプログラミング工学からの素晴らしいプレゼント で,ソフトウエアの開発が1000倍も楽になるという実感をもっている. 数値計算に則してオブジェクトを説明すると次のようになる.オブジェクト とはデータ構造とそのデータの演算規則のペアのことである.例えば, ベクトルというオブジェクトは,数の並びという,ベクトルデータとベクトルの スカラー倍,加減算という演算規則のペアである.また,区間という オブジェクトは区間[a,b]というデータと区間演算という演算規則のペア である.プログラムの中で

z=x+y

< というステップがあると,オブジェクト指向言語では,x,yがどういう オブジェクトか調べて,もしx,yがベクトルであれば,ベクトルの加算を してzに代入する.一方,もしx,yが 区間であれば,区間演算の加算を する.さらに,もしx,yが成分として区間をもつベクトルであれば, xyの成分を足すときは区間演算で行い, zには,成分毎に区間演算 で計算された値が,ベクトルとして代入される.

このようにオブジェクト指向言語では,一つの命令が,中身が何であるかに よって様々に解釈される.これを演算子多重定義という.したがって, ガウスの消去法のプログラムを普通の行列に対して書いておけば, 演算子多重定義の機能により,その行列が成分として区間をもって いるときには,そのプログラムは自動的に区間行列 のガウス消去法のプログラムになる.このように,オブジェクト指向言語 を用いると,プログラムの自動変換が実に容易に実現できる. これは,まるで夢のようである.

では,自動微分をオブジェクト指向で実現するにはどうしたらようので あろうか.それには微分というオブジェクトを作ればよい.微分 オブジェクトは微分型というデータとその演算規則からなる。 微分型データを(d,d') と定義する.dは関数値,d'は微分値である. 微分型の演算規則として 、まず、四則演算を

(d,d')+(e,e')=(d+e,d'+e')

(d,d')-(e,e')=(d-e,d'-e')

(d,d') * (e,e')=(de,d'e+de')

(d,d')/(e,e')=(d/e,(d'e-de')/e2)

で定義する.これには、関数の四則演算の微分則を利用している。 次に,微分型データを初等関数に代入する演算規則を, 例えばsin, cosについて定義すると

sin (d,d')=(sin d,(cos d) * d'),

cos (d,d')=(cos d,(-sin d) * d')

となる.これには合成関数の微分則であるチェインルールが利用されている。 このとき,入力xに対して関数値f(x) を計算するプログラムにおいて xを微分型データ(x,1)とすれば, そのプログラムは(f(x),f'(x)) を出力する.こうして,関数値f(x)を計算するプログラムを 関数値f(x)xにおけるfの微分値f'(x) を同時に計算するプログラムに自動変換 できることがわかる.すなわち,オブジェクト指向言語の演算子多重定義 の機能を使って,微分オブジェクトを用いれば,自動微分が実現できる.この ような自動微分の実現法を別の観点からみてボトムアップ方式と呼ぶ. トップダウン方式の自動微分,両者を巧みに組み合わせた,高速自動微分など も考えられている.[5]

こうして,オブジェクト指向言語を用いれば,関数値を計算するプログラムを 自動変換して,微分値も同時計算できるプログラムとする,自動微分が 容易に実現できることがわかる.このような自動微分の実現例からも わかるように,自動微分は数値微分のような誤差を含まないのと,手計算 や数式処理の記号微分を必要としないので,手計算によるバクが 入る心配もなく,正確,高速,かつ容易に微分値が自動微分により計算できる. これも革命的な数値計算の新技法である.

例えば,平均値形式に出て来るf'([a,b])という項は, 関数値f(x)を計算する プログラムにおいてx([a,b],1)とすれば, そのプログラムは (f([a,b]),f'([a,b]))を出力するので簡単に計算できるようになる.

3.精度保証付き数値計算と構成的数学

(目次に戻る)

3.1ニュートン法と区間解析

(目次に戻る)

精度保証付き数値計算という実用技術を生み出すもう一つのブレークスルーは ニュートン法と区間解析の結合から生まれた.ニュートン法はニュートンが プリンピキアの中で3次方程式の解法として与えたのがその起源である. 以後,極めて多数の非線形方程式の数値解法が考えられたが,結局ニュートン 法に優るものは見つかっていないといって良いであろう[6].

ニュートン法はこのような3次方程式の解法として出発したが,より一般の 非線形方程式の解法として,連立非線形方程式から非線形関数方程式まで に拡張された.今,簡単のため,f を実数の上で定義された1回微分可能な 実数値関数とする.f(x)=0の解x* を求めるために,ニュートン法では

k(x)=x-L-1f(x)

を考える.ただし,cを近似解としてLf'(c)の近似とする.L-1 が存在するとき,k(x)=xf(x)=0は同値となる. したがって,k(x)=x をみたすxを求めれば良い.このようなxk の不動点という.数値計算 法としてはx0=cとして逐次代入

xn+1=k(xn), n=0,1,2,...

を行い,nが無限大となるときx_n=x* を期待するのがニュートン法である. このような反復が収束するための十分条件は,数学では不動点定理を利用して 示される.その一番簡単で強力なのがバナッハの縮小写像原理と呼ばれる 不動点定理である.今,I=[a,b] ,c=(a+b)/2としよう. 区間[a,b] 上でkが縮小写像原理の成立条件をみたすとは

k([a,b])[a,b]の部分集合となること

および

|k'([a,b])|<1

が成立することである。このとき、kの不動点 x*が存在して、任意の 区間[a,b]に含まれるxに対して kn(x)は不動点x*に収束する。 ここで,

k(x)=c+L-1(L(x-c)-f(x))=c+L-1(L(x-c)-f(c)-f'(q)(x-c))

に注意する.ただし、pは区間[a,b]の要素で ,平均値の定理によりその存在が 保証される。cが良い近似解でf(c)が零に近く, Lf'(c)の良い 近似であればL-1(L(x-c)-f(c)-f'(q)(x-c))は零に近い.すなわち, L(x-c)f'(q)(x-c)とでかなりの打ち消し合いが 生じるはずである.そこで, k([a,b])を平均値形式で評価するというアイデアが生まれる.

K([a,b])=k(c)+k'([a,b])([a,b]-c)

がニュートン作用素の平均値形式であるが,これを最初に考えたクラフチック [7]にちなんでクラフチック作用素と呼ぶ.明らかにk([a,b])K([a,b])の部分集合となる. k(x)を計算するプログラムはf(x) を計算するプログラムから 簡単に作れる.これを利用して,k(c)が計算できる. 一方,自動微分と 区間演算により,上で説明したオブジェクト指向言語の演算子多重定義を 利用すればk(x)を計算するプログラムから自動的に k'([a,b])を計算 できる.よって,kの不動点が区間[a,b] 内に一意に存在する十分条件

K([a,b])は区間[a,b]の部分集合

|K'([a,b])|<1

が成立するか否か,kを計算するプログラムを利用して, 計算機で自動的に 判定できることになる.もし,この十分条件が満たされていれば, f(x)=0の 真の解x*,/MATH>が区間[a,b]内に 唯一つ存在することが縮小写像の原理により わかる.真の解の存在とその誤差範囲[a,b]が保証されるので, このような 数値計算による結果を精度保証付き数値解析という.

経験的に,縮小写像の原理が成立するか否かを判定する区間として

[a,b]=c+[-2|L-1f(c)|,2|L-1f(c)|]

を選ぶと良い.

3.2 構成的数学との関係

(目次に戻る)

ここで,構成的数学との関連について論じる.上記の議論によって f(x)=0の 真の解x*が区間[a,b] 内に唯一つ存在することが縮小写像の原理により わかったとする.このとき,縮小写像の原理により

[a_n,b_n]=Kn([a,b])[a,b]の共通集合

は,真の解を含み,nが無限大となるとき, anbnはともに x*に収束する ことが示される.[an,bn] を計算するプログラムはfを計算する プログラムから自動変換によって作れるから,結局,x* は計算可能実数となり,そのx* を任意精度で計算するプログラムも具体的に作れることになった. 上に述べた構成的数学では,任意精度で計算するプログラムをもって,実体 と考えることから,その立場でf(x)=0の厳密解 x*が得られたことになる. すなわち,f(x)=0の解を求めるという 連続数学の問題の厳密解が計算機を 利用して求められたことになる.マルコフ流の構成的数学の立場からは このようなことができたとしても不思議ではないという人がいるかもしれないが, 計算機を利用して積極的な側面で構成的数学を動かしている点が新しいと 筆者は考えている.なにより,``計算機では連続数学の問題の近似解は得られる としても厳密解は求まらない''という一般常識が迷信に過ぎなかったことを 示している点で,以上のような議論は画期的ではなかろうか.

4. 精度保証付き数値計算で微分方程式を解く

(目次に戻る)

確かに,連続数学の問題の厳密解が精度保証付き数値計算を利用して求まる ことはわかった.でも,一次元非線形方程式なんてたいしたことはない. 理工学の連続数学の問題としては非線形常微分方程式や非線形偏微分方程式 を解くことが重要なんだ.多分,このような印象をもたれた人が多いと思われる. しかし,カントロビッチや占部実[8],[9],[10]といった先人が既に示しているように, ニュートン法はバナッハ空間などの関数空間の上で定義できる.また, 縮小写像の原理はバナッハ空間上の非線形作用素の不動点を求めるための手法であった.彼らは,縮小写像の原理に準じる原理を利用して,バナッハ空間上の ニュートン作用素の収束定理を作った.これと,区間解析をバナッハ空間などの 関数空間に拡張したものを結合して非線形常微分方程式 や非線形偏微分方程式の境界値問題の厳密解も精度保証付き数値計算を利用 して作れるようになってきた[11],[12].すなわち,これらの問題と同値な関数空間上の ニュートン作用素の近似解が与えられたとき,その近傍で縮小写像の原理が 成立するか否か計算機で判定でき,その十分条件が満たされるとき, 解となる関数を任意精度で近似する関数を計算するプログラムが作れるように なってきた.この場合,解が正則であればこのような検証はうまくやれば 必ず成功することもわかる.ここに解x* が正則であるとは,問題をf(x)=0 とすると,解x*におけるfの微分 f'(x*)が1:1で全射となることである. また,分岐点の存在検証など正則でない解を 求めるためには,特異点の解消を行って正則化して議論することとなる. ここで,分岐点とは考えている系がパラメータを含み, そのパラメータのある値で今までとは 性質の異なった解が突然現れることをいう.x* を分岐点とするとf'(x*)は 1:1で全射とはならない.このようなx*を特異点ともいう.

4.1 ダフィング方程式の周期解の分岐

(目次に戻る)

精度保証付き数値計算による微分方程式の解析の例をあげよう[11]. 次のダフィング方程式を考える.

d2x/dt2+0.1dx/dt+x3-B cos t=0.

ただし,Bはパラメータである. ダフィング方程式は奇対称性非線形性をもつ広い範囲の系のモデルとして 有用であることが知られている。電気回路のモデルとして電気の分野でも深く 研究されてきた.特に,京大上田教授はその周期解の分岐現象を実験的に 追求し,不規則遷移現象という新概念を抽出した.これがカオスの発見 であった.このように,ダフィング方程式の分岐現象については十分に信頼で きる数値実験結果が多量に蓄積されている.しかし,そのような分岐現象の内, 数学的に証明されたものは,そのほんの一部であると思われる.これは 数値実験結果が極めて複雑で豊かな分岐現象を示唆していることによる. ここでは,最も簡単な場合であるが,このような周期解の分岐現象の存在の 証明を精度保証付き数値計算を通して行った例を示す.

ダフィング方程式においてはx(t)が周期解となるとき, y(t)=-x(t+\pi)も周期解となるという性質がある. これを対称性という. 変換S

Sx(t)=-x(t+\pi).

によって定義する.周期解x(t)Sx=x を満たすものを対称周期解, そうでないものを非対称周期解という.

ダフィング方程式の周期解に関して、Bを変化させたとき、 Bが小さい範囲では図1の概念図

に示されるような分岐現象が発生することが知られている (徳島大学川上博教授学位論文(1973))。 そこで,Bを 0から増やしながら各種周期解の枝をガレルキン法で近似的に計算したところ、 図2の結果を得た。

そこで,精度保証付き数値計算により図2に示された 対称周期解の族の存在証明を行ってみた. この近似解の枝の近くに真の対称周期解の枝が存在することと存在範囲を精度保証付き 数値計算で確定した所,図3のような結果を得た.

この計算においてはダフィング方程式を

f(x)=0

と書き換えている.ただし,

f(x)=d2x/dt2+kdx/dt+x3-B cos t

写像$f$は関数空間X={L2(0,2\pi)に含まれるx|Sx=x}の 中に定義域をもち,Xの中への写像と考えると,対称周期解は f(x)=0 の解となる.ただし,L2(0,2\pi)は区間 [0,2円周率]上の2乗可積分な関数 の作るヒルベルト空間とする.後は一次元非線形方程式のときと同様, fニュートン作用素kを考え, 縮小写像の原理の十分条件を数値的に 検証することによって,kの不動点の存在を数値的に検証し, 存在範囲を 確定している.ただし,今の場合,対称周期解の枝(一次元多様体となる)その ものの存在証明をするので,陰関数定理のニュートン法による証明に習って, 陰関数定理が成立することを図3の各点で証明し,さらに陰関数定理の成立範囲が 隣どうしで互いに重なるようにして,周期解の枝そのものの存在証明としている. この計算にはワークステーションを動かし続けて約3ヶ月要したが,時間が かかるとはいえ,周期解の枝そのものが存在証明されるとともに, それを任意精度で計算するプログラムも作られたので,周期解の枝そのものの厳密解が 得られたと考えることができる.

図3からわかるように,Bが小さいとき、 対称周期解は折り返し点を2つ.そこで、次に、 この折り返し点の存在の証明精度保証付き数値計算で行った. 折り返し点は分岐点であり,f'(x)-1 が存在しない特異点である. そこで,

d2x/dt2+kdx/dt+x3-B cos t=0,

d2y/dt2+kdy/dt+3x2y=0,

l(y)=0,

という拡大方程式を考える.ただし,x,yXに属すとし, lは規格化のための線形汎関数である. 折り返し点はこの方程式の正則解となる.この方程式の解の存在を精度保証付き 数値計算によって証明することにより,前述の対称周期解の枝が実際にB が小さい方で2つの折り返し点をもつことが証明された.

Bを大きくしていくと図1のように 対称周期解の枝から非対称周期解が分岐する. その分岐点は対称性破壊分岐点と呼ばれ,やはり適当な拡大方程式を考えることに よってその存在検証が精度保証付き数値計算によってできる. 対称周期解の 枝はBを大きくしていくと更に伸びていく.一方,非対称周期解は Bを 大きくして行くと,2回折り返し点を通過した後, 1/2分数調波解が分岐する.これも一種の 対称性破壊分岐として捕らえることができるので,精度保証付き数値計算に よってその存在証明ができる.こうして,ダフィンング方程式は異なる周期解を Bが小さい方からの適当な区間で少なくとも 1,3,1,3,7,3,7 個もつことがわかる.

4.2 ホモクリニック分岐

(目次に戻る)

前の例は局所的分岐現象の存在検証の例であった.次に,大域的分岐現象の 存在検証例として,ホモクリニック分岐の存在検証例を示そう[12]. ホモクリニック軌道とはtが-無限大で平衡点を出て,t が無限大で,その平衡点に戻ってくる軌道のことである.系に含まれるパラメタを 変化させたとき,ある値でそれまで存在していなかったホモクリニック軌道が 突然出現する現象がホモクリニック分岐である.シルニコフの定理や退化 ホモクリニック軌道の研究などから,カオスを生成する力学モデルである スメールの馬蹄力学系が埋めこまれることがホモクリニック軌道の存在から 適当な条件の下で示唆されることが知られている.この意味でホモクリニック 分岐はカオス理論などで重要となるが,その分岐現象の存在証明と存在範囲 の精細な評価は極めて困難なことが知られている.これは,数値計算的には ホモクリニック分岐点を求める問題が悪条件問題となることとして現れる.

次の非線形常微分方程式を考える:

dx/dt=y,

dy/dt=x-x2+r y+0.5xy

平衡点は(0,0)(1,0)にある. この方程式においてはr=-0.5 でスーパークリティカルホップ分岐が(1,0)で発生する.この値より パラメータrを増加させると,発生した周期解がr のある値 rcでサドル(0,0) のホモクリニック軌道になり,それより大きな rの値でホモクリニック軌道も周期軌道も消滅するというホモクリニック 分岐が起こる.ここではホモクリニック分岐点rcと ホモクリニック軌道 の数値的存在検証を行い,その値を包みこむことにする.

近似解はチェビシェフ多項式を用いる占部の方法[13]を 用いて求めた. 近似解はxについて90次のチェビシェフ多項式で計算した.この 結果rcの近似値として

rc=-969926592/2258238379=-0.429505848904...

を得た.また,ホモクリニック軌道の近似として図4のような軌道を得た.

この値及びホモクリニック軌道の精度保証を行った所, ホモクリニック分岐点の存在が証明でき,上式で与えられた ホモクリニック分岐点の誤差が高々0.0000000022で与えられることが 示された.したがって,ホモクリニック分岐点の値rc は区間[-0.429505852,-0.429505846]に含まれることがわかった.

5. むすび

(目次に戻る)

ソリトンやカオスの発見がそうだったように,計算機を用いて数値実験により これは確からしいという現象が発見される.すなわち,現代科学では 良い近似解はすでに多く発見されている.では,その近似解は本当の解であるの かゴーストか.それを見抜くのが天才であった.しかし,精度保証付き数値計算 を利用すれば,良い近似解の近くに真の解が存在することが計算機で証明できる 時代となった.すなわち,精度保証付き数値計算による厳密解の求解法は, 人間が数値実験を経て直観的に正しそうだと感じた良い近似解の存在を 前提として,厳密解を計算機で求めるという,現代科学の方法論において 不足していた部分を埋める手法であるともいえよう.

解の存在証明及び厳密解の構成という数学の論理で従来行っていた部分を, ニュートン法,縮小写像の原理,現代区間解析,構成的数学といった いわばシンプルな原理を利用し,後は精度保証付き数値計算を計算機の 能力一杯行うというような力ずくに置き換えるのである.この場合, 良い近似解をもってくるという所が人間の直観によっており,そこから このような方法によってとてつもない結果生み出されるかもしれないという 期待が生じる.このような方法論はつい最近になって実用的と認識され始めた ばかりである.このような新しい方法論によって,科学技術の新発見が次々に 成されることを期待したい.

参考文献

(目次に戻る)

[1] 山本哲朗:``エッセイ:マジソン便り'',Computrol,12,pp.5-8 (1985).
[2] 中尾充宏:``関数方程式の数値的検証法'', 数学, 42 (1990), pp.16-31.
[3] M.J.Pour-El, J.I.Richards:``Computability in Analysis and Physics'', Springer-Verlag (1989).
[4] A.Neumaier:``Interval Methods for System of Equations'', Cambridge University Press (1990).
[5] 伊理正夫,久保田光一:``高速自動微分(I),(II)'', 応用数理, Vol.1, No.2, pp.17-35, No.2, pp.53-63 (1991).
[6] 山本哲朗:``Newton法とその周辺'', 数学, 37, pp.1-15 (1985).
[7] R.Krawczyk:``Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken'', Computing, 4, (1969) pp.187-201.
[8] M.Urabe:``Galerkin's Procedure for Nonlinear Periodic Systems'', Arch. Rational Mech. Anal., 20, (1965) pp.120-152.
[9] M.Urabe:``An Existence Theorem for Multi-Point Boundary Value Problems'', Funkcialaj Ekvacioj, 9 (1966), pp.43-60.
[10] 篠原能材:``数値解析の基礎'',日新出版(1978).
[11] S.Oishi:``Numerical Verification of Existemce and Inclusion of Solutions for Nonlinear Operator Equations'', J. Computational and Applied Math., 60 (1995).
[12] S.Oishi:``Two Topics in Nonlinear System Analysis through Fixed Point Theorems'', IEICE Trans. Fundamentals, E77-A, pp.1144-1153 (1994).
[13] M.Urabe:``Numerical Solution of Multi-Point Boundary Value Problems in Chebyshev Series Theory of the Method'', Numerische Mathematik, 9 (1967), pp.341-366.

(目次に戻る)


理工ジャーナル目次へ

御意見、御質問は、oishi@mn.waseda.ac.jpまで。

本文は岩波書店より発行された雑誌``科学''1996年6月号 pp.437-445に掲載された記事をHTMLに変換したものです.