応用数学

同義語:applied mathematics

概要

応用数学(applied mathematics)とは、数学の外から来た問題を扱う数学の総称で、現象を方程式で表すモデル化、方程式の解を調べる解析、解を数値で求めて誤差を保証する計算の三つの作業とその理論からなる。代表的な結果として、Newton の冷却則の解と観測からの係数の決定、Lipschitz 条件のもとで Euler 法の誤差が刻み幅 $h$ に比例する評価、刻み幅が大きいと減衰する解の近似値が発散する安定性の条件がある。微分方程式の解法、積分変換、数値解析、最小二乗法、次元解析などの分野を横切って、一つの問題を最後まで解くときの手順の地図として使う。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: 常微分方程式, Taylorの定理, 行列の階数

熱いコーヒーを机に置くと、温度は初めは速く、やがてゆっくり室温へ近づく。10 分後に何度になったかを測れば、30 分後の温度や、飲みごろになる時刻を予測できるだろうか。この問いに答えるには、まず現象を方程式に書き直し、方程式から結論を引き出し、必要なら計算機で数値を求め、最後にその数値が現実と合うかを確かめる。応用数学はこの一続きの作業を支える数学である。
この往復は次の図のようにまとめられる。上の段は現実の側、下の段は数学の側であり、矢印の一つ一つに数学の道具が対応する。
$$ \begin{CD} \text{現象} @>{\text{モデル化}}>> \text{方程式}\\ @A{\text{検証}}AA @VV{\text{解析・計算}}V\\ \text{予測} @<{\text{解釈}}<< \text{解} \end{CD} $$
つまり、現象を方程式に移すのがモデル化、方程式から解を得るのが解析と計算、解を現象の言葉に戻すのが解釈、予測をデータと比べるのが検証である。どの段にも、条件を外すと結論が崩れる箇所がある。この記事では、応用数学の範囲を定め、分野の見取り図と代表的な公式の表を示したうえで、冷却の問題と数値解法の誤差を例に、公式が正しいことと手順が正しい答えを出すことを証明する。

定義と範囲

応用数学

応用数学とは、数学の外から来た問題を扱うための数学の総称であり、次の三つの作業とその理論からなる。
(i) モデル化:現象を変数・パラメータ・方程式で表す。
(ii) 解析:方程式の解の存在・一意性・性質を調べ、閉じた形の解や近似式を求める。
(iii) 計算:解を数値で求める手順を作り、その誤差と停止の条件を評価する。

ここで言う「数学の外」は物理学・工学に限らない。化学反応の速さ、薬の血中濃度、生物の個体数、通信路の雑音、金融の価格、画像の復元など、量の関係を式で書ける問題はすべて対象になる。逆に、扱う問題が数学の内部から来ていても、上の三つの作業の形をとる研究は応用数学と呼ばれることが多い。応用数学は扱う道具では決まらず、問題への向き合い方で決まる分野だと言ってよい。
理論の数学との違いは、主張の正しさの基準ではない。応用数学の公式も、仮定のもとで証明されて初めて使える。違いは、仮定そのものを現象から選ぶ点にある。例えば「温度の変化の速さは室温との差に比例する」という仮定は数学からは出てこない。仮定を選んだ時点で、結論はその仮定の範囲でだけ意味をもつ。したがって応用数学では、定理の仮定が現実のどの状況に当たるかを確かめる作業が、証明と同じ重さをもつ。
もう一つの特徴は、答えが数値で求められることが多い点である。「解は存在する」だけでは、30 分後の温度は分からない。解を閉じた形で書けない方程式では、数値解法で近似値を求める。そのとき近似値と真の値の差、つまり誤差がどれだけかを保証できなければ、数値は予測として使えない。誤差の評価は応用数学の中心的な定理の一つの型である。

分野の見取り図

応用数学の主な分野を、典型的な問いと代表的な道具で並べる。右の列は、それぞれの主題を詳しく扱う記事である。

