高安 亮紀

スライドショーには JavaScript が必要です。

講演者：高安亮紀（早稲田大学）
題　目：解析半群を用いた半線形放物型方程式の解に対する精度保証付き数値計算法とその応用
日　時：2015年1月15日(木) 18:00-
場　所：早稲田大学西早稲田キャンパス63号館数学応用数理会議室（102号室）

オープニング

Ａ：お疲れさまです。今日は山中君がいらっしゃらないので、最初の司会だけ僕がやります。今日は3回目で、なかなかいいペースなのか、ちょっと速いのか分かんないですけど、3カ月で3回目になります。1回目のやつがそろそろWebにかたちが出来上がると思いますので、また見てください。
　今日は3回目で、早稲田大学の高安先生による『解析半群を用いた半線形放物型方程式の解に対する精度保証付き数値計算法とその応用』ということです。よろしくお願いします。

講演

高安：ご紹介ありがとうございます。早稲田大学の高安亮紀です。本日はお足元の悪い中お越しいただきまして、誠にありがとうございます。
　第3回ということで、少しだけ数学的な方も良いかなと思いまして、ゴリゴリに数式が今日はいっぱい出てきますのでご容赦ください。
　解析半群を使うということで、これまでの精度保証付き数値計算法とはちょっと違った見方から、数値計算法を確立していこうというのが今回やりたいことです。その方法をご紹介します。

　自己紹介のスライドを作ったんですが、皆さんは顔見知りですので（笑）、今回は飛ばします。ただ、一応「僕、専門分野は数値解析です」と言いたいんですが、特に偏微分方程式とカッコを付けることでこだわりを見せています。

　今日やりたいことは、計算機を用いて非線形の偏微分方程式を解析することです。この解析という意味は数値計算を使って数学的に厳密に議論を展開したいという意味合いを込めて、解析アプローチと銘を打っています。
　物語はとてもきれいで、そこに向かう道筋もどうやらあるっぽいんだけど、「その中に入ってったら大変なことになったよ」という感じのお話を今日はしたいと思います。僕はよくこのスライドを研究に見立てて使うんですが、すぐ内容が曇ってしまって分からなくなるんですが、今日はゆっくりゆっくりと一個ずつお話しさせていただければと思います。

　今日のお話なんですが、共同研究者が3名いらっしゃいます。一人目は早稲田大学の大石研究室の学生さんであるD2になりました水口さんです。後ろにいらっしゃってます。あと、筑波大学の久保先生、そして早稲田大学の大石先生です。この3名とのディスカッションがないとこの話は成り立たなかったので、今日は僕一人で今喋ってますが、この3名との共著論文となっております。

　問題です。半線形放物型方程式を以下に書き下します。 $\Omega$ は2次元平面上の有界多角形領域とします。その領域上で定義された( $P_J$ )という方程式を考えるんですが、いま $u$ が時間と空間、 $t$ と $x$ に対する関数だと思ってください。 $\partial_t$ は $u$ の $t$ 方向の偏微分を表します。 $A$ はあとで明示しますが、空間方向の微分作用素です。非線形項を $f(u)$ と書き、境界条件と初期値をセットしてあげると、放物型の偏微分方程式として初期値境界値問題となります。ここで、 $J$ を時間区間とします。本研究で考える $J$ は $(t_0,t_1]$ というかたちをしているときもありますし、 $J$ が $0$ から $\infty$ 、無限大まで考えるときもあります。本講演ではこの2パターンを考えます。

非線形項 $f$ はあまりにも一般的にすると解析ができなくなってしまうので、2階のフレッシェ（Fréchet）微分がとれる非線形項と本講演では仮定します。それに対して、 $u = 0$ on $\partial\Omega$ がディリクレ（Dirichlet）境界条件を表し、初期関数 $u_0$ はこちらで与え、 $H^1_0$ の滑らかさを持つような初期関数とします。
　今後、 $L^p$ は $p$ 次のルベーグ（Lebesgue）可積分が可能なバナッハ（Banach）空間を表します。 $H^1$ は $L^2$ 空間に対する1階のソボレフ（Sobolev）空間と、 $H^1_0$ は $H^1$ のうちディリクレ境界条件が課されている関数空間にです。すなわち、今考えている初期値境界値問題の空間方向は常に $H^1_0$ と思ってください。以上の問題設定の中で考えていくのですが、微分作用素 $A$ だけちょっと書き切れなかったので次のスライドで書きます。

　 $A$ はあるドメイン( $D(A)$ )、このドメインは領域によって変わります、このドメインから $L^2$ への線形作用素とし、 $A$ はこのようなかたちをしている空間方向の微分作用素とします。ラプラシアン（Laplacian）を若干一般的にしたものと思ってください。ラプラシアンは2階偏微分作用素ですが、 $A$ は1階微分の所に $a_{ij}$ という行列関数を掛けてあげて微分する微分作用素になっています。これが $a_{ij} = a_{ji}$ をみたす、つまり対称な行列関数で、関数がそれぞれ対称であることが今回の条件です。

各 $a_{ij}$ は $W^{1,\infty}$ 、つまり $L^\infty$ の1階のソボレフ空間の元とします。これは1階微分しても有界ということを示しています。 $a_{ij}$ がこの条件を満たしているとするとし、さらに $\xi$ を $\mathbb{R}^2$ の任意のベクトルとして次の不等式をみたす $\mu>0$ があるとすると、 $A$ という作用素が必ず正の方向に持ち上がっている微分作用素となります。この状況下で初期値境界値問題 $(P_J)$ を考えます。

