微分方程式の数値解析パッケージ

Butcherニュートン・コーツ求積法
ガウスの求積法

NumericalDifferentialEquationAnalysisパッケージは，ブッチャー(Butcher)ツリー，ガウス求積法，ニュートン・コーツ(Newton–Cotes)求積法を使い，微分方程式を解析するための機能を統合したものである．

パッケージをロードする．

Butcher

ルンゲ・クッタ法はある種の常微分方程式を数値的に解くのに有効である．しかし，高次のルンゲ・クッタ(Runge–Kutta)法を導き出すのは決して簡単なことではない．この理由はいくつかある．まず，いわゆる次数条件を見付けるのが困難であるということである．次数条件は，任意の整数 の 次のメソッド（ここで は刻み幅である）における誤差を作るために満たされなければならないメソッドの係数の非線形方程式である．2番目に，この方程式を解くのが困難であるということがある．非線形であるだけでなく，これには一般に一意解が存在せず，多くのヒューリスティックや簡約が仮定される．最後に，組合せ的爆発の問題がある．12次のメソッドでは次数条件は7813にもなる．

このパッケージはこのうちの最初のタスクを行う．つまり満たすべき次数条件を見付ける．結果は未知の係数 ，， で表される． これによると， から へ進む 段のルンゲ・クッタ法は以下のようになる．

ここで，次が成り立つ．

行列の行にある要素の和は と に課せられた条件で繰り返される．これを認め，また，表記上での便宜のために，係数 を導入し，定義を以下のようにするのが一般的である．

この定義は行和の条件と呼ばれる，一連の行簡素化条件の最初のものである．

すべての について であるなら，このメソッドは陽的である．つまり，各 は前もって計算された値によって定義される．行列が厳密に下三角行列ではなかったら，この方は陰的で，各時間ステップの（一般に非線形）連立方程式の解が必要となる．対角線的に陰的な法ではすべての において である．

次数条件を表現する方法はいくつかある．段数 が正の整数で指定されていれば，次数条件は陽的項の和で表される．段数が記号として指定されていれば，次数条件には記号的な和が関係してくる．段数が指定されていなければ，次数条件は段数に無関係なテンソル表記で表される．行列 とベクトル および に加え，この表記にはベクトル も含まれる．これはすべての1で構成されている．この表記法には，段数 に無関係であることと特定のルンゲ・クッタ法と無関係であることの2つの際立った利点がある．

この理論についてより詳しく知りたければ，末尾に付記した文献を参照されるとよいだろう．

ai,jメソッドの の式における の係数
bjメソッドの の式における の係数
fi の便宜的表記
eベクトルの便宜的表記
tルンゲクッタ法での連続出力で使われる記号

ブッチャーの関数で使われる表記法

RungeKuttaOrderConditions[p,s]s 段 p 次のすべてのルンゲ・クッタ法が満足しなければならない次数条件のリストを与える
ButcherPrincipalError[p,s]s 段 p 次のルンゲ・クッタ法の，誤差のテーラー級数展開における，次数 p+1の項のリストを与える
RungeKuttaOrderConditions[p], ButcherPrincipalError[p]段数に依存しないテンソル表記の結果を与える

ルンゲ・クッタ法の次数条件と関連する関数

ButcherRowSum次数条件のリストに の行和条件を明示的に含むかどうかを指定する
ButcherSimplifyブッチャーの行と列の仮定の簡約を適用するかどうかを指定する

RungeKuttaOrderConditionsのオプション抄

次は，次数10までの各段階の次数条件の数を与える．組合せ的爆発に注目のこと．
次は，陽的に行和条件を含んだ三段一次ルンゲクッタ法で満足されなければならない次数条件を与える．
次は，陽的に行和条件を含んだ三段一次ルンゲクッタ法で満足されなければならない次数条件を与える．係数は形式文字で表され，明示的な表記では下付き文字になる．
次はすべての三段二次ルンゲ・クッタ法で満足されなければならない次数条件である．ここには行和の条件は含まれていない．

次数条件の左辺に含まれる和は記号的なまま残され，段数が記号引数として残された場合は展開されない点に注意したい．これにより，高次多段階法の結果が大幅に簡素化される．段数を全く指定しなければ結果はもっとコンパクトになり，答はテンソルの形で与えられる．