分野典型的な問い代表的な道具記事
数理モデル現象をどんな方程式で表すか保存則、比例の仮定、コンパートメント数理モデル
次元解析方程式を立てる前に関係の形を絞れるか次元の行列、Buckingham の $\Pi$ 定理次元解析
微分方程式の解法時間・空間の変化を解く求積法、Laplace 変換、変数分離法Laplace変換による微分方程式の解法, 変数分離法(偏微分方程式)
積分変換微分を掛け算に変えるLaplace 変換、Fourier 級数・変換Laplace変換, Fourier級数の計算
複素関数の応用実の積分・級数の値を求める留数、Jordan の補題、主値留数による定積分の計算
数値解析近似値と誤差の保証Euler 法、Newton 法、数値積分誤差(数値解析), Newton法
線形代数の応用データに模型を当てはめる最小二乗法、特異値分解正射影と射影行列
確率・統計ばらつきのある量を扱う確率分布、推定、乱数による模擬乱数とシミュレーション
最適化最もよい選び方を求める勾配法、Lagrange の乗数、線形計画線形計画法(高校数学)
漸近解析・摂動法小さい・大きいパラメータで近似する漸近展開、Laplace の方法Stirlingの公式

この表の分野は互いに独立ではない。冷却の問題を例にとると、方程式を立てる段で数理モデル、パラメータを決める段で最小二乗法、方程式を解く段で微分方程式の解法、解を数値で求める段で数値解析を使う。一つの問題を最後まで解くと、いくつもの分野を横切ることになる。表は道具の置き場所の地図として使い、問題の流れは後の「計算の手順」で述べる。

代表的な公式と使える条件

応用数学の公式は、使える条件と組で覚える必要がある。次の表は、この記事と関連する記事で扱う代表的な公式を、条件と根拠とともに並べたものである。根拠の列の label はこの記事の定理を、記事名はその記事で証明された定理を指す。