　 $L^2$ 空間の内積とか $H^1_0$ のノルム（norm）とか、 $H^1_0$ の共役空間を $H^{-1}$ と書きますよとか、厳密にお話をするべきですが、ここでは記号の定義だけですませます。特に何も工夫はしていないノルムと内積をここで導入します。

　我々がやりたいことは、精度保証付き数値計算をやりたい。特にこの講演では、初期値境界値問題 $(P_J)$ の弱解を考えます。時間区間 $J$ の中で弱解を $u(t):=u(t,\cdot)$ と書きます。 $u(t)$ が $H^1_0$ を満たしていて、かつ、初期関数が $u_0$ とする。 $t = t_0$ とすると、 $u(t_0)=u_0$ になるという意味です。

さて弱解は元の方程式の解の領域次第で、古典的な解に滑らかさを回復することができたりしますが、ここでは詳細は話しません。代わりにこの弱解の存在を数値的に検証しようというのがここでやりたいことです。

$a$ は双線形形式です。 $a:H^1_0\times H^1_0\to\mathbb{R}$ という内積のようなものです。さらに $a(u,v)$ は、さきほど定義された $a_{ij}$ を使って書けます。さっき $a_{ij}$ が正の関数で定義されているので、 $a(u,v)$ は内積になります。それを記述する方法として、 $a(u,u)$ が下から $\mu$ ...さきほどの $\mu>0$ と、この $\mu$ は一緒です...と抑えられているとします。 $a_{ij}$ は1階微分しても上から抑えられますので、上に有界です。
ということで、この弱解の存在と一意性を計算機を利用して証明するというのが、我々がやりたいことです。

$X$ を $J\times\Omega$ 上のあるバナッハ空間とします。この $X$ に対して数値解 $\omega$ というのを自分で作ります。そして $X$ 上で中心が数値解 $\omega$ 、半径 $\rho$ の閉球の中に数学的に厳密に包み込むことを考えます。これが精度保証付き数値計算で時間発展方程式を抑え込むための一番最初のアイデアです。

　先行研究は、日本の中尾先生、木下先生、木村先生による2012年の結果があります。2012年、2013年、これは2013年と書いてありますが2014年出版の論文がそれぞれ1本ずつ、目標は一緒です。目標はスライド9ページにあるような、解をあるBanach空間の閉球の中に包み込むというのをやりたい。先行研究では、バナッハ空間 $X$ が時間方向に $L^2$ 関数で空間方向に $H^1_0$ 関数とそれぞれ弱い関数空間でとります。弱い関数空間ということは条件が弱いですから精度保証がしやすいです。

一方、我々は、空間方向に関しては $H^1_0$ の条件はそのままに、時間方向は $L^\infty$ で考えます。ほとんど至るところ有界になっているような時間方向で解を追っていったらどうだろうと。少しだけ先行研究に比べて条件を強めてあげて解の検証を行っていくかたちです。

　先行研究の3本に対して、実は力学系の観点から見た先行研究があります。これだけではないんですが、代表的なのを2つ抜き出してきました。

こちらは領域を矩形領域に限定してしまって、フーリエ級数展開します。そうすると、無限個の項を持つ非線形方程式が出てくるので、精度保証付き数値計算のフレームワークで解くことをしています。ただ未知係数が無限個あるので、無限をどうするかというところに新規性があります。

もう一つは70ページぐらいある論文で、空間方向の微分可能性を解析してあげれば、シャープな誤差がえられるでしょうと。方程式がdissipative PDEを考えていて、空間方向に4階微分を課します。すると、ものすごい勢いで解が拡散していくんです。このような解に関してうまくいった手法を記述したら70ページになっちゃったということですね。すごくタイトに誤差が抑えられていて、数学者が本気で問題を解くならこういう感じになるのかな。とてもタイトな評価が得られています。

　さて我々はちょっと力学系からは外れて、日本の藤田宏先生の論文で勉強をはじめました。その中で紹介されているのは、離散半群（半群を離散化したもの）を使った数値スキームの解析結果で、これらがまとまった本が2001年に出ています。この本で我々はたくさん勉強をしました。
じゃあ離散半群じゃなく、もとの半群を用いて解の検証をしてあげればうまくいくんじゃないか、というのがこの研究のスタートです。

　まずは時間区間 $J=(t_0,t_1]$ を1つの区間に限定します。すると、 $(P_J)$ は以下のようになります。これに対して近似解 $\omega$ 、解の指針となる中心の解を自分で作ります。 $V_h$ を $H^1_0$ の有限次元部分空間とします、つまりコンピュータで構成できるような部分集合を考えます。

その $V_h$ から2点、任意に2つ持ってきます。これはどう持ってきてもいいですが、 $\omega$ がよい近似解ならば精度保証付き数値計算も成功しやすいです。この2点を使って、 $\omega(t)$ を $\omega(t)=\hat u_0\phi_0(t)+\hat u_1\phi_1(t)$ と構成します。この $\phi_0(t)$ 、 $\phi_1(t)$ は、 $t_0$ 、 $t_1$ における線形のラグランジュ基底（Lagrange basis）です。

この基底2つを使ってあげて2点を張ってあげると、 $\hat u_0$ と $\hat u_1$ 次第で近似解の精度が決まります。