次はすべてのs段二次メソッドで満足されなければならない次数条件である．
sを3で置き換えるとRungeKuttaOrderConditions[2, 3]と同じ結果が得られる．
次は，すべての二次法で満足されなければならない次数条件である．ここではテンソル表記が使われている．ベクトルeは長さが段数である1のベクトルである．
テンソル表記も同様に，完全な条件を与えるように展開することができる．
これらはすべての三次法における第一誤差係数である．c2はベクトルであり，その成分はベクトルcの成分の二乗に等しい．
これは が0に近付く際の限界にあるすべての三次法の局所誤差の境界で，常微分方程式の影響を除去するために正規化してある．
次は一段四次のすべてのルンゲ・クッタ法が満足しなければならない次数条件である．この次数条件を満たすことはできない．すべての条件を満足するためには未知のものが十分にあるように，より上の段階である必要がある（第2引数をもっと大きくしなければならない）．
解は存在しない．
RungeKuttaMethod次数条件を求めているルンゲ・クッタ法のタイプを指定する
Explicit次数条件が陽的ルンゲ・クッタ法のものになるように指定するオプションRungeKuttaMethodの設定
DiagonallyImplicit次数条件が対角陽的ルンゲ・クッタ法のものになるように指定するオプションRungeKuttaMethodの設定
Implicit次数条件が陰的ルンゲ・クッタ法のものになるように指定するオプションRungeKuttaMethodの設定
$RungeKuttaMethod値をExplicit，DiagonallyImplicit，Implicitのいずれかに設定できる大域変数

RungeKuttaOrderConditionsおよび関連関数でのルンゲ・クッタ法のタイプのコントロール

RungeKuttaOrderConditionsや関連する特定の関数にはRungeKuttaMethodオプションがある．このオプションのデフォルト値は$RungeKuttaMethodである．一般的には， $RungeKuttaMethodをImplicit，DiagonallyImplicit，Explicitのいずれかに設定することで，ルンゲ・クッタ法が考慮されるように決めたいところだが，このオプション設定を指定したり，個々の関数についてデフォルト値を変更したりすることもできる．

次は，二段三次の対角陰的なすべてのルンゲ・クッタ法が満たさなければならない次数条件である．
対角陰的なメソッドを得るこれの代りの（しかしこれほど効果的ではない）方法に，上三角要素を0で置換してaを強制的に下三角にするものがる．
次は二段三次のすべての陽的ルンゲ・クッタ法が満足しなければならない次数条件である．この次数条件の矛盾点は，段数が次数よりも小さいときにすべての陽的ルンゲ・クッタ法に当てはまる結果を持つメソッドが不可能であることを示唆している．
ButcherColumnConditions[p,s]s 段 p 次，あるいはそれまでの条件を簡約する列を与える
ButcherRowConditions[p,s]s 段 p 次，あるいはそれまでの条件を簡約する行を与える
ButcherQuadratureConditions[p,s]s 段 p 次，あるいはそれまでの求積条件を与える
ButcherColumnConditions[p], ButcherRowConditions[p] 等 段数に依存しないテンソル表記で結果を与える

ルンゲ・クッタ法の位数条件に関連した他の関数

Butcherは，次数条件の数と複雑さはいわゆる仮定の簡約を適用することで高次において非常に簡単にできることを示している．例えば，十分な数の行の仮定の簡約と列の仮定の簡約，それに求積タイプの次数条件を適用するとこのような簡約が行える．RungeKuttaOrderConditionsのButcherSimplifyオプションはこれらを自動的に決定するのに使うことができる．

次は四次までの列を簡約する条件である．
次は四次までの行を簡約する条件である．
次は四次までの求積条件である．

ブッチャーの形式では，木は基本的なものである．木によってルンゲ・クッタ法のベキ級数展開における導関数と係数に対するこれに関連した次数条件の両方が与えられる．このパッケージはブッチャーツリーに関連した多くの関数を提供する．

