他の例: NDSolve

常微分方程式：基本的な用法

このコマンドは，0から3までの各点xで1次導関数と等価で，xがのときの値がである関数の数値近似を求める．

NDSolveはyをInterpolatingFunctionオブジェクトで置換する規則を返す．

In[1]:=

Out[1]=

Mathematicaはこの興味深い微分方程式の閉じた形の解を純関数として表現する．

In[2]:=

Out[2]=

これを普通の代数式に書き直す．この解はで特異点を持つ．

In[3]:=

Out[3]=

次に同じ式をDSolveの代りにNDSolveを使って近似的に解く．閉じた形式の解から分かるように，この警告は適切である．

In[4]:=

Out[4]=

数値に評価される他の関数と同じようにInterpolatingFunctionオブジェクトを使うことができる．境界条件を検証する3つの方法がある．

In[5]:=

Out[5]=

In[6]:=

Out[6]=

In[7]:=

Out[7]=

解のプロットも簡単にできる．次に示すのは導関数の逆数となる関数の近似値のグラフの見方である．

当然のように，このプロットは平方根関数とよく似ている．逆方向に解こうとしたらどうなるだろうか？

In[8]:=

Mathematicaは高次数の式も扱うことができる．次は3階の方程式の解のプロットである．

解を一意に決定するためには十分な初期条件を指定する必要がある．

In[9]:=

連立方程式の解も同様に求められる．2元連立方程式の場合は，位相プロットが解を視覚化するのに適した方法であることが多い．これは徐々に減速している振り子の動きを表している位相プロットである．

In[10]:=

In[11]:=

常微分方程式：オプション

この方法で余弦関数の微分方程式を解こうとしても区間が長過ぎて解けない．

In[12]:=

Out[12]=

MaxStepsオプションの設定を大きくする必要がある．予想されたように区間の最後では解が近くになる．

In[13]:=

Out[13]=

Methodオプションを使うと解の近似にNDSolveが使うメソッドを指定できる．デフォルトはAutomaticで，硬さによって自動的にBDF（後退微分法）かAdams多段階法に切り換える．Methodオプションにはさまざまな設定を行うことができる．方程式の種類によって適切なオプションは異なる．一般に硬くない方程式に使うことができ有効でもあるメソッドは，いろいろな次数の陽的ルンゲクッタ法を使うExplicitRungeKuttaである．問題によってはルンゲクッタ法の方が少ないステップ数で解を求めることができる．

Rösslerの方程式の場合は，ルンゲクッタ法はデフォルトの半分程のステップで済む．

In[14]:=

常微分方程式の非常に精密な解を得る方法に，WorkingPrecisionオプションの値を十分大きくすることが挙げられる．AccuracyGoalとPrecisionGoalはデフォルトでWorkingPrecisionの値の半分になっているので注意のこと．

WorkingPrecisionの値を大きくし過ぎないように注意する必要がある．というのも作業精度が上がるに伴って各演算操作や関数評価にかかる時間が長くなるばかりでなく，一般にステップ数も指数関数的に増大するからである．また，問題の係数が厳密な数であるか，少なくともWorkingPrecisionの設定値と同じくらい高くなっているように注意する必要がある．

これらのコマンドはWorkingPrecisionの値を変えてLorenzの方程式の解を比較する．その際，外挿法による硬さの切換えが使われる．これは，特に高精度の解の計算に有効である．

In[15]:=

In[16]:=

Out[16]=

In[17]:=

Out[17]=

In[18]:=

Out[18]=

最後の2つの解が指定されたWorkingPrecisionの桁数で与えられたからといって，それがその精度を持つというわけではない．実際，AccuracyGoalとPrecisionGoalのデフォルト値はWorkingPrecisionの半分であるが，この特定の問題に誤差を拡大する性質があるため（この問題はこれによって有名なのであるが）全体としての解の精度はこれほど高くはない．次は外挿法を使って実際の誤差を計算時間と対数のスケールの関数として表したグラフである．