　どのように解を包み込んでいくかという話をするとき、まず近似解と初期関数との誤差が $\varepsilon_0$ で抑えられるという条件が必要です。 $u_0$ は与えられた初期関数ですから、この誤差は簡単に計算することができます。よって、この $\varepsilon_0$ を基に $L^\infty(J;H^1_0(\Omega))$ という関数空間の中で解を捉えます。これは下のノルムでバナッハ空間になります。すなわち、近似解 $\omega$ を中心に半径 $\rho$ の閉球の中で解を包み込む。これはさっき説明した目的と全く同じことです。

　 $A$ は空間方向の微分作用素でした。作用素 $A$ の弱形式を $\mathcal{A}$ とします。そうすると、 $-\mathcal{A}$ は $e^{-t\mathcal{A}}$ という半群を生成します。この半群は、各時刻において $H^{-1}$ から $H^{-1}$ への作用素となります。ただし時刻が変わるので、時刻が変われば作用素も変わる。すなわち作用素の族となります。この $\mathcal{A}$ の性質が良いため、 $\mathcal{A}$ から生成された半群（ $e^{-t\mathcal{A}}$ ）は解析半群となります。

ここで一番重要なのは、次のようなかたちのCauchy問題に対して、これは $\mathcal{A}$ をラプラシアンに変えれば熱方程式、初期値を $u_0$ にするとこの半群を使って解を一意に書くことができます。つまり、基本解のような役割を半群はしているんだということです。この事実は以下からわかります。詳しくは各種の教科書を参照ください。

　この半群を用いて、今日は次のような定理をご紹介します。これが、我々が今回導き出した定理。はじめに初期関数と近似解との誤差が $\varepsilon_0$ で抑えられているとします。次に半群を用いて次のような積分評価、実はこの積分は時刻 $t$ による関数で、時間方向に $L^\infty$ とった評価が $\delta>0$ で抑えられるとします。ただし積分の中は、 $\omega$ は近似解ですから、この部分は元の方程式に $\omega$ を入れたような残差を表しています。すなわち、残差に半群を掛けて、 $t_0$ から $t$ まで積分をしてあげたものの $L^\infty$ 評価です。これが $\delta$ で抑えられる。

さらに、この非線形項 $f$ に対して、任意の考えたいボールの中から $\varphi$ と $\psi$ を持ってきてあげて、引き算してあげると上から定数 $L_{\rho_0}$ で抑えられるという条件を課します。local Lipschitz conditionといわれるやつです。この $L_{\rho_0}$ があるならば、次の不等式をみたすかどうかチェックします。 $\mu$ とか $M$ は問題を決めれば一意に決まる値で、 $\varepsilon_0$ は最初の誤差です。 $\tau$ は時間ステップの幅です。 $L_{\rho}$ はここの $L_{\rho_0}$ で、 $\delta$ はここで計算可能な $\delta$ です。これら全て計算可能になりますから、左辺よりも右辺の方が大きいということを逐一チェックします。これを満たすような $\rho$ を一つ決定します。

すると、元の初期値境界値問題 $P_J$ の解 $u(t)$ が、我々が考えた、中心が数値解そして半径が $\rho$ のボールの中に一意存在するということがいえます。

　証明の概要を述べます。
$z$ を各 $t$ に対して $H^1_0(\Omega)$ の元とし、 $z(t)=u(t)-\omega(t)$ というふうに書けるとおきます。そうすると、実は任意の部位に対して $z$ が満たす式は、右辺がこのように書けます。ここを $g(z)$ と置いてあげると、 $S(z)$ という非線形写像、 $L^\infty(J;H^1_0(\Omega))$ から $L^\infty(J;H^1_0(\Omega))$ への非線形写像が定義できます。

　この非線形写像が $Z$ という関数空間の中で不動点を持つかどうかをチェックします。
この $Z$ は $L^\infty(J;H^1_0(\Omega))$ のノルムの意味で半径 $\rho$ で、中心 $0$ の閉球です。
この $Z$ の中に必ず $z = S(z)$ を満たす不動点が存在することになり、すると $z(t)=u(t)-\omega(t)$ より解が一意に存在するということがいえます。
実は不動点 $z$ が存在することと、解が存在することが同値であることからどちらか一方を示すことで存在証明は完了します。

　 $S(Z)\subset Z$ をチェックします。ごめんなさい、この $\zeta$ は $u_0$ です。そして上からこの値で抑えられます。

　さっき定義した $g(z)$ をリプシッツ（Lipschitz）連続で抑えられそうな項と残差の項で分けます。
それに対してリプシッツ連続の項で抑えていくんですが、特筆すべきは、この半群の評価ができるということが一番今回重要な点の一つです。今まで困っていた部分が計算可能な定数で抑えられることになりました。

　 $g_1(s)$ の部分は上からリプシッツ連続の項だけが残るというようなかたちになって、ここら辺が半群の名残なんですけど、半群の名残は計算可能な定数になりました。

　リプシッツ連続の項はこれで抑えられるし、もう一個の残差の項はそのまま $\delta$ で抑えられていますから、 $S(z)$ のノルムは上から $\rho$ で抑えられるならば、 $S(Z)\subset Z$ が成立します。

　さてバナッハの不動点定理を使いたいので、次に一意性の条件をちょっとチェックしないといけません。
$z_1$ 、 $z_2$ を $Z$ から持ってきて評価してあげると、こんな式が出てきます。実は、これはさっきの条件が成り立つならば $1$ 未満であるということが自動でいえるので、 $S$ は縮小写像であるということが証明できます。