fブッチャーツリーの表現で使われる基本的シンボル
ButcherTrees[p]p 次ルンゲ・クッタ法のための木のリストを次数で分けて与える
ButcherTreeSimplify[p,η,ξ]p 次までの求積条件，η 次まで行の簡約条件，ξ 次まで列の簡約条件がすべて当てはまると仮定して，ブッチャーの仮定の簡約によって約分されていない p 次までの木の集合を与える．結果は最初に消えない木から次数によってまとめられる
ButcherTreeCount[p]p 次までの木の数のリストを与える
ButcherTreeQ[tree]木あるいは木のリスト tree が有効な関数のシンタックスを持っている場合にはTrueを，それ以外ではFalseを与える

ブッチャーツリーの構成と列挙

これはすべての三次法に必要な木である．木は基本的なシンボルfを使って関数的形で表されている．
次は2つの木のシンタックスの有効性を検証する．ブッチャーツリーは関数fの乗算，ベキ乗，あるいは展開で構成されなければならない．
次は10次までの各次数における木の数を評価する．結果はOut[2]に等しい．しかし次数条件や木を実際に構築しないので，計算ははるかに効率的である．
10次までの各次数において必要とされる木の合計数を計算するのに，前の結果を使うことができる．

行と列の簡約の仮定を使ったメソッドの条件の数は段数によって決まる．これらの仮定が当てはまるとき，ButcherTreeSimplifyは既約ではないブッチャーツリーを与える．

これは，四次までの求積条件と，一次の行と列の簡約の仮定が当てはまるとき，四次法に必要な追加的な木を与える．結果は四次の単一の木である（これは単一の四次条件に対応する）．

木や森をグラフィカルに可視化できると便利なことが多い．例えば木を描画することで洞察が得られ，これによってルンゲ・クッタ法の構築が助けられることがある．

ButcherPlot[tree]木 tree のプロットを与える
ButcherPlot[{tree1,tree2,…}]森{tree1,tree2,…}にある木のプロットの配列を与える

Butcherツリーの描画

ButcherPlotColumns木のリストのGraphicsGridプロット中の列の数を指定する
ButcherPlotLabelプロット中の分岐点のラベルとして使うプロットラベルのリストを指定する
ButcherPlotNodeSizeプロット中の木の分岐点のスケーリング因子を指定する
ButcherPlotRootSizeプロット中の各木の根をハイライトするためのスケーリング因子を指定する．ゼロ値は根をハイライトしない

ButcherPlotのオプション

これは四次までの木をプロットし，これにラベルを付ける．

ブッチャーツリーを生成・描画することに加え，これらを計測・操作するための多くの関数が提供されている．これら関数の重要性に関する詳しい説明は，ブッチャー [1]をご覧いただきたい．

ButcherHeight[tree]木 tree の高さを与える
ButcherWidth[tree]木 tree の幅を与える
ButcherOrder[tree]木 tree の次数，すなわち頂点の数を与える
ButcherAlpha[tree](m,n)が辺ならば m<n となるように，完全に秩序立ったラベルの集合で木 tree の頂点にラベル付けする方法がいくつあるかを与える
ButcherBeta[tree]根にはラベルが付かず，頂点には1つおきにラベルが付くように，ButcherOrder[tree]-1種類の異なるラベルで木 tree にラベル付けをする方法がいくつあるかを与える
ButcherBeta[n,tree]葉にはすべてラベルが付き，根にはラベルが付かないように，n 種類の異なるラベルで頂点のうちの n 個にラベル付けをする方法がいくつあるかを与える
ButcherBetaBar[tree]tree 根を含むすべてのノードにラベルが付くように，ButcherOrder[tree]ラベルでラベル付けをする方法がいくつあかるかを与える
ButcherBetaBar[n,tree]すべての葉にラベルが付くように，頂点のうち n 個に n 種類の異なるラベルでラベル付けをする方法がいくつあるかを与える
ButcherGamma[tree]木 tree の密度を与える．密度の逆数は tree によって課せられた次数条件の右辺である
ButcherPhi[tree,s]木 tree の重みを与える．重みΦ(tree)は tree によって課せられた次数条件の左辺である
ButcherPhi[tree]テンソル表記を用いてΦ(tree)を与える
ButcherSigma[tree]木 tree のそれ自身の同型写像の対称性群の次数を与える

ブッチャーツリーに関連するその他の関数