精度が高くなるにつれ，時間が指数関数的に長くなっている．

In[19]:=

In[20]:=

常微分方程式：境界値問題

単純な線形境界値問題は，境界値の方程式を適切に構築することで解くことができる．方程式の次数がの場合，何らかの関数の組合せの値，あるいは個の点における導関数を知らなければならない．

この正規化された方程式は光学媒体の右端(x=)における波動効果を表している．

In[21]:=

これは，関数の値が3点で分っている3次方程式である．

In[22]:=

In[23]:=

微分代数方程式

微分代数方程式(DAE)とは同時に解かれるべき微分方程式と代数方程式が混ざった系のことである．Mathematicaは多くの種類の微分代数方程式系を解くことができる．

Mathematicaはデフォルトにより導関数を記号的に解いての形の1次の系を得ようとする．以下の例では平方根関数の2つの分岐により解が2つ与えられている．

In[24]:=

Out[24]=

しかし，導関数についての方程式の解は付近で十分な連続性がないので，どちらの分岐の解も問題を含んでいる．NDSolveを使って導関数に記号解を使わせる方法の代りに，オプションを使うことができる．この場合，NDSolveは方程式を満足させるために代数条件として効果的に単に数値的に解く．

In[25]:=

Out[25]=

この場合，NDSolveは区間全体に渡って解くが，数値解を求める際は1つの分岐しか追っていかないので解は1つしか与えられない．これやその他の微分代数方程式を解くために，MathematicaはLawrence Livermore国立研究所で開発されたIDA法を用いる．

もちろんについて微分することで方程式の解を求めることも可能である．こうすると大変なじみ深い単調和振動の式が得られる．

In[26]:=

Out[26]=

一般に，微分代数系の指数は微分方程式の系を得るために必要な微分回数に等しい．通常IDA法は指数の系を扱う．しかし，系を微分することで解の特徴を失うこともあれば数値誤差を導入してしまうこともあるので，注意が必要である．

連続する2つの解で代数条件がどのように満足されているかを比較すると興味深い．

In[27]:=

赤い方は代数条件からの直接の解を表している．この方法が許容率内に収まっているのは明らかである．一方微分した解の方では代数方程式は間接的にしか満足されておらず，偏差は時間が経つにつれて範囲から外れている．

NDSolveに与えた式が完全に代数的だった場合（つまりその式に導関数を含む項がなかった場合），NDSolveは化学反応速度論の例のようにIDA法を直接使って系を解く．

In[28]:=

Out[28]=

In[29]:=

In[30]:=

偏微分方程式：基本的な用法

最初の3つの例は1次の熱方程式を考える．固定境界条件を持つ，両端の温度を固定した絶縁ロッドでの放熱モデルである．

このコマンドは左端()を温度に固定し，右端()を温度に固定し初期熱の断面がxの2次方程式で与えられる式を解く．

In[31]:=

Out[31]=

常微分方程式だけの場合，解は従属変数にInterpolatingFunctionオブジェクトを代入する規則の形で返される．唯一目立つ相違点はInterpolatingFunctionオブジェクトが2次方程式であるという点である．つまり，NDSolveコマンドの変数と同じ次数の評価する引数が2つある点である．

境界条件はディリクレ(Dirichlet)とノイマン(Neumann)型の条件の一次結合となり，時間依存型の可能性がある．例えば，これらのコマンドを入力すると，左端で温度が正弦波状に変化する解をプロットする．

In[32]:=

偏微分方程式系の解を計算することもできる．これらのコマンドは放射線状と双曲線状の混合した系を解き，各従属変数の等高線プロットを別々に作成する．

In[33]:=

最終結果が大きくなるかもしれないというエラーメッセージが出される．これは空間近似の点の数が初期条件に基づいているためかもしれない．誤差推定は，これがの場合はAccuracyGoalオプションとPrecisionGoalオプションが典型的に満足されるように計算される．これの偏微分の場合のデフォルト値はなので，誤差に含まれるという余分な因子までが誤差が1%より小さいことを示している．