　ただし、いつでもあるというわけではなくて、大前提としてこれらがちゃんと数値的に確認できないといけない。従ってこいつは十分条件になっています。この十分条件が成り立てば解が一意存在するということになります。だから、近似解の精度が悪いと残差が大きくなってしまって、「こんな $\rho$ はない」と言われることもありますし、問題によってはこのリプシッツ項が大きいので、「解がありません」と失敗するときもあります。

　次に初期関数から時刻を進めていくにあたって、次の時刻、次の時刻、次の時刻……とこの定理を何回も適用したいわけです。そうなると、 $\varepsilon_0$ として前回の誤差評価を使います。前回の誤差評価がここに現れるので、いわゆるこれはラッピングエフェクト（Wrapping effect）といいます。前回の誤差が伝播するわけです。従って過剰に誤差評価をしていくことになり、結果検証が失敗する。なかなかそこまでうまくいくものでもない。万能ではないということです。

　さて一回、解を捉えることができれば、次のような誤差評価が成り立ちます。先ほどは $t_0$ から $t_1$ までの $J$ という時間区間の中で $L^\infty$ でした。しかし、こちらは時刻を $t=t_1$ に固定した $u(t_1)-\hat u_1$ の誤差評価が計算できます。先ほどのやつにちょっと $e^{-\tau\lambda_{\min}}$ が付いたやつと、 $\delta$ が $\tilde\delta$ になります。

実は $\rho$ というのは時間の中で一番最悪ケースをとってくるんですけど、時刻 $t=t_1$ に縮めてあげると誤差が縮むんです。少し分かりづらいんですが。
よって、 $t_0$ から $t_1$ まで最悪ケースを考えて誤差評価し、 $t=t_1$ における出口をちょっと考えて、これは水口君の言葉なんですが「出口調査」をしたいと。すると、そのときの誤差評価が小さくできるということになります。

　 $1$ 個の時間区間で評価できましたから、次はそれを何回も何回も繰り返してあげる。すると、先ほどの出口調査の誤差評価 $\varepsilon_1$ が $\varepsilon_0$ というふうに置き換わって、もう一度主定理を適用してあげてチェックする。またチェックする、チェックする……と繰り返します。

　そうすると、どんなふうに近似解を持ってきてもいいんですが、ひとまずこの後の例で使っているので、Backward Eulerを使いましょう。Backward Eulerの例は、このようなもともとここが時刻 $u$ の時間微分でしたけれども、 $u$ の時間微分を $1$ 階のいわゆる差分法で近似する。そしてこんな感じで近似的に次の時間ステップを解く。これで各時刻における近似解が出てきますので、それらを線形でつないであげれば、それを $\omega(t)$ として計算していく。

　ということで、ちょっと文章で言うよりも絵で描いた方が分かりやすいので、こんな図を描きました。

横軸が $t$ で、縦が今考えたい関数の位相です。だいたい時刻に $L^\infty$ 、空間に $H^1_0(\Omega)$ という感じです。 $u_0$ という点があって、これは与えます。これに対して、 $\hat u_1$ 、 $\hat u_2$ 、 $\hat u_k$ は近似解ですのでどうやってやってもいいですけど、ひとまずBackward Eulerで計算したものとします。
それらを線形でつないで、最初の区間、 $0$ から $t_1$ までの間でまずは $\rho$ を計算する。すると、この中に必ず解が一意に存在する。こんな感じのボックスの中に必ず解があるということを保証します。さらに端点での出口調査をしてあげて、それからまたもう一度同じ操作をやる、さらにもう一度、このように繰り返していくことで、時刻 $t_k$ までの誤差評価を各ステップごとに評価します。

　ここまでのまとめです。実はここまでで30分の予定だったんですけど、だいぶ大幅に遅れてます。 $J_k$ という時間区間において、問題 $(P_J)$ を考えると、その解の存在と一意性を逐次的に数値解の近傍に包み込むことができるようになりました。逐次的にやるとどうしても前回の誤差が伝播します。伝播する誤差をどうにか抑え込まないといけないというのが、今後の課題になってきます。

一応、数式で書くと、元の初期値境界値問題の解は各 $k=1,2,\cdots,n$ について、実は $B(\omega)$ というチューブのまわりに一意存在することが計算機援用証明できるということになっています。

　ちょっと例をいくつか示してみようと思います。僕は藤田宏先生の本から勉強し始めたので藤田型方程式が大好きなんです。この藤田型方程式、べき乗の非線形項を持つような方程式に対して解がどんな振る舞いをするのかをチェックします。
$u_0$ はこのようなかたちの問題を考えます。実は藤田型方程式というのは初期関数が大きければ解が爆発するということが既に知られています。Large data blow-upといって、この $\gamma$ を大きくしていけば、いずれ解が爆発する点が存在します。ただし $\gamma$ が小さい場合は時間大域的に存在することも分かっているので、その $\gamma$ のギリギリまで攻めていきたい。そうゆう問題設定です。
今、 $h$ はメッシュサイズ、空間方向に有限要素法の $2$ 次の基底、時間方向にはBackward Euler methodを使って近似解を構成して解を追っていきます。

　すると、これが結果になっているんですが、0から0.0039062という小さな区間で「まずは $\varepsilon_k$ がこうなって、 $\rho_k$ がこんなでした」というのを並べています。ちょっと分かりづらいので図にします。