手法公式・結論使える条件根拠
Newton の冷却則$T(t)=T_a+(T_0-T_a)e^{-kt}$$k>0$、室温 $T_a$ が一定prop-am-cooling
冷却の速さの決定$k=\dfrac1{t_1}\log\dfrac{T_0-T_a}{T_1-T_a}$$T_1-T_a$ と $T_0-T_a$ が同符号で $\lvert T_1-T_a\rvert<\lvert T_0-T_a\rvert$prop-am-cooling
Euler 法の誤差$\lvert y(t_n)-y_n\rvert\le\dfrac{Mh}{2L}\bigl(e^{L(t_n-t_0)}-1\bigr)$$f$ が $y$ について Lipschitz 定数 $L>0$、$\lvert y''\rvert\le M$thm-am-euler-error
Euler 法の安定$y'=-\lambda y$ で $y_n\to0$ $\iff$ $0< h<2/\lambda$$\lambda>0$prop-am-euler-stability
Newton 法誤差が 2 乗で減る単根、$C^2$ 級、初期値が根に近いNewton法
最小二乗法$A^{\mathsf T}A\,x=A^{\mathsf T}b$$A$ の列が 1 次独立正射影と射影行列
Laplace 変換の微分$\mathcal L[f'](s)=sF(s)-f(0+)$$f$ が絶対連続、$f,f'$ が指数オーダーLaplace変換
Buckingham の $\Pi$ 定理$\Pi_0=F(\Pi_1,\dots,\Pi_{n-r})$量が正の値、関係が単位の取り替えで不変次元解析

表の「使える条件」の列を落とすと、公式は誤った答えを出す。その実例は後の「使える条件と反例」で数値とともに示す。最小二乗法の公式は 正射影と射影行列 の定理「最小二乗解と正規方程式」で、Newton 法の速さは Newton法 の定理「局所収束と2次収束」で証明されている。この記事では、それらを使う側の手順に集中する。

モデルから解へ:冷却の問題

冷却の問題を最初から最後までたどる。物体の温度を $T(t)$、室温を一定の $T_a$ とする。Newton の冷却則は「温度の変化の速さは、室温との差に比例する」という仮定であり、式では
$$ T'(t)=-k\bigl(T(t)-T_a\bigr),\qquad k>0 $$
と書ける。右辺の負号は、物体が室温より熱ければ冷え、冷たければ温まることを表す。$k$ は物体の形・材質・周りの空気の流れで決まる定数であり、理論から値は出てこない。そこで観測から $k$ を決める。

冷却則の解と速さの決定

$k>0$、$T_a,T_0$ を実数とする。初期値問題 $T'=-k(T-T_a)$、$T(0)=T_0$ の解は $\mathbb R$ 上でただ一つで、
$$ T(t)=T_a+(T_0-T_a)e^{-kt} $$
である。さらに $T_0\ne T_a$ とし、時刻 $t_1>0$ での値 $T_1$ が、$T_1-T_a$ と $T_0-T_a$ が同符号で $\lvert T_1-T_a\rvert<\lvert T_0-T_a\rvert$ を満たすなら、観測 $T(t_1)=T_1$ と両立する $k>0$ はただ一つで、
$$ k=\frac1{t_1}\log\frac{T_0-T_a}{T_1-T_a} $$
である。

$u(t)=e^{kt}\bigl(T(t)-T_a\bigr)$ とおく。積の微分により
$$ u'(t)=ke^{kt}(T-T_a)+e^{kt}T'=ke^{kt}(T-T_a)-ke^{kt}(T-T_a)=0 $$
なので、$T$ が方程式を満たすことと $u$ が定数であることは同値である。$u(0)=T_0-T_a$ だから、解は $T(t)=T_a+(T_0-T_a)e^{-kt}$ に限られ、逆にこの関数は方程式と初期条件を満たす。これで存在と一意性が示された。
次に $k$ を決める。観測 $T(t_1)=T_1$ は $(T_0-T_a)e^{-kt_1}=T_1-T_a$ と同値である。仮定から比 $\rho=(T_1-T_a)/(T_0-T_a)$ は $0<\rho<1$ を満たすので、$e^{-kt_1}=\rho$ は $k=-\log\rho/t_1$ とただ一つの解をもち、$\log\rho<0$ からこの $k$ は正である。$\rho$ の逆数の対数をとれば主張の式になる。$\square$

証明の前半は 1階線形微分方程式と積分因子(高校数学) の定理「1 階線形微分方程式の解」の特別な場合でもある。ここで条件の役割を見ておく。$\rho<0$ なら、モデルの解は室温の反対側へ移ることがないので、観測とモデルが矛盾している。$\rho=0$ なら、モデルの解は $T_0\ne T_a$ のとき有限の時刻で室温に達することがないので、やはり矛盾している。$\rho\ge1$ なら、差が縮んでいないので $k>0$ と矛盾する。観測が式の仮定と両立するかを先に確かめるのは、モデルを使う手順の一部である。

コーヒーの冷却

室温 $20$ ℃ で、$90$ ℃ のコーヒーが $10$ 分後に $60$ ℃ になった。30 分後の温度と、30 ℃ まで冷める時刻を求める。
使う公式:prop-am-cooling。
答え:$k=\frac1{10}\log\frac74\approx0.0559616$(毎分)、$T(30)=20+70\left(\frac47\right)^3=\frac{1620}{49}\approx33.0612$ ℃、30 ℃ になる時刻は $t=\dfrac{10\log7}{\log(7/4)}\approx34.7723$ 分。

途中の計算を開く

手順 1:$\rho=(60-20)/(90-20)=4/7$ なので $k=\log(7/4)/10$。手順 2:$e^{-30k}=\rho^3=64/343$ から $T(30)=20+70\cdot64/343=20+640/49$。手順 3:$70e^{-kt}=10$ を解いて $t=\log7/k$。

答えの 33 ℃ は、モデルが正しい限りの予測である。実際のカップでは、表面からの蒸発や、カップ自体が熱を吸う効果があり、比例の仮定は近似にすぎない。後の例題で、複数の観測からの当てはめと、予測の検証を扱う。

解の数値計算と誤差の保証

冷却の方程式は閉じた形で解けたが、多くの方程式はそうはいかない。そのときは数値解法で近似値を求める。最も単純な数値解法が Euler 法である。初期値問題
$$ y'=f(t,y),\qquad y(t_0)=y_0 $$
に対し、刻み幅 $h>0$ を決めて $t_n=t_0+nh$ とし、
$$ y_{n+1}=y_n+h\,f(t_n,y_n) $$
で近似値 $y_n\approx y(t_n)$ を順に求める。これは各時刻で接線に沿って $h$ だけ進む手順である。
近似値を予測として使うには、誤差 $e_n=y(t_n)-y_n$ の大きさを保証しなければならない。次の定理は、刻み幅を半分にすれば誤差の上界も半分になることを示す。数値解法の「次数 1」とはこの意味である。

Euler 法の大域誤差

$f$ を $[t_0,t_0+T]\times\mathbb R$ 上の連続関数とし、ある $L>0$ に対して
$$ \lvert f(t,y)-f(t,z)\rvert\le L\lvert y-z\rvert\qquad(t\in[t_0,t_0+T],\ y,z\in\mathbb R) $$
が成り立つとする。初期値問題 $y'=f(t,y)$、$y(t_0)=y_0$ の解 $y$ が $[t_0,t_0+T]$ 上で $C^2$ 級で、$\lvert y''(t)\rvert\le M$ を満たすとする。Euler 法の近似値を $y_n$($t_n=t_0+nh\le t_0+T$)とすると、
$$ \lvert y(t_n)-y_n\rvert\le\frac{Mh}{2L}\Bigl(e^{L(t_n-t_0)}-1\Bigr). $$

誤差を $e_n=y(t_n)-y_n$ とおく。$e_0=0$ である。
1 歩でのずれ。 真の解を Euler 法の式に入れたときの残り
$$ \tau_n=y(t_{n+1})-y(t_n)-h\,f\bigl(t_n,y(t_n)\bigr) $$
を評価する。$y'(t_n)=f(t_n,y(t_n))$ なので、$\tau_n=y(t_n+h)-y(t_n)-hy'(t_n)$ である。Taylorの定理 の定理「Lagrange剰余つきのTaylorの定理」を 1 次で使うと、$t_n$ と $t_{n+1}$ の間のある $\xi_n$ で $\tau_n=\frac{h^2}2y''(\xi_n)$ となる。よって $\lvert\tau_n\rvert\le\frac{Mh^2}2$。
誤差の漸化式。 真の解の式 $y(t_{n+1})=y(t_n)+hf(t_n,y(t_n))+\tau_n$ から Euler 法の式 $y_{n+1}=y_n+hf(t_n,y_n)$ を引くと
$$ e_{n+1}=e_n+h\bigl(f(t_n,y(t_n))-f(t_n,y_n)\bigr)+\tau_n. $$
Lipschitz 条件により第 2 項の絶対値は $hL\lvert e_n\rvert$ 以下なので
$$ \lvert e_{n+1}\rvert\le(1+hL)\lvert e_n\rvert+\frac{Mh^2}2 . $$
漸化式を解く。 $\lvert e_n\rvert\le\frac{Mh^2}2\sum_{j=0}^{n-1}(1+hL)^j$ を $n$ についての帰納法で示す。$n=0$ では右辺は空の和で $0$ である。$n$ で成り立てば、上の不等式から
$$ \lvert e_{n+1}\rvert\le\frac{Mh^2}2\sum_{j=1}^{n}(1+hL)^j+\frac{Mh^2}2=\frac{Mh^2}2\sum_{j=0}^{n}(1+hL)^j $$
となり $n+1$ でも成り立つ。等比数列の和を計算すると
$$ \lvert e_n\rvert\le\frac{Mh^2}2\cdot\frac{(1+hL)^n-1}{hL}=\frac{Mh}{2L}\bigl((1+hL)^n-1\bigr). $$
最後に、すべての実数 $x$ で $1+x\le e^x$ が成り立つ($e^x-1-x$ は $x=0$ で最小値 $0$ をとる)ので、$(1+hL)^n\le e^{nhL}=e^{L(t_n-t_0)}$。これで主張を得る。$\square$

この定理の上界は $h$ に比例するので、$h\to0$ で近似値は真の値に収束する。ただし係数 $e^{L(t_n-t_0)}$ は時間とともに指数的に大きくなりうる。長い時間を計算するときに誤差の保証が急に弱くなるのはこのためである。一方、上界は最悪の場合を見積もったものであり、実際の誤差はもっと小さいことが多い。
誤差が小さいことと、近似値が真の解のふるまいを再現することは別の問題である。次の命題は、刻み幅が大きすぎると、減衰するはずの解に対して近似値が発散することを示す。

Euler 法の安定性の条件

$\lambda>0$ とし、$y'=-\lambda y$、$y(0)=y_0\ne0$ に Euler 法を刻み幅 $h>0$ で使う。このとき $y_n=(1-h\lambda)^ny_0$ であり、
(i) $0< h<2/\lambda$ なら $y_n\to0$、
(ii) $h=2/\lambda$ なら $\lvert y_n\rvert=\lvert y_0\rvert$ が続き、
(iii) $h>2/\lambda$ なら $\lvert y_n\rvert\to\infty$
となる。真の解 $y_0e^{-\lambda t}$ はどの場合も $0$ に収束する。

$f(t,y)=-\lambda y$ なので Euler 法は $y_{n+1}=y_n-h\lambda y_n=(1-h\lambda)y_n$ であり、帰納的に $y_n=(1-h\lambda)^ny_0$。$q=1-h\lambda$ とおくと、$\lvert y_n\rvert=\lvert q\rvert^n\lvert y_0\rvert$ である。$\lvert q\rvert<1$ は $-1<1-h\lambda<1$、つまり $0< h\lambda<2$ と同値であり、このとき $\lvert q\rvert^n\to0$。$h\lambda=2$ なら $q=-1$ で $\lvert y_n\rvert=\lvert y_0\rvert$。$h\lambda>2$ なら $q<-1$ で $\lvert q\rvert^n\to\infty$。$y_0\ne0$ なので結論を得る。$\square$

この命題は thm-am-euler-error と矛盾しない。$h$ を固定して $n\to\infty$ とする極限と、時刻 $t_n$ を固定して $h\to0$ とする極限は別のものだからである。$\lambda$ が大きい方程式、つまり速く減衰する成分を含む方程式では、精度の要求より厳しい $h<2/\lambda$ という条件が刻み幅を決める。このような方程式は硬い方程式と呼ばれ、陰的な解法を使う理由になる。

計算の手順

応用の問題を解く手順を、冷却の問題に沿って番号付きで述べる。各段で使う公式と、手順が止まる条件を添える。

  1. 量と仮定を書き出す。 変数(温度 $T$、時刻 $t$)、パラメータ($k$、$T_a$)、仮定(室温が一定、差に比例)を列挙する。仮定ごとに、それが崩れる状況を一つ書いておく。
  2. 次元をそろえ、可能なら無次元化する。 式の両辺の次元が一致するかを確かめる。$T'=-k(T-T_a)$ では $k$ の次元は時間の逆数である。次元解析 の定理「Buckingham の Π 定理」で、パラメータの数を減らせることがある。
  3. 方程式を解析する。 解の存在と一意性を確かめる(常微分方程式 の定理「Picard–Lindelöf の定理」)。閉じた形で解けるなら解く(prop-am-cooling)。平衡点と、その安定性を調べる。
  4. パラメータを決める。 観測が 1 つなら prop-am-cooling の式、複数なら最小二乗法(正射影と射影行列 の定理「最小二乗解と正規方程式」)を使う。列が 1 次独立でないなら、パラメータはデータから決まらないので、観測を増やす。
  5. 数値で解く。 閉じた形で解けないときは数値解法を使う。刻み幅は prop-am-euler-stability の安定の条件を満たすように選び、thm-am-euler-error の上界か、刻み幅を半分にした計算との差で誤差を見積もる。
  6. 検証し、解釈する。 モデルの予測を、パラメータを決めるのに使っていないデータと比べる。合わなければ手順 1 に戻って仮定を直す。
    手順 5 の「刻み幅を半分にした計算との差」について補足する。次数 1 の解法では誤差がおよそ $Ch$ なので、$h$ と $h/2$ の計算の差はおよそ $Ch/2$ であり、細かい方の誤差とほぼ同じ大きさになる。この見積もりは $h$ が十分小さいという仮定に頼っているので、差が刻み幅とともに半分ずつ減っているかを、3 段以上の刻み幅で確かめるのが安全である。
    手順 6 で使うデータを、手順 4 で使ったデータと分けておくことは重要である。パラメータの数を増やせば、手元のデータにはいくらでもよく当てはまる。しかしその模型が新しい状況を予測できるかは、当てはめに使わなかったデータでしか分からない。応用数学の結論が「この仮定のもとで、この精度で、この範囲の予測ができる」という形をとるのはこのためである。
    場合やること結果の形
    方程式が線形で係数が定数積分因子・Laplace 変換で解く閉じた形の解
    平衡点のまわりの小さな変化線形化する指数関数の和で近似
    閉じた形で解けない数値解法、刻み幅の条件を確かめる数値の表と誤差の見積もり
    パラメータがデータから決まらない観測を増やす、モデルを簡単にする決まるパラメータだけの模型

例題

Euler 法の誤差の減り方

$y'=y$、$y(0)=1$ を Euler 法で $t=1$ まで解き、真の値 $e$ と比べる。
使う公式:thm-am-euler-error($L=1$、$M=\max_{[0,1]}\lvert y''\rvert=e$)。
答え:$y_n=(1+h)^{1/h}$。$h=0.1$ で $y_{10}\approx2.593742$、誤差 $\approx0.124539$(上界 $\approx0.233539$)。$h=0.05$ で誤差 $\approx0.064984$、$h=0.025$ で誤差 $\approx0.033218$。刻み幅を半分にするたびに誤差はほぼ半分になる。

途中の計算を開く

Euler 法は $y_{n+1}=(1+h)y_n$ なので $y_n=(1+h)^n$、$t=1$ では $n=1/h$。上界は $\frac{eh}2(e-1)$ で、$h=0.1$ では $0.233539$。誤差の比は $0.124539/0.064984\approx1.916$、$0.064984/0.033218\approx1.956$ で、$2$ に近づく。

図: Euler 法の誤差と刻み幅 図: Euler 法の誤差と刻み幅
図は、$y'=y$ を $t=1$ まで解いたときの誤差と、定理の上界を、両対数の目盛で描いたものである。どちらも傾き $1$ の直線に乗り、誤差が $h$ に比例することが読み取れる。上界は実際の誤差のおよそ 2 倍で、最悪の場合の見積もりとしては十分に鋭い。

硬い方程式と刻み幅

$y'=-50y$、$y(0)=1$ を Euler 法で $t=1$ まで解く。
使う公式:prop-am-euler-stability($\lambda=50$、条件は $h<0.04$)。
答え:$h=0.1$ では増幅率 $1-5=-4$ で $y_{10}=(-4)^{10}=1048576$。$h=0.01$ では増幅率 $0.5$ で $y_{100}=2^{-100}\approx7.8886\times10^{-31}$。真の値は $e^{-50}\approx1.9287\times10^{-22}$。

途中の計算を開く

$y_n=(1-50h)^n$。$h=0.1$ では $n=10$、$h=0.01$ では $n=100$。$h=0.01$ の値は真の値と桁は大きく違うが、どちらも $0$ にきわめて近く、減衰するというふるまいは再現されている。

冷却のデータへの当てはめ

室温 $20$ ℃ で、時刻 $t=0,10,20,30$(分)に温度 $90,61,44,34$(℃)を観測した。$\log(T-20)=a-kt$ の形に最小二乗法で当てはめ、$t=40$ の温度を予測する。
使う公式:prop-am-cooling の解の形の対数、正射影と射影行列 の定理「最小二乗解と正規方程式」。
答え:$a\approx4.249369$、$k\approx0.0536383$。予測は $T(40)=20+e^{a-40k}\approx28.1976$ ℃。

途中の計算を開く

手順 1:$Y_i=\log(T_i-20)=\log70,\log41,\log24,\log14$。手順 2:$A$ の行を $(1,-t_i)$ とし、正規方程式 $A^{\mathsf T}A(a,k)^{\mathsf T}=A^{\mathsf T}Y$ を解く。$A^{\mathsf T}A=\begin{pmatrix}4&-60\\-60&1400\end{pmatrix}$ は正則。手順 3:$e^a\approx70.0612$ は初期の差 $70$ にほぼ一致する。

この当てはめでは、温度そのものではなく対数をとった値の二乗誤差を最小にしている。観測の誤差が温度に一様に乗るなら、対数をとると室温に近いデータの誤差が拡大されるので、厳密には重みを付けるか、温度のまま非線形の最小二乗法を使う。どちらを選ぶかはデータの誤差のモデルによる。応用数学ではこのように、計算の手順の選び方そのものが仮定に依存する。
ex-am-cooling では観測 1 つから $k\approx0.0560$ を得たが、ここでは 4 つの観測から $k\approx0.0536$ を得た。二つの値の差は、観測に含まれる誤差と、比例の仮定のずれの両方を反映している。どちらの $k$ を使うかで 40 分後の予測は約 $0.73$ ℃ 変わる(観測 1 つの $k$ では $T(40)=20+70\left(\frac47\right)^4\approx27.4636$ ℃)。予測にこの程度の幅があることを添えて答えるのが、応用の問題の答え方である。

使える条件と反例

外す条件反例成り立たなくなること
$f$ が $y$ について $\mathbb R$ 全体で Lipschitz$y'=y^2$、$y(0)=1$ を $[0,1.5]$ で解が $t=1$ で発散して区間全体で存在しない。Euler 法は有限の値を返し、誤差の定理は使えない
刻み幅の安定の条件$y'=-50y$ を $h=0.1$ で減衰する解の近似値が $(-4)^n$ で発散する
列が 1 次独立列がすべて等しい $A$ への当てはめ最小二乗解が一つに定まらない
$T_1-T_a$ と $T_0-T_a$ が同符号冷却で $T_1-T_a$ が $T_0-T_a$ と逆符号両立する $k>0$ がない
観測の差が縮む$T_a=20$、$T_0=90$、$T_1=95$$\rho=75/70>1$ で $k<0$ になり、$k>0$ と両立しない

1 行目の反例を確かめる。$y'=y^2$、$y(0)=1$ の解は $y=1/(1-t)$ であり、$t\to1$ で無限大になる。$f(t,y)=y^2$ は $\mathbb R$ 全体では Lipschitz 条件を満たさないので、thm-am-euler-error の仮定が成り立たない。それでも Euler 法の漸化式 $y_{n+1}=y_n+hy_n^2$ は計算でき、$h=0.1$ で $t=1.5$ まで進めると約 $1.62505\times10^4$ という有限の数を返す。数値解法は解が存在しないことを知らせてくれない。$t=0.9$ では真の値は $10$ だが、Euler 法の値は $h=0.1,0.01,0.001$ で $4.289186$、$8.278593$、$9.777491$ と、ゆっくりにしか近づかない。爆発の近くでは $M=\max\lvert y''\rvert$ が大きく、上界が役に立たないからである。
2 行目は ex-am-stiff で確かめた。方程式の右辺は Lipschitz 定数 $L=50$ をもち、誤差の定理の仮定はすべて満たされている。それでも $h=0.1$ では近似値が百万を超える。定理の上界 $\frac{Mh}{2L}(e^{50}-1)$ は天文学的に大きく、この刻み幅では何も保証していない。仮定が満たされていても、結論の上界が使える大きさかどうかは別に確かめる必要がある。
3 行目では、例えば $A=\begin{pmatrix}1&1\\1&1\\1&1\end{pmatrix}$、$b=(1,2,3)^{\mathsf T}$ とすると、$x=(x_1,x_2)$ は和 $x_1+x_2$ だけを通して $Ax$ に入る。残差は $x_1+x_2=2$ のすべての点で同じ最小値 $\sqrt2$ をとり、$(2,0)$ と $(0,2)$ はどちらも最小二乗解である。モデルの言葉では、二つのパラメータが同じ効果しかもたず、データから区別できないということである。これは 正射影と射影行列 の例「反例:列が 1 次独立でないと最小二乗解は一つに定まらない」と同じ現象である。
4 行目は、モデルと観測が矛盾する場合である。$T_a=20$、$T_0=90$ で $T_1=15$ を観測したなら、温度が室温を越えて下がったことになり、Newton の冷却則のどの $k>0$ でも説明できない。この反例では $\lvert T_1-T_a\rvert=5<70=\lvert T_0-T_a\rvert$ で差は縮んでおり、破れているのは同符号の条件だけである。この場合は計算を進めず、室温が一定という仮定や、測定そのものを疑うべきである。
5 行目は、同符号のまま差が広がる場合である。$T_1=95$ なら $\rho=75/70>1$ で、$e^{-kt_1}=\rho$ の解は $k=-\log\rho/t_1<0$ となり、冷却を表す $k>0$ とは両立しない。

よくある誤り

誤り正しくはなぜ
数値解法が値を返したので解は存在する存在は解析で確かめる1 行目の反例
誤差の上界が小さいので長い時間も安心上界は $e^{L(t-t_0)}$ で増えるthm-am-euler-error
精度が足りれば刻み幅は任意安定の条件 $h<2/\lambda$ も必要prop-am-euler-stability
パラメータを全部データで決められる列の 1 次独立性を確かめる3 行目の反例
当てはめに使ったデータで検証する別のデータで検証する当てはまりの良さは予測の良さではない

補足

応用数学の各分野には、それぞれ独立した理論がある。微分方程式の解の存在と一意性、Fourier 級数の収束、数値解法の収束と安定性の関係、最小二乗法の統計的な意味などである。数値解法については、1 歩でのずれが小さいことと、解法が安定であることの二つがそろうと収束が得られる、という一般的な原理が知られている。Euler 法の定理はその最も簡単な場合であり、Lipschitz 条件が安定性の役割を果たしている。
モデル化の段については 数理モデル で、保存則からの方程式の立て方、平衡点の安定性、パラメータの決め方を扱う。方程式を立てる前に関係の形を絞る方法は 次元解析 で扱う。現実の問題から方程式を立て、解析し、計算し、解釈するという往復の考え方は、微分方程式の入門書でも最初に述べられることが多い(例えば LebDQ26 §0.2.3)。

関連項目

参考文献

Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する