In[34]:=

In[35]:=

周期的な境界条件は数値解に便利なことがしばしばある．他の独立変数のすべての値に対する従属変数の解が領域の右端と左端で同じになるように指定することでNDSolveで周期的境界条件を付けて解くことができる．

例えば，周期解の興味深い解に非線形シュレーディンガー(Schrödinger)方程式がある．これらのコマンドは周期的初期条件（この場合はたまたまソリトン）を設定し，解を計算し，解の法と実数部，虚数部を示すのプロットを作成する．

次のセルを評価する．次にセルグループを閉じ，画像をダブルクリックしてアニメーションをスタートさせる．画像が入っているセルのセルブラケットを選択し「セル」メニュー項目の「グラフィックスのアニメーション化」を選んでもアニメーションを見ることができる．

In[36]:=

In[37]:=

偏微分方程式：オプション

偏微分方程式の解は線の数値解析法を使って解かれる．まず，時間変数として1つの変数が選ばれる．この変数の初期条件が指定されている必要がある．これ以外は空間変数として選ばれる．次に，空間変数について離散化が求められる．空間変数の離散化の結果はNDSolveが持つ常微分方程式を解くための多くの方法で解くことができる（一般に大きな）常微分方程式系である．一般にNDSolveに直接与えるオプションが時間変数の常微分方程式を解くために適用される．しかし，次のオプションは常微分方程式の空間的離散化と時間の系の解の両方に影響する．

オプションの中には空間的離散化特有のものもある．そのようなオプションには上記のオプションと一致して，解の時間的局面と空間的局面で別々の値を指定できるものもある．

これはデフォルト設定の可能なテンソル積グリッドオプションを示している．

In[38]:=

Out[38]=

これらのオプションの使い方に関する詳細はNDSolveの「拡張ドキュメント」を参照のこと．

離散化の局面で使われるデフォルト値は4次有限差分である．場合によってはこれ以外を選んだ方が，効率がよかったりより正確な解を与えたりすることもある．例として正弦Gordon方程式の，xとy方向への周期境界条件付きの2次元の一般化が挙げられる．

In[39]:=

Out[39]=

これはでの解のプロットである．

In[40]:=

このように周期境界条件付きの滑らかな解の場合，実際には最大空間差分の次数を制限するいわゆるpseudospectral導関数を用いると大変効率がよい．時間がの因数近くまで短くなっていることに注意．

In[41]:=

Out[41]=

In[42]:=

偏微分方程式：限界

NDSolveは線の数値解析法と呼ばれる偏微分方程式の解法を使う．これは1変数に分離して常微分方程式系を作り，この系をNDSolveに組み込まれた常微分方程式のメソッドを用いて解く．

この方法が作用するためには1つの変数に対しては初期関数が，もう1つの関数に関しては境界条件が指定されていなければならない．初期関数は常微分方程式系の初期条件を見付けるために使われる．矩形最大3辺で境界値と初期値を指定することができる．

この方法にはかなり大きなクラスの方程式が解けるという利点がある．しかし，これが解けない種類の方程式もある．

楕円形の問題は境界値が領域のすべての辺で指定されないと条件が不十分になるので，この方法では楕円形の問題の解は見付けられない．古典的な例題にラプラス(Laplace)の方程式がある．このコマンドを入力すると，この種のことをしようとしたときに何が起るかが分かる．一般的に条件が整わないと数値的な不安定さとなって現れ，作成されたプロットのスケールを見ると（警告）メッセージにもっと注意を払うべきであったことが分かる．

In[43]:=

現在はこの方法で解けないこれとは別の問題に，解に特異点が含まれるものがある．離散化は有限差分法を使って行われるので，正面が誤って解かれたり，完全になくなってしまったりする．例えば，バーガース(Burgers)の方程式は衝撃波形成を含む気体力学のある種の特徴のモデルである．これらのコマンドを入力すると有限差分法と正面のインタラクションに典型的に見られる振動を示すプロットが得られる．

In[44]:=

In[45]:=