　 $\gamma= 1$ のときに初期誤差はここです。横軸が $t$ で、縦軸が各 $\rho$ です。すなわち、 $t$ が $0$ から $t_1$ のときに誤差がこれだけ、 $t_1$ から $t_2$ のときに誤差がこれだけと、定数関数で誤差をプロットしています。そうすると、ちょっと時間がたつと小さくなっていくということが分かります。

$\gamma= 30$ とか $50$ と上げていくと、 $30$ のときは頑張って戻ってきてくれるんですが、 $\gamma= 50$ のときは頑張りきれずにそのまま誤差が爆発していきます。これは誤差評価が爆発してしまうだけで実際の近似解は実はちゃんと収束しているので、blow-upを示しているわけではない。ですが、割と大きな $\gamma$ に対して計算ができるようにしていくことでこの精度保証の方法を改善するいい指標ができています。

　非線形項がべき乗だけだとちょっと寂しいので、 $u-u^3$ みたいにちょっぴり変えてみます。これはアレン＝カーン（Allen-Cahn）方程式の一種ですが、これに対して初期値がこれで計算してみます。パラメータ $h$ は空間方向のメッシュサイズ、 $\tau$ は時間方向の時間区間幅で、これらを変えて計算します。

　はじめに $h$ をまあまあ小さくとってあげて、 $\tau$ を変えてます。そうすると、 $\tau$ が大きいと誤差は大きいのですが、 $\tau$ が小さくなればなるほど誤差が小さいのかと言ったらそうではない。赤い所が実はこのラインなんですけど、赤い所と緑の所が拮抗している。いい勝負してます。ですが、 $\tau$ を小さくとり過ぎると全体的に誤差は大きくなってしまうんです。

これは $\tau$ を小さくとれば精度保証する回数が単純に増えます。よって、前の誤差が伝播し、蓄積してしまうんです。従って、なるべく $\tau$ は大きめに。でも大きくとり過ぎるとダメなんだけど、「大きくとり過ぎない程度に」とると良い精度で計算できます。数値解析の授業とかでEuler法を学ぶ時に「 $\tau$ は小さくとり過ぎない方がいいんだよ。でも、大きくとり過ぎてもダメなんだよ」と聞くと思うのですが、それと同じ状況です。

　次に $\tau$ を十分小さくとってあげて $h$ を小さくとってあげると、きれいに落ちていきます。もちろん $\tau$ の値によってで打ち止まりが来ますので $\tau$ 次第ということになりますが、 $\tau$ と $h$ のバランスは設定が難しいです。

　ここまで来ると、手法としては何回も何回もステップを繰り返して検証を繰り返すことで、時間方向にも弱解の存在を追っていくことができるようになりました。ここで僕はだいぶ満足したのですが、大石先生が「これは時間無限大にはいかないのかい？」とおっしゃられまして、「いきません」とは言えませんでした。「こんなに誤差が収束しているんだから、きっと誤差が収束することが証明できるはずだ」という大石先生のご指摘から、解析をしてみました。それが次のパートです。

　ある程度時刻がたってしまえば、その後はもう解は定常解に落ち着くだろうというのが、感覚として我々がスタートした部分です。今回は解のうち $H^1_0(\Omega)$ に属し、 $t$ が $(0,\infty)$ で存在する解を時間大域解と呼びます。この時間大域解の存在が計算機援用証明できないかなと考えます。

$t'$ をある時刻とします。この時刻は精度保証できるところまでというイメージです。そうすると、 $t'$ から $\infty$ までの解が存在することがいえればいい。存在する範囲を示すには、どのような関数空間に存在するかが重要になります。 $t'$ までなら数値計算した結果がありましたので、数値計算した結果のまわりに解が存在するということがいえました。だけど、 $t'$ から $\infty$ までの間は数値計算できないので、指針がありません。じゃあ、どうしようか。僕らは定常解のまわりに存在することをいえばいいとなりました。定常解はいわゆる楕円型方程式の解。これはこれまでの精度保証付き数値計算のテクニックを使えば、楕円型方程式の解は精度保証付き数値計算ができます。そこで得た定常解のまわりに $t'$ から $\infty$ で近づいていくような時間大域解が存在ことを数値検証します。

まとめると、ある時刻 $t'$ までは先ほどのスキーム、Concatenationで連続でつなげていき、その後定常解まわりの解の存在を示します。結論として、 $L^\infty(0,\infty)$ の空間で時間大域的に存在する解をこの関数空間で抑えます。

