はじめに
木更津高専アドベントカレンダー19日目の記事です。何書くか悩んだんですが、いや、今(2021年12月18日16時48分)も悩んでいるんですが、まあ適当に書いていきます。
間違いがあるかもしれません。また、大学生や数学科の方には冗長かもしれませんが(というか、この記事を見られるの恥ずかしいですね)、高専って集合記号とか有界とかそういう単語とか使わないで大学工学部1～2年次の数学をするので、分かりやすいように少し補足などを加えました。
どんな方程式を扱うの？
以下の微分方程式を考えます。ただし、$x : [0,T] \rightarrow \mathbb{R}，f : [0,T] \times \mathbb{R} \rightarrow \mathbb{R}，$ ある正定数$T$と$A$に対し， $0\leq t \leq T，|x - x_0| \leq A$とします。
まあでも、微分方程式とかやったとき、有限個の点が発散してる関数とか無限遠考えたりしますよね。そこはこう、考える空間を制限したり、極限を取ったりしてください。
ほんとは式をenumerateで書きたかったのですが、なんか書けなかったので最初に宣言します...。
\begin{align}
\frac{dx}{dt} &= f(x(t),t) \ \ \ (t \in [0,T]) \\
\\
x(0) &= x_0
\end{align}
この微分方程式の解は次の条件を満たすことを仮定します．適当な$x$，$x^*$を持ってきたときにある$L$が存在して
|f(x,t) - f(x^*,t)| \leq L |x - x^*|
が成り立ちます。これはLipschitz連続とかいうやつで、これがあると嬉しいです。なぜなら、この条件を満たしてると初期値を与えたときに解が一意に定まってくれるからです。
これを満たしてないと解がたくさんあって嬉しくないです。数値計算でも重要になってきます。あとで使います。
高専ってこういう記号とか使った解説しないし、リプシッツ連続もなんだかよく分からなかったんですが(だからなんだよ。どう判定すんねんとかいう気持ち)、例えばこういうのを考えるとイメージしやすいかも。
例
適当な区間$I \subset \mathbb{R}$において、関数$f(x)$は微分可能で、その微分$f'(x)$は$I$上有界であるとします。このとき、$f(x)$はLipschitz連続です。
証明
平均値の定理を用います。$f(x)$が微分可能かつ$f'(x)$が有界なので、ある$x_1$と$x_2$が与えられたときに
f(x_1) - f(x_2) = f'(c) (x_1 - x_2)
を満たす$c \in (a,b)$が存在します。$f'(x)$は有界であることを仮定してるので、どんな$x$を選んでも勝手な正定数$M$で抑えられます。
|f(x_1) - f(x_2)| = |f'(c)| |x_1 - x_2| \leq M |x_1 -x_2|
これより、Lipschitz連続であることが示されました。
Lipschitz連続の嬉しさについて話せよ
はい。
さっきチラっといった通り、$f(x,t)$がある閉領域$\Omega = [a,b] \times [c,d]$内で連続かつ、$x$についてLipschitz連続だと、初期値を与えたときに解が一意に定まります(小泉進次郎)。嬉しいですね。
この時、次のような方法で解を構成できることが知られてます。これをピカールの近似列と言います。考える微分方程式と初期値は上の通りです。
\begin{align}
x_0 (t) &= x_0 \\
x_n(t) &= x_0 + \int_{0}^{t} f(t,x_{n-1}(t)) dt \\
x(t) &= \lim_{n \rightarrow \infty} x_n (t)
\end{align}
証明は長いので書きません。微分方程式の教科書でも見てください(例えば、笠原「微分方程式の基礎」(朝倉書店)の2章とか)。
で、Lipschitz連続でないと、この近似列の収束が保証されないことが知られています。別に$x_n (t)$を定義するだけならLipschitz連続は仮定しなくてもいいんですが、収束の保証はされません(保証されないというだけで、収束するような$f(x,t)$もあります)。
Lipschitz連続だとガッチリ一意性が保証されて嬉しいんですが、強い条件なので、条件を弱めようと頑張ってる人がいたらしいです。具体的には、岡村の条件というのもあるらしく、Lipschitz連続でない場合の一意性の条件を決めてしまうとてもすごい定理らしいですが、本記事では説明しません。Osgoodの条件というのもありますが、資料によって十分条件だったり必要十分条件だったりとか書かれてて訳が分かりません。
あとは、こんなのもあります。
Cauthy-Peanoの定理(ざっくり)
$f$が定義された領域$\Omega = [a,b] \times [c,d]$内で連続関数なら、点$(x_0,t_0) \in \Omega$の近傍で(局所)解が存在する。
ゆるすぎ。バケモンか？
点$(x_0,y_0)$近傍の解曲線は、領域$D$から出ない範囲でいくらでも延長できることも知られているらしいです。詳しくは、京大の偏微分方程式の講義ノートなどを見てください。
数値計算
この記事って数値計算したいんじゃなかったのか。
はい。そうです。
具体的にどういう手法があるのかな～っていうのを書いていきます。コードも書く気でいましたが、その辺に転がってるし、普通に書く気力がなくなってきました。数学の記事書く人ってすごいんだな。
それはそうと、書いていきましょう。さっきと添え字の意味が変わります。考えてる時刻$[0,T]$を$N$個に分割して，$t_n = n \Delta t，x_n = x(t_n)，\Delta t = T/N，n = 0,...N$としてみましょう。
また、$f$がLipschitz連続であることも仮定しちゃいます。豪華ですね。別に仮定しなくても数値計算はできると思いますが、一意性は担保されません。で、
\begin{align}
\frac{dx}{dt} = f(x(t),t) \ \ \ (t \in [0,T])
\end{align}
を解くんでしたね。やや天下りっぽいですが、両辺を$[t_{n},t_{n+1}]$で積分してみましょう(差分を取ってもいいです)。
\begin{align}
\int_{t_{n}}^{t_{n+1}} \frac{dx}{dt} dt &= \int_{t_{n}}^{t_{n+1}} f(x(t),t) dt \\
x_{n+1} &= x_{n} + \int_{t_{n}}^{t_{n+1}} f(x(t),t) dt
\end{align}
右辺の積分をどう処理するかによって、名前が変わります。
前進オイラー法
今、$\Delta t$は小さく取るので、次のような近似を考えます。
\int_{t_{n}}^{t_{n+1}} f(x(t),t) dt \approx f(x_{n},t_n) \Delta t
そして、$x_{n+1}$の更新式を
x_{n+1} = x_{n} + f(x_{n},t_n) \Delta t
とします。代入するだけで次の値が出るのでうれしいですね。こういうのを 陽解法 と言います。
精度は後進オイラー法より1次くらい劣ることが知られています。
更新則がLipschitz連続であること
$ f(x_{n},t_n)$ がLipschitz連続であることは、
|f(x,t) - f(x^*,t)| \leq L |x - x^*|
を仮定しているので、これから明らかです。
後進オイラー法
\int_{t_{n}}^{t_{n+1}} f(x(t),t) dt \approx f(x_{n+1},t_{n+1}) \Delta t
として、$x_{n+1}$の更新式を
x_{n+1} = x_{n} + f(x_{n+1},t_{n+1}) \Delta t
とします。右辺に未知の値$x_{n+1}$を含むので、この部分を数値計算しないといけないですが(こういうのを 陰解法と言います)、この手法自体は前進オイラー法より精度がいいことが知られています。いい話。
複雑な方程式じゃなかったらこの部分を手計算で$x_{n+1} = ...$の形にしてしまい、前進オイラー法っぽく計算するのもアリ。
更新則がリプシッツ連続であること
前進オイラー方程式と同じです。ただし、前に仮定したとおり、$t_{n+1}$及び$x_{n+1}$は考えてる領域内にあるとします。
Heun法
導出は略しますが(離散化誤差が$O(h^3)$となるようにうまく右辺の２つの中の$f$の値と、$f$の係数を決める感じです)
x_{n+1} = x_{n} + \frac{\Delta t}{2} \{ f(x_{n},t_{n}) + f(x_{n} + f(x_n, t_n)\Delta t,t_{n} + \Delta t) \}
とします。
更新則がLipschitz連続であること
F(x,t,h) = \frac{1}{2} \{ f(x_{n},t_{n}) + f(x_{n} + f(x_n, t_n)\Delta t,t_{n} + \Delta t) \}
とします。$|F(x,t,h) - F(x^*,t,h)|$を調べます。
|F(x,t,h) - F(x^*,t,h)| = \frac{1}{2} | f(x,t) + f(x + f(x, t) \Delta t,t + \Delta t) - f(x^*,t) - f(x^* + f(x^*, t)\Delta t,t + \Delta t) | \\
ここで三角不等式を用います。
\frac{1}{2} | f(x,t) + f(x + f(x, t) \Delta t,t + \Delta t) - f(x^*,t) - f(x^* + f(x^*, t)\Delta t,t + \Delta t) | \\
\leq \frac{1}{2} \{ | f(x,t)- f(x^*,t) | + | f(x + f(x, t) \Delta t,t + \Delta t)- f(x^* + f(x^*, t)\Delta t,t + \Delta t) | \}
Lipschitz連続を仮定してるので、この式は上から抑えられます(あってるよね...？)。
なんでLipschitz条件考えたんだよ
嬉しい定理の証明に必要だから。
収束しない方法って使っちゃいけないと思うんだ
まあそうじゃん？というわけで、次の定理を紹介します。
局所誤差の上界
$e_n$を数値解と真の解(解析解)の誤差$e_n = x_n - x(t_n)$で定めます。この時、
\max_{0 \leq t_n \leq T} |e_n| \leq \frac{C}{L} \{ e^{LT} -1 \}(\Delta t)^p
が成立します。ずっと$\Delta t$で書いてたので$\Delta t$で書いたけど、本では大体$h$で書かれてます。つまり，$t_{n} = t_0 + hn$って感じです。多分。ただし、$L$はリプシッツ定数(上から抑えるやつ)、$p$は精度の次数とします。
実際にはこんなに上手にはいかず、マシンイプシロンを考慮した$\Delta t$の決め方とかしないといけません。悲しいですね。マシンイプシロンによる誤差と丸め誤差(公式誤差)が同程度だと嬉しいです。
筑波大学の佐野先生の講義ノートに大雑把な評価の仕方は載っていますが、なんかどっかで厳密な決め方の論文か何かが転がってた気がします。
証明を書こうとしたが、力尽きた
助けてくれ。もう1日分アドカレ書かないといけないし、まあその分だと思ってそのうち追加します。
Stiff equationの話をしたかった
$\Delta t$を小さく取らないとかなり軽率に発散するテスト方程式
\frac{dx}{dt} = \lambda x
を使って色々調べてみるようなことを書きたかった。助けてくれ。
多段法の話も書こうとしたが力尽きた
助けてくれ。
お前こんな投げやりの記事で何が言いたかったんだよ
高専数学や高専でやる数値解析ってだいぶ端折った事言ってるしこういう原理の説明とか0なので、裏にこういう話があるんだよというのをしたかったが、レポートと課題の山で1日で書くのは無理だった。
参考図書
菊地文雄，斎藤宣一，数値解析の原理，岩波図書，2016
笠原晧司，微分方程式の基礎，朝倉書店，1982