これは木f[f[f[f] f^2]]の次数を与える．
これは木f[f[f[f] f^2]]の密度を与える．
これはs段階法にf[f[f[f] f^2]]で与えられる基本的な加重関数を与える．
下付き表記はフォーマットを整える装置であり，下付き記号は単に指標を付けられた変数NumericalDifferentialEquationAnalysis`Private`$iに過ぎない．

Solveとその関連関数を使って次数条件の解を得ることもできる．Sofroniou [6]にはこのパッケージを使ったルンゲ・クッタ法の構築に関連するさまざまな問題が載っている．この論文にはButcher.mに使われているアルゴリズムの詳細とその応用についても記されている．

ガウスの求積法

Wolfram言語の関数NIntegrateはかなり高度なガウス・クロンロッド(Gauss‐Kronrod)則をメソッドの一つとして用いている．NumericalDifferentialEquationAnalysisで提供されるガウスの求積法では，通常のさほど高度ではないガウス求積法の背後にある理論を簡単に学ぶことができる．

ガウスの求積法の基本的な概念は，被積分関数の値の線形結合としての積分が特定の点で評価した場合は値を近似するというものである．

選ぶことのできる の自由パラメータ（横座標 と重み の両方）があり，積と和の両方が線形操作であるので， 以下の次数のすべての多角形にあてはまる公式を作ることができると考えられる．選択的な横座標と重さが何かを知ることに加え，近似誤差がどれくらい大きいのかを知っていることも望ましい．このパッケージを使うことでこの2つの問いの両方に答えることができる．

GaussianQuadratureWeights[n,a,b]区間 a から b において，ペアのリストを求積法の機械精度まで与える
GaussianQuadratureError [ n , f , a , b ]
機械精度の誤差を与える
GaussianQuadratureWeights [ n , a , b , prec ]
ペアのリストを精度 prec まで与える
GaussianQuadratureError [ n , f , a , b , prec ]
誤差を精度 prec まで与える

ガウスの求積法の公式を求める

これで区間における5点のガウスの求積法の横座標と重みが与えられる．
次はこの公式における誤差である．残念なことに未知の点における の10次導関数が含まれているので，誤差自体が何かは分からない．
区間が長くなると誤差が急速に小さくなるのが分かる．

ニュートン・コーツ求積法

Wolfram言語の関数NIntegrateはひとつの方法としてかなり高度のガウス・クロンロッド則に基づいたアルゴリズムを用いる．積分の公式には他のタイプのものもあり，それぞれ独特の利点を有する．例えば，ガウス公式では，奇数間隔の横座標における被積分関数の値を用いる．これを用いて等間隔の横座標における表形式で表された関数を積分しようとしてもあまりうまくいかない．そのような場合にはニュートン・コーツ公式を用いるとよい．

ニュートン・コーツ公式のもとになっている基本的な考え方は，等間隔の点で被積分関数を評価した値の線形結合としての積分の値を近似しようというものである．

これに加え，和に終点を含めるかどうかという問題がある．終点を含んだ求積式は閉じた式と称される．終点を含まない求積式は開いた式と称される．開いた式では最初の横座標をどこに置くかという点が曖昧になる．このパッケージでの開いた式では，最初の横座標は下の終点から半ステップのところに置かれる．

選択可能な 個の自由パラメータ（重さ）があり，積分と和の両方が線形操作なので，およそ 未満の次数のすべての多項式について正しい式ができると予想できる．重さがどのくらいかを知ることに加え，近似における誤差がどのくらいになるかを知ることも望ましい．このパッケージを使うと，両方の答が与えられる．

NewtonCotesWeights[n,a,b]区間 a から b における求積のために，のペアが n 個含まれるリストを与える
NewtonCotesError[n,f,a,b]式の誤差を与える

ニュートン・コーツ積分公式を求める

オプション名
デフォルト値
QuadratureTypeClosed求積の型．OpenあるいはClosed

NewtonCotesWeightsとNewtonCotesErrorのオプション

次は区間における5点の閉じたニュートン・コーツ公式の横座標と重さである．
次は公式の誤差である．未知の点における の六次導関数が関連しているので，誤差それ自体がどんなものなのかは分からない．
区間が長くなるにつれて誤差は急速に小さくなる．
次は区間における5点の開いたニュートン・コーツ公式の横座標と重さである．
次はこの式の誤差である．