　先行研究を調べたら、実は早稲田の木村拓馬先生からS. Caiさんが2012年に学位論文でこういうことをしていると教えていただきました。そこでは反応拡散系を考え、あるクラスの定常解に関する精度保証付き数値計算法を示す学位論文でした。その中の一部に $(t',\infty)$ で定常解まわりに大域的に存在する範囲を解析半群を用いれば計算できると言ってます。ここで考えた関数空間は、 $L^\infty(\Omega)\times L^\infty(\Omega)$ 。
要素が $2$ つある方程式なので、空間方向に $2$ つの直積で考えたバナッハ空間上で半群を生成しないといけない。学位論文では確かに定常解を抑えることができ、この定常解のまわりに解があることもいえるというフレームワークを示しています。
しかし、ある時刻 $t'$ までのつなぎ方の記述はありませんでした。単に時刻 $t'$ から $\infty$ までの定常解まわりの時間大域存在が数値検証できるという話をしています。我々の話も同じような話です。実はこれは我々の話ができてから後で教えられて、「おお、一緒のようなことをしてますね」という状態なんです。

　では問題を考えます。今度は $(0,\infty)$ で考え、大域存在する解を考える。ここでは $A =-\Delta$ で固定します。

　ここから先は、ある時刻 $t'$ において $\|\eta-\hat u_n\|_{H_0^1}\le\varepsilon_n$ が成り立っていると仮定します。その仮定の下で問題 $(P_G)$ 考えます。これは $(t',\infty)$ で定義された初期値境界値問題です。
それに対して、 $\phi$ というのをこの方程式の定常解とします。そうすると、この $\phi$ のまわりに解が一意存在することが示せるようになることがここでの目標です。

