A1 波動方程式（一次元）パート2
パート1ではソフトウェアのコマンドや簡単なコードで解が得られるものを紹介した。パート2では、なぜそのような解が得られるのか？数学的に教育的に追ってみた。ここまでくるとLive math makerでは力不足なので、Mapleで仕上げることになった。最終的な数式は式（61）のすぐ上に書いてある。
----------------------------------------------------------------------------------
----------------------------------------------------------------------------------
だいぶページが嵩張った。ODEの計算はMapleコマンドを使ったけど、ここまでの数学的な論理を通してアプリケーションの結果が得られる事を、一度確認しておくことも良いでしょう。
参考になった図書は
・小出氏の図書（パート1で紹介）
・川村氏の図書
・フーリエ級数はクライツィグ氏の図書ですね。
偏微分方程式（以下PDE）
独立変数が複数になるので、ODEに比べ複雑さは桁違いだ。それでも適用アプリケーションが存在するので、できるだけ紹介していきたい。
方程式には1階、2階という階数があるが、自然界の法則を表現するPDEには2階が非常に多い。何故多いのかは上手く説明できないが、私的には2次関数の考え方をアナロジー的に利用しているからと思っている。これまで3次関数の延長線での応用など聞いた事が無いし、線形代数の公式も3次になると途端に面倒になるし、2から3へのステップアップは「数学の専門分野」の話になるのかもしれない。まずは整理されている2階線形からスタート。
2階線形PDEの分類
現在、様々な参考書では次のように扱われている。
2階線形PDEは、独立変数が２つの場合は、
式（1）で表せる。ここでuはxとyの関数で、A～Gはxとyの関数または定数である。最終的には表1のように大きく3種類に分類される。
|
判別式 |
根 |
型 |
代表PDE |
|
|
B 2-4AC > 0 |
2実根 |
双曲型 |
|
|
|
B 2-4AC = 0 |
重根（実根） |
放物型 |
熱伝導方程式 |
|
|
B 2-4AC < 0 |
共役な2複素根 |
|
||
表1
これ以降、3種類のPDEについて、例を挙げていく。ODEの説明では解析解と非解析解に分けたが、PDEは具体的に初期条件などを与えて全て非解析解（級数解、数値計算解）とし、出来るだけ視覚的にも工夫を凝らして説明を進めていく。
A1 波動方程式（一次元）パート1
一次元波動方程式を調べていくと、無限に長い弦の振動から、２つの任意関数を含む一般解「ダランベールの解」にたどり着く（式（2））。
更に工学的利用を考えて初期条件を考慮すると、以下の「ストークスの波動公式」が得られる。
初期条件は、
ここで、αとβはxを独立変数とする任意関数である。数学的な振る舞いを調べるにはこれらの形が都合良いらしい。
工学的には初期条件だけでなく境界条件も必要な場合が多い、これにアプローチするにはフーリエの方法と言われ、変数分離やフーリエ級数などを用いて進められる。ここまでは解析解と言って良い。
計算例
（1）数式処理ツールMapleの計算結果を以下に示す。PDEなのでpdsolveコマンド類を使っている。初期条件はu(x,0)=?、1回微分はD[2](u)(x,0)=?の様に入力する決まり（1番目の引数で微分する場合はD[1]、2番目の引数で微分する場合はD[2]）になっている。#を先頭に付けると、それ以降は計算実行に影響を与えない説明文になる。unionは初期条件と境界条件を連携させるために入れてある。
解plotを見ると初期条件は弦の真ん中を上に引張ったイメージである事が解るだろう。MapleにはPDE数値計算手法が色々と用意されているが、時間tの2階微分を用いる波動方程式の解法は残念ながら有限差分法のみである（それでpdsolve(,,,numeric)を使った）。それでもプログラミングを利用するまでもなく解を可視化できる所が凄い。
------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------
（2）次はWolfram CloudのMathematicaだ。
今回は計算に必要なメモリーが多く、無料でクラウド計算できる範囲を超えてしまい結果が得られなかった。取り敢えず、果敢に取り組んだ事だけ見てもらおう。NDコマンドを利用した。Mathematicaを使用できる環境にいる方は是非試してほしい。
-------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------
（3）次はFlexPDEという名前のツールで、PDE Solutions Inc.のHPからダウンロードできる。その名の通りPDE数値解析専門ツールだ。Liteバージョンなら、制限はあるものの試用可能だ。計算原理は有限要素法（以下FEM）なので、機械系の解析がメリット多い感じだ。FEMに乗せられれば、構成方程式（数学記号に似ている）を選ばない。その構成方程式と共に初期条件や境界条件を入力する。様々なPlotコマンドがあるので、最後に試行錯誤しながら解をPlotしていく。数値解析ツールなので、得られるのは数値データだけで数式は得られない。
初めての紹介なので、プログラムコードの補足を加える。赤い字は約束された記号である。TITLEは解く問題の名前。SELECTでは計算誤差を1/10000以下にするように指示できる。VARIABLEではその名の通り変数を指定する。FlexPDEは時間の2階は積分できないので、本問題では二つのPDEに分けて連立して解くしかない。DEFINITIONSでは構成式中の定数などを指定する。本問題ではcとLそして計算スタート時間t0だ。INITIAL VALUESでは初期値を指定する。本問題ではvはuを時間tで微分したものである。EQUATIONSでは、その名の通り構成方程式を指定する。本問題では先に説明したように時間の2階微分を1階に落とすために連立式にした。BOUDARIESは境界条件を指定する。本問題では大きさはy方向、時間軸はx方向にして二次元で表現する。FEMはこの二次元平面を要素分割して計算を進めていく。start(0,L)は計算スタート地点の座標で、このスタート地点から二次元平面を反時計回りに計算していく習わしである。最後のLine to closeで二次元平面を閉じている。座標間に初期条件や境界条件を指示する。valueはディリクレ形境界条件で、naturalはノイマン型境界条件だ。TIMEは計算する時間範囲を指定している。PLOTSは様々なplotコマンドを入力するエリアだ。最後はENDでプログラムを閉じる。
プログラムコードを下に、続けて下から5行目のhistory(u)の結果を載せる。
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------------
このplot（少し歪んでいるが）とMaple出力plotを見比べて欲しい。x=0.5の値が時間経過で振動している姿がソフトウェアが異なっても見事に再現されている。
今回、参考にした図書を挙げておこう。
・小出氏の参考書は具体例が多く重宝する。
Mathematicaに関しては、藤井氏のテキストを参考にさせてもらった。記述は懇切丁寧で図も多い。
FlexPDEにに関しては、下のGunnar氏のテキストを参考にさせてもらった。この本はPDE solutions Inc.のHPにリンクがあるので、辿ると購入可能だ。
BB非解析解
1. 数値解析解
AAでは解が陽的に数式で得られるケースを挙げたが、自然界の振舞を表すODEはそれだけではない。非線形現象のためむしろ数式で表せない方が圧倒的に多い。そんな時は数値解析手法に頼る事になる。つまりAAで頑張っても解けない場合は数値解析しましょうという訳である。このブログの性格上、コーディングは最終手段なので、できる限りソフトウェアの搭載機能を利用していきたい。
数値解析手法は大きく次の３つに分けられる。陽解法は数値解析の基礎である。とにかくODEをどのようにコンピュータに計算させるのかを理解するために通らなければならない入口だ。本ブログではルンゲクッタ法を使った計算例に限定する。陰解法は非線型ODEがターゲットである。ODEの性質により著しく解きにくい（stiff）場合などは予測子—修正子法が使用される。
（i）陽解法
・オイラー法：微分方程式の数値解法の一番手は、どのテキストを眺めてもオイラー法から始まる。解りやすいのが一番だ。しかし一次精度のため計算刻みを相当ゼロに近づけないと計算誤差が累積していき、近似解精度が保てなくなる。理解しやすいがコスパが悪すぎる。
・ルンゲクッタ法：2次は修正オイラー法とも言われる（2次のテイラー展開利用）。一般的なのは4次のテイラー展開を利用したもので、よくRK4と表示される。5次のテイラー展開利用したものも存在する。当然次数の高い手法ほど時間がかかるが精度は高くなる。
(ii) 陰解法
・後退オイラー法：一次精度なので、通常は使用する必要はない。
・クランク・ニコルソン法：二次精度。
(iii) 予測子—修正子法：予測子に陽解法を、修正子に陰解法を用いる。計算回数が多くなると振動を起こす。
・Milne
・Adams-Moulton法
計算例
（1）数式処理ツールMapleの計算結果を以下に示す。
大きく二つに分けてある。一つはグラフ的解法。withコマンドでパッケージDetoolsをロードして、DEplotを走らせると解の流れが見られる。計算しなくてもODEの解の全体像が俯瞰できる、凄いよね？ 初期条件を入力することで、当然それを通過した解曲線がプロットされる。ここで大事なことは、初期条件により解の結果が異なる、つまり初期値に敏感なODEが存在する事を認識することだ。
二つ目はMapleに搭載された様々なルンゲクッタ法を使用した例である。rk3、rk4、rkf45、ck45、dverk78、rosenbrock。計算結果はリスト表示しているが、その違いは判別できないくらい小さい。コーディングする前に色々試せるということ。
------------------------------------------------------------------------------------
図1
------------------------------------------------------------------------------------
（2）次はWolfram CloudのMathematicaの計算結果だ。NDSolveコマンドでODEの数値解析ができる。計算結果はオブジェクトInterpolatingFunctionに格納される。そのデータをEvaluateコマンドを使用してplotする。
二つ目はMethodにExplicitRungeKuttaを指定したケースである。マニュアルにはImplicitRungeKuttaも記載されている。
両プロットとも違いは分からないほどだ。
------------------------------------------------------------------------------------
図2
--------------------------------------------------------------------------------------
（3）最後はLiveMathによる計算結果だ。
数値解析する場合はまず計算したいODE全体をセレクトする。画面上部の計算メニューからTabular->微分方程式の数値解法をセレクトする。すると図4の様なODEソルバー画面が現れるので、画面上部のボックスからルンゲ-クッタ法をセレクトする。また独立変数の範囲や計算間隔を入力すると、計算結果の個数を合計欄に表示してくる。間隔を小さくすれば、それだけ計算合計が増えるのでプロットは綺麗になる（時間は掛かる）。初期値は初期値画面で入力するが、従属変数の値しか入れられない（独立変数の初期値は範囲の最初の値になる）。okボタンを押すと計算結果が表（図3のプロットのすぐ上にある行）で格納される。そのデータをプロットしたものが図3である。
-------------------------------------------------------------------------------------
図3
図4
図5
-------------------------------------------------------------------------------------
今回の3種類のソフトウェアはコマンドこそ異なるが、計算結果が非常に近いことが理解できるはずだ。
実は、AAの解析解グループにはグリーン関数法・たたみ込み積分法などの手法も存在するが、どのような場合にこれらを使用しなければならないのか、不勉強で残念ながら僕では伝える力が無い。また数値解析の前段階には近似解法として、テイラー級数やフーリエ級数を利用した手法もある。詳細は参考書を挙げるだけで勘弁願いたい。
べき級数解は、
フーリエ級数解は、
数値計算（ルンゲクッタ法含む）については素晴らしい成書が随分あるが、新しい本として、以下を紹介したい。
AA解析解
3. 線形ODE（ラプラス変換）
2階定数係数ODEを、積分変換法の一つであるラプラス変換（17式）を用いて解こう。
ここで、f(t)は0 ≦ t < ∞で定義された連続関数、tはC級（実数および複素数）である。
ラプラス変換は積分変換を導入することによって、微分演算を代数的変換（変換空間内で）に置き換えたテクニックである。関数 f(t)に制限は有るものの、機械的に解が求められる手法に感動すら覚えるのは僕だけでは無いだろう。特に自動制御の分野では、時間空間 f(t)を複素空間 F(s) に写像し演算子sで解析する手法が確立されている。
解法の大まかな流れは、「ODE—＞ラプラス変換—＞代数方程式—＞代数方程式の解—＞ラプラス逆変換ー＞ODEの解」となる。
ラプラス変換やラプラス逆変換の公式は適宜、手持ちの教科書を参考にされたい。
計算例
（1）数式処理ツールMapleの計算結果を以下に示す。dsolve(?,?,method=laplace)を使うと、一発で回答してきた。初期条件はu(0)=?、2回微分はD(u)(0)=?の様に入力する決まりになっている。
------------------------------------------------------------------------------------
----------------------------------------------------------------------------------
（2）次はWolfram CloudのMathematicaの計算結果だ。見ての通りだいぶコマンドが異なるが、計算手順は判りやすい。%は直前の計算結果を利用することを意味する。最後のExpandは展開のコマンドである。
----------------------------------------------------------------------------------
-----------------------------------------------------------------------------------
（3）次は、オープンソースの数式処理ツールwxMaximaによる計算結果だ。今回はメニュー画面が使用できなかったのでコマンド入力とした。assume(?)はデータの形を仮定するコマンドである。LY1:は式の名前を指定している（今回は見易くするため）。以降は:の前の名前のみで進められる。%はMathematicaと同じように直前の計算結果を意味する。計算過程の理解しやすさもMathematica同様だ。
----------------------------------------------------------------------------------
--------------------------------------------------------------------------------
（4）最後はLiveMathによる計算結果だ。まずはラプラス変換パッケージをコピーする。ODEをラプラス変換。y(0,0)は初期値、y(1,0)は1回微分の初期値であり、それを代入する。部分分数展開を行う（結構これが大変で、結局はMapleの結果を引用させてもらった）。ラプラス逆変換を行う。求められたyの形をプロットしてみた。
-----------------------------------------------------------------------------------
-----------------------------------------------------------------------------------
ラプラス変換は微分の公式に陽に初期条件（特に時間経過を表現する tの場合はコーシー条件とも言われる）が現れるので、初期条件が分かっていれば、途中で代入（式はシンプルになる）しても、最後に代入（全体像が把握できる）しても構わない。微分方程式の数値解法の一つ有限要素法の定式化の途中においても、このように「陽に境界条件が現れる」という現象がある。興味深いものだ。
ラプラス変換は減衰する指数関数で分解される、つまり色々な時定数で減衰する関数
の重ね合わせで表現される。これは積分範囲が -∞では発散することを示唆するので、よってラプラス変換は、積分範囲の片側が有限な領域で定義された関数についてのみ適用になることに注意してほしい。
ラプラス変換などとは異なる手法にミクシンスキー演算子法なるものがある（ラプラス変換できない関数に役立つかもしれない）。
部分分数展開について、思っていたより難儀した。それでも整理されたテキストを見つけたので参考まで紹介する。部分分数に展開したい式の因子の種類によって大きく4つに分けられる様だ。
・非重複因子sーa
・重複因子(s-a)^m
・非重複複素因子(s-a)*(s-abar)：aとabarは共役複素数。
・重複複素因子{(s-a)*(s-abar)}^2：aとabarは共役複素数。
参考になったのは、E.Kreyszigの「技術者のための高等数学」第3巻「フーリエ解析と偏微分方程式」のラプラス変換の章だ。
今回新たに参考にした図書を以下に載せるので、興味ある方はどうぞ！
ラプラス変換公式が教科書にない場合は、これを参考にしています。
AA解析解
2. 線形ODE
のうち、応用例の多い2階定数係数ODEを考える。式（4）を2階に、またその係数を定数に置き換えると、
ここでa、bは定数である。
（i）初めに次のような2階定数係数同次方程式から始める（右辺＝0）。
式（5）の基本解が三角関数や指数関数で表せるパターンが多いことが知られているから、式（6）のy(x)にe^λxを予想して代入してみる。ここでλはC級（複素数全体）である。
整理すると、
式（8）が成り立つためには、
式（9）（特性方程式という）が成り立てば良い（e^λxは非ゼロだから）。ここで見覚えのある2次方程式の解の公式の判別式D（平方根内部）を借用すれば、その値によって２つの特性解が異なってくるはずである。式（6）の一般解をまとめてみると、
・D > 0（特性解が異なる実数λ1、λ2）
・D = 0（特性解が１つの実数で、重解λ1）
・D < 0（特性解が複素共役で、λ1＝ p + qi、λ2＝p - qi）
(ii) 次に2階定数係数非同次方程式を解く事を考える（定数変化法）。
式（13）の一般解は、同次方程式（6）の基本解２つy1(x)、y2(x)を用いて、次のように表せる。
計算例
（1）数式処理ツールMapleの計算結果を以下に示す。2階になって難しさが増えたはずだが、同次・非同次に関係なく、ほぼ一発で回答してきた。非同次EQを解く方法は大きく三つある。積分因子法、未定係数法、定数変化法である。Maple内部はよく解らないが、右辺の式で解法パターンを分けているはずだ。本例題の右辺は多項式なので、おそらく未定係数法だろう。
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
（2）次はWolfram CloudのMathematicaの計算結果だ。Mapleと同様に一発回答だ、素晴らしい。気を付けなければならないのは積分定数だ。記号は同じでも、同次式と非同次式では中身が異なる。
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
（3）次は、オープンソースの数式処理ツールwxMaximaによる計算結果だ。今回の積分定数は%k1、％k2になってきた。%cとどのように使い分けているかは不明。
--------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------
（4）最後は教育的な（？）LiveMathによる計算結果だ。積分定数は最終的にA1、A2に置き換えてある。ページ数が大きいので、適宜分けて掲載する。
右辺の関数形に限定されない定数変化法でトライしたが、手計算だとこんな感じになってしまう。やはり、右辺の関数を見て適切な未定係数法を選択すべきだね。ロンスキー行列式や不定積分の公式は説明なしに使っているので、気になる人は教科書見てね。
とにかく、上に紹介したツールがどの位優秀なのかが理解出来るはずだ。
------------------------------------------------------------------------------------
同次（homogeneous）
非同次（inhomogeneous）はページが膨大なので分けて掲載する。
part1
part2
part3
---------------------------------------------------------------------------------------
今回学ぶ事は「ODEの手計算チェックは最後にしよう」だ。先見事項は一発回答してくれるフリーのツールを複数試して同じ回答になることを確認することだ。
手計算する際に、様々な公式が必要になる。今回は公式集を紹介する。
微分方程式は大きく２つ、独立変数が1個の常微分方程式と、独立変数2個以上の偏微分方程式に分けられる。読み方を簡単に、常微分方程式をODE、偏微分方程式をPDEと表記する。以降この表記で進める。
自然界の現象を微分方程式で表現しようとすると、多くがPDEになる。つまり独立変数が単一では無いケースがほとんどという事だ。よく考えればその通りだ。座標空間3次元（x,y,z）、これだけで2変数関数のPDEになる．時間経過も要因に入れると、3変数関数のPDEだ。1変数関数y = f(x)のような結果であれば平面にプロットできる。z=g(x,y)のような2変数関数だと3Dプロットになる。3変数関数h(x,y,t)などは、3Dアニメーション表示するしかない。
まずはODEからスタート。
AA解析解
このタイプは結構頻繁に出くわす。g(y)がゼロでない仮定のもとに、一般解は式(2)のようになる。ここでCは積分定数。
例題として放熱現象を挙げる。沸かし終えた風呂の水温の時間変化率は、その時の浴室温度との差に比例する。つまり浴室温度が10℃の時、時刻tにおける風呂水の温度T℃は、
のようなODEになる。ここでkは正定数とする（本例題では、浴室の構造や風呂桶の容量・材質などの細かい条件がkに含まれていると仮定する）。
（1）数式処理ツールMapleの計算結果を以下に示す。流石に一発で答えてきた。一々変数分離とか指示しなくていい。
-----------------------------------------------------------------------------------
----------------------------------------------------------------------------------
（2）次は、Mapleと2大巨頭のMathematicaの計算結果だ。Wolfram Cloudならばフリーで使えるのは前にも書いた通り。ソフトウェアにより入出力表現は若干異なる。これも一発で答えてきた。当然c1は積分定数だ。
----------------------------------------------------------------------------------
------------------------------------------------------------------------------------
（3）次は、オープンソースの数式処理ツールの老舗wxMaximaによる計算結果だ。思っていたよりインターフェースが洗練されていてお勧めだ。
------------------------------------------------------------------------------------
-------------------------------------------------------------------------------
（4）最後はLiveMathによる計算結果だ。上で紹介したツールと違い、自分で教科書を見ながら実行していくタイプだ。計算過程が全て残るので教育的にはいい。積分定数は-c2-c1を-c11に置き換えてある。
-------------------------------------------------------------------------------
----------------------------------------------------------------------------------
アプリケーションにより一長一短はあるが、LiveMathのようなツールは手計算に近く、論理的思考の記録として価値がある。技術的な報告書のオモテにはあまり出ないが、報告書フォルダには検討資料（手計算の殴り書きメモもある）として必ず入れている。
解りやすい参考書を一冊挙げておく。
今や、線形代数の応用分野は理系だけではなく、おそらくコンピュータの使用される分野全てと考えていいだろう。メイン分野は以下だ。
・構造解析
・流体解析
・制御系
・暗号理論
教育機関で特に数学系以外なら、ベクトル、行列、行列式、連立一次方程式、固有値までの説明は受けるはずだ。しかもSarrusの方法が使用できる2次、3次までがメインだ。4次以上、つまり一般化に近づくと、Sarrusの方法が使用できない事や、様々な定理や公式が見辛くなって理解を妨げる事、計算コストを抑えるために特有の数値計算手法を使用せざるを得ない事が、その理由と考えられる。エッセンスだけなら3次まででいいだろう、て考えてしまうのも仕方ない。当然、教科書も学部教養レベル・数学系・大規模計算では異なり、他の数学分野に比べ、一般産業と親和性が高く、未だに進化し範囲を広げ続けているのが線形代数だ。僕も数年に一冊は新しいテキストを購入している。
という事で「線形代数に関して、理系ツールの使用例を少し載せよう」と意気込んでいたが、自身があまり明るくない分野である事、先達の素晴らしいブログが存在する事、を理由にこの程度の文章だけになった、グスン。
マトリックスの計算で有名なソフトウェアはMATLABだ、オープン系ならScilabやGNU Octaveが存在する。
最後に手元にある新し目のテキストを紹介しておく。