　また記号を定義します。 $B$ と書いたらボールです。時刻は $(t',\infty)$ 。ある時刻から先の全部を考える。その時刻においてフレッシェ微分が上からこんな値で抑えられるというのが、今回新たに加える仮定になります。ただし、 $f$ が $2$ 階のフレッシェ微分可能ならばこの値は存在するので、条件としてはあまり変わっていません。

　バナッハ空間 $X$ というのを $X_\lambda$ とします。 $X_\lambda$ は何かというと、先ほどは $L^\infty((t',\infty);H^1_0(\Omega))$ というようなバナッハ空間だったんですけど、そのノルムに $e^{(t-t')\lambda}$ という重みを付けてあげた関数空間です。これは解析手法でよく使う重みです。この $X_\lambda$ は、このノルムの元でバナッハ空間になります。

　ということで定理のかたちで紹介すると、まずは解が時刻 $t'$ まで包含されていることが第一条件。それに対して定常解 $\phi$ もちゃんと近似解近傍で存在していることが第二の条件。
近似解のまわりに必ず $\phi$ があるということは、精度保証付き数値計算を行えば分かります。
さらに $\lambda$ を $1$ 個、これは自由度が出てきてしまうんですが、0から $\lambda_{\min}/2$ の間で1個値を決めてあげる。そして次のような不等式が $1$ 個成り立てば、 $U_\phi$ という部分集合……これはボールです。 $\phi$ を中心として半径 $\rho$ のチューブの中で必ず解が一意存在します。 $t'$ から無限大まで必ずこの $\rho$ から出ないということがいえます。

$\rho$ から出ないということが一個重要な部分、実はこの $X_\lambda$ はexponentialの重みを付けてますから、exponentialの重みを外してあげると実は $\phi$ に向かって収束する。実はこの定理はこのようなこともいえます。ただ、今の興味はちゃんと時間大域存在するかどうかということなので、時間大域存在がいえる条件が、この不等式となり、これを計算機でチェックしてあげればOKということになるわけです。

　さて、また証明です。でもほぼ一緒です。 $f(z)$ というのをこのかたちで定義してあげれば、上から抑えていき、バナッハの不動点定理を使って計算するというかたちです。

　一点だけ注意しないといけないのが、こちらの項。定常解が $\phi$ なので、その $\phi$ と初期関数 $\eta$ の誤差評価なんてできるのかなと。ここは近似解を使ってどうにか評価を計算したいという話にします。つまり出口調査した結果と $\hat u_n-\hat\phi$ の項、そして定常解の誤差評価にわけて評価してあげる。そうすると、出口調査の部分の評価は $\varepsilon_n$ で抑えられるし、定常解の誤差評価は $\rho'$ という値で抑えられる。あとはもう両者とも近似解ですから、単に引き算して評価してあげれば評価はできる。よって、どうやらConcatenation schemeと、もうひと条件チェックしてあげることで、 $t'$ から $\infty$ まで全部チェックできるとなりました。

　では数値例を出しますが、藤田型方程式で頑張ろうと思います。これは水口君の計算結果を紹介します。関数空間 $V_h$ は $(0,1)^2$ の正方領域で考えてますから、フーリエ基底で張って考えます。水口君すごくて、Backward EulerじゃなくてCrank-Nicolsonで近似解を作ってくれたので、そいつを使って計算しました。

　 $\gamma$ の値は初期関数に与えたパラメータ。それを大きくしてあげます。すると、 $\lambda= 0.01$ のときはConcatenation schemeの反復 $n = 5$ 解でで時刻 $t' =0.046875$ から先は必ず時間大域存在しました。
$\gamma$ を大きくします。すると時間大域存在を示せる時刻 $t'$ が先になります、その $t'$ が $\gamma= 5$ ぐらいまで大きくすると、 $n = 47$ ぐらい反復して、時刻 $0.375$ で $\rho= 0.17656$ という値で評価でました。

よって、 $X_\lambda$ の関数空間でとっているので、 $u(t)$ のノルムは $\rho e^{-(t-t')/40}$ の値で、 $t'$ から $\infty$ で $\phi=0$ に向かって収束するということになります。

実はこの藤田型方程式はこの収束先が零関数なんです。これは解析結果で、ラージデータを与えてあげるとblow-upして、スモールデータを与えると0にいくというのがあります。あとおそらく $1$ 個だけ神のみぞ知る軌跡があって、実は不安定定常解にいく不安定多様体の方向から近づけてあげるといく軌跡があると思うのだけど、通常の数値計算ではできないので、通常は爆発していくか、0に収束していくのどちらかになります。

　次に定常解が非零であるものを考えます。 $u^2$ はそのままに、そこに関数を与えてあげます。すると、定常解が非零な関数に近づきます。計算すると定常解の近似解 $\hat\phi$ はこんな形をしています。この状況では前の例よりもフーリエモードがたくさん必要になりますが、 $10$ 個のフーリエ級数を用いて計算すると $\rho=0.04035$ で、そのとき $t'$ は $0.2578125$ でした。このとき定常解のまわりに時間大域解の存在が証明できました。

　まとめます。一応、解析半群を使うことが今回の手法で新しいところです。解析半群を用いる精度保証付き数値計算方法を今回紹介しました。これはある一つの定理を何回も何回も繰り返しチェックすることによって実現します。Concatenation schemeという名前を今回付けました。さらに、精度保証付き数値計算を用いた時間大域解の存在証明もできるじゃないかということで、定常解のまわりに包み込む方法を紹介しました。

今後の課題は、今、半線形の放物型方程式をやっていますので、それをいろいろな反応拡散方程式に適応していきたい。特にS. Caiさんがやっているような連立している放物型の方程式系、反応拡散系とよばれているものをやっていきたいと思います。

さらに、波動方程式も解析半群が生成するということが分かっていますので、波動方程式の解析半群を用いた精度保証付き数値計算もおそらくできます。波動方程式ができればシュレディンガー（Schrödinger）方程式もできるって大きなことを言ってたら、大石先生が渋い顔をされていました。「あまり大きなことは言うんじゃない」という話ですね。

もう少し夢を見ると、無限次元力学系と関連がとてもあるので、力学系の観点から見た、ただ単に不動点があるということではなくて、「不動点にどのような近づき方をするのか」という振る舞いを記述できれば、もっと精度保証付き数値計算に幅が出るんじゃないかなと思い、無限次元力学系とも関連していきたいと思っています。

近々では藤田型方程式の爆発時刻、いったいいつ爆発するのか。「今でしょ」じゃないんですけど、爆発時刻の精度保証付き数値計算法を確立したいと思います。

長々と話しましたが、これでお終いです。どうもありがとうございました。（拍手）

質疑

Ａ：質問等あればお願いします。
Ｂ：定常解のときって、広くとることもあえてできますよね。
高安：評価式が成り立てばいいので、右辺はなるべく小さくしたい。小さくなくてもいいんだけど、この部分が大き過ぎると、失敗します。よって $\eta-\phi$ の値がだいぶ小さい値になるように作りたい感じです。
Ｂ：定常解のところの一意存在がいえてる範囲にとれないんですか？
高安： $U_\phi$ を作るところがポイントです。 $U_\phi$ を作ると、問題によって評価が決まるんです。そして、時刻 $t = t'$ になったときに、この問題の解の包含が $U_\phi$ の中に入っているかどうかをチェックしている。つまり、 $U_\phi$ を作成した時点で、時刻 $t = t'$ の時点のチェックだけでいいんです。
Ａ：さっきの九大の人がやっている的が小っちゃいっていうのは、左側のとこがすごく小っちゃいっていう……
高安：ここです。はい。
Ｂ：あっちが小っちゃいのか。
Ａ：そうです。確かにそれは厳しいかも。
高安：この大きさに対してこの $\varepsilon_n$ が入ってればOKなんですけど。現状、なかなか難しいですね。
Ｂ：そこ、あえて広げることはできますよね。
高安：この $\lambda$ はパラメータなので、ここを操作すると広げることはできます。つまり、 $\lambda$ を $\lambda_{\min}/2$ まで近づけていくと、この $\rho$ は大きな値になります。すると、不等式は満たしやすくなるかもしれないし、 $\rho$ がどんどん大きくなるから評価が難しくなるかもしれない。どちらともいえません。
Ｂ：たぶん今はできる限り小さく小さくというのをやろうとしてると思うんですけど、実は $\rho$ をあえて1に近づくようにとっていくとたぶんそこの幅が広がるので、なんかとりやすくなるような気がしたんです。
高安：1に近づくように……
Ｂ：近づくように、かつ、速く無条件で。
高安：リプシッツの定数が1に近づくようにすると大きくなるのかな。
Ｂ： $\rho$ から……もうちょっと考えてみる。
高安： $\rho$ は、時間大域解を示すにあたっては小さければいいというもんでもない。なるべく大きくとってもいい。
Ｂ：一意存在ならば、なお大きければいいという話ですよね。
高安：そうね、そうね。
Ｂ：あえて大きくという……。
高安：小さくとろうという誤差評価の感覚とはまたちょっと違いますね。
Ｂ：一意で分けるので。
Ａ：最初の定義の一意存在っていうのは、ボールの中で一意って言ってるんですか？　全体で一意？
高安：一応、この一つのボールの中で一意。こっちのボールの中でも一意なので、ここが大きくなっていけば、絶対にこの中に一意。
Ａ：ごめんなさい。その外にはあるかもしれないっていうこと？
高安：そうです。この外には軌道があるかもしれない。
Ａ：だからうまく解が得られるんだ。一意存在の範囲は $\rho$ をあえて大きくとれば、もっと一意存在の幅が広げれるっていう……
Ｂ：そうもいえるよね。
Ａ：だから最小の $\rho$ と最大の $\rho$ を2個作ると、解はこの狭い範囲にあって、一意の幅はもっとでかい範囲っていうのを作れる。
Ｂ：たぶんいえるんでしょうね。なので、その広い方がっていう。
Ａ：なるほど。
高安：そういう意味で広い方？
Ｂ：はい。
高安：ああ、なるほど。
Ａ：一意の幅が広ければ広いほどいいから。
高安：そういう意味では、広い方はまだとってない。小さい方しか使ってないです。
Ｂ：高安さんのたぶん小さいってのは存在の方で、一意の幅はメチャメチャ広くあります。残差は小さいし。
高安：そうなのかな。分かりました。もう少し調べる必要があります。ありがとうございます。
Ａ：時間の刻み幅を小っちゃくしていったときに誤差のバウンド $\rho$ が大きくなっちゃうっていうのは、精度保証を繰り返しやるからっていう感じだと思うんですよね。近似解としてはいいものになってるはずなんですか？
高安：近似解としては……
Ａ：そうとも限らない？
高安：いや、いいものになっているはず。
Ａ：じゃあ、近似解を作るときの刻み幅と精度保証するときの刻み幅って、厳密に一緒にする必要はないですよね。
高安：ないです、ないです。
Ａ：だから、近似解は持ってるけど、あえて飛ばし飛ばしでやった方が狭めれる可能性はある？
高安：そっちの方がいいかもしれません。
Ａ：なんかうまく組み合わせて、どの刻み幅でやったらいいかっていろんなので試して、最小のやつをとっていくみたいなことができるよね。
高安：そうですね。「この近似解、いいな」っていうのが分かるのが、この残差 $\delta$ の値ですね。残差の評価が良くなるならば「これはいい近似解だな」と。
Ａ： $\rho$ をチェックするコストっていうのはそんなに大きくないんですか？
高安： $\rho$ の不等式ですか？　これはもう微々たるものです。
Ａ：じゃあ、本当に毎近似解でやっていくのと飛ばしでやるのをいろいろ同時にやって、最小のやつを探しにいくみたいなことでも、あまりコストは掛からない？
高安：でも、残差の計算だけ少し時間が掛かるので、そこだけが大変……。どんな近似解がいいかがここで分かるので。
実は同じようなことをやってみたことがあります。例えばBackwardですごい時間幅を小さくして10回ぐらいいった先の近似解とはじめのやつを使って検証してみたら、そんなに残差は良くならないんです。残差の評価の仕方がうまくできてない部分があるのかも。この評価部分を少し変えると良くなるのかもしれないですけど、本質的には変わらないかも。
Ａ：残差って、両端じゃなくて間を線形につないだ残差になるから？
高安：そうですね。
Ａ：だからキツイのか。そういうこと。線形を線形じゃなくすることはできるの？（笑）それは難しいのか。
高安：今は一応、時間方向に $1$ 次近似です。例えば常微分方程式を考えるときは時間方向に $24$ 次とかとるので、圧倒的に精度が足りませんね。ただ、今はそれに耐えうる近似解を用意できないのと、もし $24$ 次の関数を持ってきたとしても、残差を計算する評価式が今はない。頑張ればできるかな。時間方向に高精度な解を得る方法は今後の課題の一つです。
Ｃ：僕が聞いてなかっただけだと思うんですけど、数値実験してるじゃないですか。あのときの方程式ってどういうかたちなんですか？
高安：どの数値実験で？
Ｃ：なんかピョンピョン、ピョンピョン……（笑）大きさが小っちゃいの。
高安：ああ、ピョンピョン、ピョンピョンやってた時。これは最後の数値実験で、このときに $k = 1$ 、 $l = 1$ のやつです。言ってなかった、ごめんなさい。これは1から3まで全部足し合わせたやつをこの実験ではやっているんだけど、 $k = 1$ 、 $l = 1$ の場合です。
Ｃ： $4$ は付いてる？
高安：うん、 $4$ は付いてる。
Ｃ：っていうことは、定常解はゼロ解に対応するか。安定解？
高安：そうです。
Ｃ：すみません、あともう一個なんですけど。だいたい誤差が最初上がって、下がってったじゃないですか。相対誤差はやっぱり上がっていくのが観察できるんですよね。
高安：上がっていく？
Ｃ：要するに、相対誤差が上がっていくのと解のノルムが下がっていくのを掛け合わせると、なんかターニングポイントというかがあって、そこから下がっていくみたいな？
高安：そのとおりです。解のノルムが小さくなるので。例えば解が大きくなってた時は誤差評価は大きくなっていたと思います。
Ｃ：相対誤差ってだいたいどのぐらいのオーダーで上がっていくと、線形とか2乗ぐらいのとか、だいたいグラフとか使って……。
高安：いや、グラフは作ってないんですけど、logスケールで直線で上がっていくような感じ。
Ｃ：ああ、やっぱり基本、指数オーダー？
高安：うん。そこは解像化が必要です。
Ａ：他はないですか。大丈夫ですか。じゃあ、ありがとうございました。
高安：ありがとうございました。（拍手）

プロフィール

早稲田大学理工学術院総合研究所　次席研究員

1984年千葉県出身。2008年、早稲田大学教育学部理学科数学専攻卒業。10年、同大大学院基幹理工学研究科修士課程修了（理学）。12年、同大学大学院基幹理工学研究科博士課程修了（理学）。主な研究テーマは微分方程式の解の精度保証付き数値計算。日本応用数理学会、日本数学会、日本シミュレーション学会、各会員。