関数近似パッケージ

これでパッケージをロードする.
In[1]:=
Click for copyable input

パデ近似を使ったケチ表現による有理近似

EconomizedRationalApproximation[f,{x,{xmin,xmax},m,k}]
から の範囲に有効で次数を持つ,f のケチ表現による有理近似を与える

ケチ表現による有理近似

パデ近似は,展開の中心部付近では非常に正確であるが,中心から離れるにつれて誤差が急速に増大する.中心付近でのよいフィットを多少犠牲にすることをいとわなければ,指定の区間で全体的によくフィットさせることが可能である.これが他のタイプの近似で行われるものである.

ケチ表現の有理数近似とは,パデ近似から始めて,誤差における主係数を小さくするようにチェビシェフ多項式でパデ近似を摂動するという方法である.摂動により,消えた項が再現する.誤差の大きさは展開の中心部付近でわずかに増加し,このわずかな増加分は遠隔部の誤差が減少することによって埋め合せられる.EconomizedRationalApproximationでは展開の中心ではなく,近似を行う区間を指定する.区間の長さがゼロになるような限界において,ケチ表現の有理数近似はパデ近似に近付く.

における次ののパデ近似である.
In[105]:=
Click for copyable input
次は,パデ近似と実際の関数との差分をプロットする.近似は展開の中心部付近では非常によいが,中心から離れるにつれ誤差が急速に増加する.
In[109]:=
Click for copyable input
区間における の次数のケチ表現の有理数近似である.
In[6]:=
Click for copyable input
Out[6]=
実際の関数とケチ表現の有理数近似との差分を与える.一般に,ゼロになることのない小さい定数項を得ることさえあり得る.
In[7]:=
Click for copyable input
Out[7]=
これは上のパデ近似と比べると,誤差が全体的に散らばっている.
In[8]:=
Click for copyable input
Out[8]=
終点 における誤差は特に小さいわけではないが,PadeApproximantが与えるものよりもかなり小さくなっている.
In[9]:=
Click for copyable input
Out[9]=

有理補間とミニマックス近似

次の有理関数とは,m 次の多項式の k 次の多項式に対する割合のことである.有理関数は初歩の算術演算しか使わないので,簡単に数値的に評価することができる.分母の多項式を使って有理特異点を持つ関数を近似することができる.このため,有理関数はしばしば数値演算で関数を近似するのに役に立つ.

このような近似の方法はたくさんある.方法は近似の優秀性という概念の解釈の仕方によって異なる.どの方法も一定のクラスの問題で有効である.このパッケージを使って一般的な有理補間とミニマックス近似を計算することができる.関数PadeApproximantにはPadé近似を実行する機能が含まれている.

これに関連したクラスの近似問題で,近似関数でデータ点のセットを補間あるいはフィットするというものがある.このような場合は,組込み関数のFitInterpolatingPolynomialInterpolationを使うとよい.詳細に関しては「データの数値操作」を参照されたい.

RationalInterpolation[f,{x,m,k},{x1,x2,...,xm+k+1]}
次の有理補間を与える
RationalInterpolation[f,{x,m,k},{x,xmin,xmax}]
自動的に選ばれた点 で有理補間を与える

有理補間

指定の関数を有理関数で近似する方法のひとつに独立変数の値1セットを選び,この値において指定の関数と一致する有理関数を構築するというものがある.RationalInterpolationが行うのはこれである.

RationalInterpolationには2通りの使い方がある.独立変数の区間のみを指定した場合は,選択した次数の近似を有効に行えることを保証するような値が自動的に選ばれる.使用する一連の値のリストを明示的に与えることもできる.この場合は,次近似を求めるなら,独立変数の 個の値のリストを指定しなければならない.

次は0から2までの等間隔にある7点で,次の有理補間を与える.
In[36]:=
Click for copyable input
Out[36]=
次は関数と近似の相違である.終点に近付くにしたがって誤差が大きくなるのが分かる.
In[3]:=
Click for copyable input
Out[3]=
次では,Mathematica が自動的に補間点を選んでいる.
In[4]:=
Click for copyable input
Out[4]=
補間点は区間の終りに集まっている.通常このような場合は最大誤差が小さくなる.
In[5]:=
Click for copyable input
Out[5]=
オプション名
デフォルト値
WorkingPrecisionMachinePrecision使用する精度桁数
Bias0補間点を自動的に選択する際のバイアス

有理補間のオプション

RationalInterpolationx の値の領域を指定すると,補間点が自動的に選択される.オプションのBiasを使うと補間点を右または左に偏向させることができる.Biasの値はからまででなければならない.この値が正のときは点が右に偏向する.デフォルト値はBias->0で,この場合点は左右対称に置かれる.

点の位置を右側に偏向させると,右側で誤差が小さく,左側で誤差が大きくなる.
In[6]:=
Click for copyable input
Out[6]=
次はバイアスの影響を示したものである.
In[7]:=
Click for copyable input
Out[7]=

RationalInterpolationを使うと指定の関数と複数の点で一致する有理関数が得られる.これにより有理関数が,ある意味で,指定の関数に近いことが保証される.有理近似でのより強い必要条件は,すべての区間で有理関数が指定の関数に近くなることだろう.この種の有理近似はMiniMaxApproximationを使って得ることができる.この近似は,近似と指定の関数との相対誤差の最大値を最小にするためにこの名が付けられている.つまり,指定の関数 のミニマックス近似 は,検討中の区間におけるの最大値を最小にする指定の次数の有理関数になる.ミニマックス近似という語は,ここでのような相対誤差ではなく絶対誤差を最小化する有理関数を指すこともあるので注意のこと.

MiniMaxApproximation[f,{x,{xmin,xmax},m,k}]
から の区間で,次数f のミニマックス近似を与える
MiniMaxApproximation[f,approx,{x,{xmin,xmax},m,k}]
反復アルゴリズムを approx で始めてミニマックス近似を与える

ミニマックス近似

MiniMaxApproximationは反復スキームを使う.まず最初にRationalInterpolationを使って有理近似を構築する.この最初の近似は,次にRemesのアルゴリズムに基づくスキームを使ってより精密な近似を生成するために使われる.新たな近似の生成は相対誤差が小さくなるように外挿点を選別することで行われる.

MiniMaxApproximationは,最大誤差が生じる点のリストと有理近似と最大誤差の値からなるリストの2つの部分からなるリストを返す.この追加情報は,ユーザの便宜のためというよりも,最初からやり直さずに手順が再開できるという機能を与えるために提供されている.アルゴリズムが反復的で,MaxIterationsに達する前に収束しないと警告とともに不完全な答えが返されるので,この機能は便利である.

以下は最大誤差の点と望ましい補間,誤差の値を含むリストを返す.
In[8]:=
Click for copyable input
Out[8]=
次は有理近似を抽出する.
In[9]:=
Click for copyable input
Out[9]=
次は区間における近似中の相対誤差のプロットである.任意の極値の誤差を小さくしようとすると,他の極値における誤差が大きくなる.
In[10]:=
Click for copyable input
Out[10]=

MiniMaxApproximationは「相対」誤差の最大値を最小にしようとするので,問題としている区間に零点が含まれる関数のミニマックス近似は求められない.有理近似は関数の零点で正確にゼロにならなければならない.さもなければ相対誤差が無限大になってしまう.このような関数を扱うことは可能ではあるが,まず関数を割って零点を除外してから,掛け直して有理関数に戻さなければならない.

で相対誤差が無限大になる.
In[11]:=
Click for copyable input
Out[11]=
(x-Pi/2)で割ると零点が相殺され,近似計算が問題なく行える.
In[12]:=
Click for copyable input
Out[12]=
近似に(x-Pi/2)を掛けるとCos[x]のミニマックス近似が得られる.
In[13]:=
Click for copyable input
Out[13]=
以下は相対誤差のプロットである.
In[14]:=
Click for copyable input
Out[14]=

MiniMaxApproximationが役に立たない場合もある.そのような場合には問題となった可能性のある箇所と問題を回避するための可能性についてのメッセージが表示される.例えば次のミニマックス近似を求めようとすると,MiniMaxApproximationは,相対誤差が符号間で振動し,極値に 回達するような有理ミニマックス近似を探す.時にはさらに多くの極値に達することもある.この問題に対処するより強力なアルゴリズムを設計することも可能であろう.しかし,現実にはミニマックス近似を異なる次数で求める方がより簡単なことが多い.

この近似を計算しようとすると警告メッセージが出される.単一の誤差ではなく,最初の部分の横軸に対応する誤差のリストになっている点に注目する.
In[15]:=
Click for copyable input
Out[15]=
これは近似を抽出し因子を約分する.
In[16]:=
Click for copyable input
Out[16]=
相対誤差には予想された6つではなく7つの極値がある.
In[17]:=
Click for copyable input
Out[17]=
希望の近似の程度を変更すると問題は解決する.
In[18]:=
Click for copyable input
Out[18]=

この他に,区間のどこかで分母で最初の有理補間がゼロになることがある.そのような場合は,異なる次数のミニマックス近似を求めるのが最も簡単である.しかし,このようなアプローチでは問題が解けないこともしばしばある.

このような場合は,より短い区間で有効な近似を使って始めるとうまくいくことが多い.ミニマックス近似は区間が長くなるにつれ,変化していくのが一般的なので,短い区間の近似をもう少しだけ長い区間における近似の始点として用いることができる.区間を徐々に伸ばしていくことで,結果的に所望の区間で有効なミニマックス近似が得られるかもしれない.

MiniMaxApproximationには,ユーザがその動きをコントロールできる複数のオプションと失敗の診断に役立つ2つのオプションがある.

オプション名
デフォルト値
WorkingPrecisionMachinePrecision使用する精度桁数
Bias0自動的に選ぶ補間点のバイアス
Brake{5,5}反復アルゴリズムに適用するブレーキ
MaxIterations20ブレーキをかけた後の最大反復回数
DerivativesAutomatic導関数に用いる関数の指定
PrintFlagFalse反復の各段階で相対誤差に関する情報を出力するか否か
PlotFlagFalse反復の各段階で相対誤差をプロットするか否か

ミニマックス近似のオプション

MiniMaxApproximationは,まず指定の関数の有理補間を求め,次に誤差が等振動となるまで係数を摂動させる.オプションのBiasは,RationalInterpolationで用いられているのと全く同じように最初の有理近似を制御するのに使われる.

MiniMaxApproximationは進むにつれて,相対誤差の極値を見付け,それらの極値の大きさが等しくなるように係数を摂動させる.ある反復から次の反復で,極値があまり移動しなければ,前回の位置を新たな位置を求める際の始点として使うことができる.反復間での移動が大きすぎる場合は,MiniMaxApproximationは追跡できなくなってしまう.MiniMaxApproximationが追跡不可能になったと分かるのは,極値の符号が変化しない,2つの極値が同じである,あるいは極値の横軸が分類順でないという場合である.

MiniMaxApproximationが追跡不可能にならないようにするためには,オプションのBrakeを設定するとよい.Brakeは,ある反復から次の反復への変化にブレーキをかけるメカニズムとして作用する.しかし,この効果は徐々に失われ,最終的には全くなくなる.アルゴリズムがほぼ収束すれば,変化はごく微小となるので,ブレーキの必要はなくなる.

Brakeの値は2つの正の整数のリストでなければならない.最初の整数がいくつの反復にブレーキ効果を及ぼすかを指定し,2番目の整数が最初の反復にどれくらいの強さのブレーキをかけるかを指定する.Brakeは高次のミニマックス近似で重要度が増す.この場合は極値の横軸間が狭まるからである.

反復スキームを使用するためには,近似する関数の最初の2つの導関数がMiniMaxApproximationに分かっていなければならない.Mathematica が導関数を解析的に計算できなければ,それを明示的に指定する必要がある.これに関連して,導関数は求めたがその計算が非常に大変で,大部分が重複しているという場合がある.例えば, のミニマックス近似を求める際に,Mathematica を評価して関数の値を求め, を評価してを第1導関数の値を求め,さらに を評価して第2導関数の値を求める必要がある.ユーザが x の各値についてこれら3つの値のリストを返す関数を指定する方がこれよりはるかに簡単である.オプションDerivativesの目的はここにある.

このオプションを使う際に留意しなければならない点が2つある.まず,x が実際に数になるまで関数を評価してはならない.さもないとこのオプションを使う意味がなくなってしまう.

次は,x が数でない場合にが評価されるのを防ぐ.
In[19]:=
Click for copyable input
MiniMaxApproximationの定義に沿って作用する.
In[20]:=
Click for copyable input
Out[20]=

無限に反復するのを防ぐために,MiniMaxApproximationにはMaxIterationsというオプションが設けられている.ブレーキの後の反復回数がMaxIterationsに達しても収束しなければ,その段階での近似とともに警告メッセージが出される.収束に時間がかかるということだけが問題であるのなら,MiniMaxApproximationの第2引数として返された近似を挿入して,その時点の近似から再び反復を開始することができる.新たな反復は今までとは異なるオプションで始めるとよいかもしれない.

質の悪い近似の例を取るために,MaxIterationsの値を小さくし,バイアスは大きくして,ブレーキを付けずにおく.
In[21]:=
Click for copyable input
Out[21]=
前回の近似の結果を第2引数として挿入し,新たな反復の始点として使う.
In[22]:=
Click for copyable input
Out[22]=

PrintFlagPlotFlagは失敗の理由を診断するのに使えるMiniMaxApproximationの便利なオプションである.この2つのオプションのどちらを設定しても計算には影響しない.PlotFlagTrueに設定すると,各反復有理近似の相対誤差のプロットが生成される.これらのプロットが任意の2つの反復間で劇的に変化するようなら,ブレーキを大きくするとよいだろう.PrintFlagTrueに設定すると,2つのことが起る.まず,極値が初めて探索されているときに,横座標に対する近似の変化のリストが出力されるということである.このリストに含まれる数は,一旦ある程度小さくなると劇的に減少する.次は,初めて極値が見付かった後で,相対誤差の極値の横座標とそれぞれの横座標における相対誤差の値からなる順序対のリストが出力されるということである.

GeneralRationalInterpolation[{fx,fy},{t,m,k},x,{t1,t2,...,tn+k+1}]
グラフが t の関数としてパラメトリックに与えられる関数 x に次数の有理補間を与える
GeneralRationalInterpolation[{fx,fy},{t,m,k},x,{t,tmin,tmax}]
を自動的に選択して有理補間を与える

一般的な有理補間

RationalInterpolationの代りにGeneralRationalInterpolationを使った方がよい近似問題がある.例えば,FindRootを使ってしか評価することができない関数の逆関数の有理補間関数を求めたいのであれば,GeneralRationalInterpolationを使った方がよい.このような場合には,RationalInterpolationを使うと評価が非常に遅くなる.

GeneralRationalInterpolationを使うと近似する関数をパラメトリックに与えることができるため,より一般的な近似問題を解くことができる.例えば,関数 のグラフは単位円の上半分だけに当たる.これは であるとしてパラメトリックに表すこともできる.このように,関数を{Cos[t], Sin[t]}と指定することによりの近似を計算することができる.

RationalInterpolationの関数をとして指定するとき,一般的に式 t の関数である.補間される関数は の区間でによってパラメトリックにグラフが与えられるような関数である.

独立変数には常にシンボルを指定しなければならないという点に注意する.パラメータ変数を独立変数として使うのは間違っている.GeneralRationalInterpolationRationalInterpolationと同じオプションを取る.

これはグラフが円の上半分になる関数の近似を与える.
In[23]:=
Click for copyable input
Out[23]=
終点付近で誤差が極めて大きくなる.
In[24]:=
Click for copyable input
Out[24]=
補間点を明示的に指定する必要はない.
In[25]:=
Click for copyable input
Out[25]=
RationalInterpolationの場合と同じように,補間点を明示的に与えない方が誤差が小さくなることが多い.自動的に選ばれる点の方がよい分布となることが多い.
In[26]:=
Click for copyable input
Out[26]=
GeneralMiniMaxApproximation[{fx,fy},{t,{tmin,tmax},m,k},x]
グラフが t の関数としてパラメトリックに与えられる x の関数の次数の有理近似を与える
GeneralMiniMaxApproximation[{fx,fy},approx,{t,{tmin,tmax},m,k},x]
区間アルゴリズムを approx で始めてミニマックス近似を与える
GeneralMiniMaxApproximation[{fx,fy,g},{t,{tmin,tmax},m,k},x]
因数 を使って誤差を計算する有理近似を与える

一般的ミニマックス近似

近似される関数はGeneralRationalInterpolationで指定されるのと同様にGeneralMiniMaxApproximationで指定される.GeneralMiniMaxApproximationのオプションはMiniMaxApproximationのオプションと同じである.

これはExpの逆関数の次のミニマックス近似を与える.
In[27]:=
Click for copyable input
Out[27]=
Expの逆関数を評価する簡単な方法があるので,この問題にMiniMaxApproximationを使うことも可能である.解における唯一の違いは極値の横座標である.
In[28]:=
Click for copyable input
Out[28]=
次はミニマックス近似を抽出する.
In[29]:=
Click for copyable input
次は区間 における近似の相対誤差のプロットである.
In[30]:=
Click for copyable input
Out[30]=

GeneralMiniMaxApproximationでオプションDerivativesを使うのであれば,パラメトリックに定義された関数の両部分に対する導関数のリストを指定しなければならない.

これはがシンボルに作用するのを防ぐ.
In[31]:=
Click for copyable input
再び同じ結果を得た.
In[32]:=
Click for copyable input
Out[32]=

GeneralMiniMaxApproximationが便利である別の例として,近似した後,デフォルトで用いられる相対誤差ではなく絶対誤差を使ってフィットのよさを見たい場合がある.誤差の測定基準を変えたければパラメトリックに定義された関数の第3引数でこれを指定することができる.関数が として与えられ,に関連付ける関数の有理ミニマックス近似であるなら,最小にされた誤差はである.が指定されていなければ,と同じものとされ,最小化されるのは最大相対誤差である.絶対誤差を最小にしたければ,にデフォルトの定数1を使えばよい.

これは最大の絶対誤差を最小にしてExpの逆関数の有理ミニマックス近似を与える.
In[33]:=
Click for copyable input
Out[33]=
近似を抽出する.
In[34]:=
Click for copyable input
これで「絶対」誤差は同じ振動になる.
In[35]:=
Click for copyable input
Out[35]=

誤差が受け入れ難いほど大きい近似となった場合は,次のようなことをしてみるとよい.分子および/または分母の次数を上げると,ミニマックス近似の誤差は通常減少し,増大することはない.分子から分母へ,またはその逆方向へ次数を移動させるだけでも誤差の大きさに影響が出る.近似が有効となる区間が非常に長い場合は,関数の漸近挙動に注目し,適切な漸近挙動を行うように分子と分母の次数を選ぶとよいだろう.例えば,大きな に対するの近似を行うのであれば,分子と分母の次数は同じでなければならない.大きな に対して関数はほぼ一定の値 を持つからである.

誤差を小さくするこの他の方法に,近似が有効となる区間を短くすることが考えられる.長い区間に対しては,1つの高次近似を行うよりもいくつかの低次の近似を行った方がよいことが多い.余分に必要となる記憶領域はそれ程大きくないし,簡単な近似の方が速く済む.

最適化されたミニマックス近似を得るためには,最終的な近似の数値的挙動を考慮しなければならない.近似に現れる変数が原点付近の値を持つように関数を定義するとよい.つまり, 区間における のミニマックス近似を求める代りに,区間における のミニマックス近似を求めた方がよいのである.

減算の桁落ちに関連した潜在的な問題をすべて回避したい場合は,近似を求めた後で変数の桁移動を行ってもよいだろう.例えば,区間における についての次数の有理ミニマックス近似は,分子に正の係数を,分母に交代符号の係数を持つ.近似における相対誤差はわずか程度であるが,指定された形式の近似を行うと相対誤差はもっと大きくなる可能性がある.それは,終点近くの桁落ちにより数ビットなくなる可能性があるからである.分子では付近で桁落ちが起り,分母では付近で起る.分子の で置換し,分母の で置換すると,分子分母双方の係数すべてが正になり, の値はの間になる.簡約は行われない.もちろん,係数自身が要求される精度を満たすように,これらすべてがより高精度で行わなければならない.

関数の零点とその漸近挙動を考慮することも大切である.結果として得られるミニマックス近似の単純さは,関数を自明に定数関数のように見せられる程度によって大きく影響される.例として区間Gamma[1/2, x, Infinity]のミニマックス近似を求めてみよう.このとき,関数Gamma[1/2, x, Infinity]Exp[x](x+4/7)を考慮することができる(AbramowitzとStegunによる「Handbook of Mathematical Functions」の6.5.31を参照のこと.4/7は経験的に選ばれている).この関数は問題としている区間でほんの数パーセントしか変化しない.この関数のミニマックス近似を求める方がもとの関数のそれを求めるよりはるかに簡単である.

最大相対誤差を最小にしたいのなら,そして問題としている区間で関数がゼロになることがあるなら,まずこのゼロを分割しなければならない.区間中に特異点が存在する場合は,ゼロを掛けるか引き算するかしてこれも除去しなければならない.

補間を使って根を見付ける

関数の根を求めるのには関数FindRoot便利である.これはかなり強力な関数で,十分に近い初期値が与えられ根が1つなら(根が複数でもオプションが適切に設定されているなら),ほとんどいつでも根を求めることができる.この強固さを持つために,FindRootは解を得るまでの関数の評価回数にはあまりこだわらないという妥協をしている.しかし,関数は大変うまくできているが,評価が,特に高精度を求める場合に,極めて高くつくということもある.そのような場合にはInterpolateRootの方がより効率的であろう.

InterpolateRootは関数の前回の評価結果,例えばを見て,点を通る補間多項式を構築する.アルゴリズムは,この補間多項式を0で評価することで関数 の根の次の近似を求める.それまでのすべてのデータを使うことが最善の方策ではないことが分かる.データ点を追加することで収束の割合は高くなるが,この割合は二次方程式よりも高くなることは決してない.しかも,使用するデータ点が増えれば増える程,アルゴリズムの強固さは失われる.InterpolateRootは前4つのデータしか使わない.それ以上使ってもあまり意味がないからである.

InterpolateRoot[f,{x,a,b}]初期値 ab の近くで関数 f の根を求める
InterpolateRoot[eqn,{x,a,b}]初期値 ab の近くで方程式 eqn の根を求める

InterpolateRootを使って根を求める

オプション名
デフォルト値
AccuracyGoalAutomatic求める根の目標確度
MaxIterations15諦めるまでの関数評価の最大回数
WorkingPrecision40算術計算で使用する最大精度
ShowProgressFalseアルゴリズムの進行に従って,中間結果やその他の情報を表示するかどうか

InterpolateRootのオプション

AccuracyGoalAutomaticとするとAccuracyGoalWorkingPrecisionより20桁小さくなる.InterpolateRootで使われるAccuracyGoalFindRootで使われるAccuracyGoalとは異なるので注意のこと.FindRootは連立方程式に使える,はるかに一般的な関数である.根自体の値の確度を調整しようとするのは難しすぎる.FindRootは関数の値が満足する程度に小さくなったところで停止するに過ぎない.InterpolateRootはずっと特化されており,単根において単変数の単関数にのみ作用し,関数が大変うまく動作していると仮定する.このような場合には,根自体の値の確度を調整するのは極めて簡単である.

評価数を決定するためにカウンタを n として,根を求めたい関数を設定する.
In[2]:=
Click for copyable input
Log[2]を求めるのにFindRootを用いる.
In[3]:=
Click for copyable input
Out[3]=
次は結果を求めるための関数の評価回数である.
In[4]:=
Click for copyable input
Out[4]=
InterpolateRootを使うと同じ結果を得るための評価回数が少なくなる.
In[5]:=
Click for copyable input
Out[5]=
近似がどのように根に収束するかを見ることができる.
In[6]:=
Click for copyable input
Out[6]=

被積分関数の補間関数を使った積分の近似

Mathematica 関数のNIntegrateは被積分関数が少なくとも数階まで滑らかであると仮定するアルゴリズムを使用する.InterpolatingFunctionオブジェクトは通常この仮定を満足しない.これは連続であるが区分的にだけ滑らかなのである.NIntegrateが使用するアルゴリズムはInterpolatingFunctionオブジェクト,それも特に数次元のものに適用されると,非常にゆっくりと収束する.NIntegrateでは,積分領域がいくつかの部分に分けられ,各部分で積分が評価されてもよい.領域の部分が,InterpolatingFunctionが滑らかである区域にある部分に対応するならば,NIntegrateは格段に速く収束する.NIntegrateInterpolatingFunctionは自動的に積分領域を分割する.

NIntegrateInterpolatingFunction[expr,{x,xmin,xmax}]
被積分関数にInterpolatingFunctionオブジェクトを含む積分の数値近似を求める
NIntegrateInterpolatingFunction[expr,{x,xmin,xmax},{y,ymin,ymax},...]
被積分関数にInterpolatingFunctionオブジェクトを含む多次元積分の数値近似を求める

被積分関数にInterpolatingFunctionオブジェクトを含む積分の数値近似

NIntegrateInterpolatingFunctionは関数NIntegrateを使うが,領域をInterpolatingFunctionオブジェクトが滑らかとなるような部分に分割する.

二次元の振動関数を近似するInterpolatingFunctionオブジェクトを生成する.
In[2]:=
Click for copyable input
Out[2]=
このリストは,積分の評価に使われる時間と積分の結果を与える.
In[3]:=
Click for copyable input
Out[3]=
NIntegrateはほぼ全く同じ結果を出すが,領域が適切に分割されていなければ収束がよくないため,ずっと長い時間がかかる.
In[4]:=
Click for copyable input
Out[4]=

InterpolatingFunctionオブジェクトの積分だけ(そのオブジェクトの関数ではなく)が必要な場合は,Integrateを使った方がよい.これはInterpolatingFunctionオブジェクトで使われた多項式近似に対して厳密な結果を与えるからである.

Order Starプロット

微分方程式を解くにあたって,数値的メソッドの安定性を解析するのは近似次数を求めるよりも難しい.このことは,安定性の定義にはさまざまな方法があるという事実に反映されている.このパッケージは数値的方法における安定領域の決定を補助するものである.

安定領域は近似解で誤差が伝播する割合を反映するので,重要である.数値解析に絶対誤差と相対誤差があるように,安定性にも絶対尺度と相対尺度がある.このパッケージはメソッドの相対的安定性を検証するのに便利なorder starを描画する.さらに,比較関数が一様に1になるように指定することで,絶対安定領域を描画することもできる.

問題に与えられた数値メソッドは,近似理論の枠組みという視点から捉えることができる.こうなると,目標は,この近似が解と比較してどれほどうまく行われるかを調べることになる.解が分かっていれば数値近似に頼る必要はないので,これは一種の逆説である.しかし,構築しようとしているのは任意のクラスのあらゆる問題に適用できる枠組みである.問題に対する解析的な解を見付けるのは困難なので,線形系に適用されたときに数値近似がどれくらいうまく動作するかを見るのが一般的である.例えば,常微分方程式の領域では,以下のような一般に非線形な f についての連立方程式

の解に興味を持つかもしれない.この問題は,以下のように解くことが可能なスカラー線形問題に置き換えるのが一般的である.

ここで は複素定数である.また,方程式を独自に決定できるように初期条件が定めてある.これで,安定性解析は単純化された微分方程式 (1)に適用されたときに,数値近似がどのように動作するかを見ることになった.方程式 (1)はしばしばスカラー線形テスト問題,あるいはダールキスト(Dahlquist)方程式と呼ばれるものである.

以下の説明ではこのパッケージの使い方に焦点を当てる.説明の中心は近似の動作になるが,根底には近似のもととなる,他の問題に適用したときに現れる数値メソッドがあることは念頭に置いて欲しい.この対応,安定性解析,order starの理論についての詳細は,セクションの末尾にある参照文献をご覧いただきたい.

OrderStarPlot[r,f]関数 rf について,となる領域を表すorder starを描画する
OrderStarPlot[r,f,z]となる複素平面 z での領域を描画する.ここで rfz の関数である
OrderStarPlot[r,f,OrderStarKind->Second]
Re(r-f)<0となる領域を表すorder starを描画する

order starの描画

パデ(Padé)近似は,すべてのパラメータが何らかの局所的な展開点で次数を最大にするように選ばれた有理多項近似である.ルンゲクッタ(Runge-Kutta)法のような数値メソッドは指数関数のパデ近似と関連がある.

次でのパデ近似が構築される.このための関数は自動的にロードされる.この近似は前進オイラー法に対応する.
In[2]:=
Click for copyable input
Out[2]=
これは前進オイラー法の相対的安定領域,すなわち第1種order starである.近似の極点がハイライトされている.
In[3]:=
Click for copyable input
Out[3]=
これは1との相対比較で得られる前進オイラー法の絶対安定領域である.
In[4]:=
Click for copyable input
Out[4]=
OrderStarInterpolationrf が同値となる点を表示するかどうかを指定する
OrderStarKindのどちらのorder starを描画するかを指定する
OrderStarLegend使用されるさまざまのシンボルを含むプロットの凡例を含むかどうかを指定する
OrderStarPoles近似および/あるいは関数の極点を明示するかどうかを指定する
OrderStarZeros近似および/あるいは関数の零点を明示するかどうかを指定する
OrderStarSymbolSize零点,極点,補間点で使用するシンボルの大きさを指定する
OrderStarSymbolThickness零点,極点,補間点で使用するシンボルの線の太さを指定する

OrderStarPlot独特のオプション

特定の機能を描画しようとしていて Mathematica がその機能を見付けられないようなときは,問題に関するより詳しい情報を伴ったメッセージOrderStarPlot::solsが出る.Solveもまた,いつ逆関数が使用されたかというようなメッセージを出す.

OrderStarPlotは独立変数が何であるかを決定するためにヒューリスティックを用いる.複雑な式の場合は使用する変数を明示的に与えて所用時間を短縮することができる.関数は数値的には評価しないので,変数の選択に曖昧な点があると,入力が評価されずに返され,適切な警告メッセージが出される.

次は使用する変数を明示し,となる点をハイライトする.関係が非代数的である場合,これは一般に可能性でないことがある.
In[5]:=
Click for copyable input
Out[5]=

TrueFalseに加え,オプションのOrderStarInterpolationOrderStarLegendOrderStarPolesOrderStarZerosも,自動的に求まらない点を指定するために座標対のリストを取ることができる.スケールした座標を与えてプロットの凡例の大きさを変えるのに加え,使用フォントのサイズやスタイルと言った凡例の情報も指定することができる.

凡例の位置はRectangleにおけるのと同じシンタックスを使ってスケールした座標で与える.フォントのスタイルと大きさも指定され,零点と極点を表すシンボルのサイズは大きくされている.
In[6]:=
Click for copyable input
Out[6]=

OrderStarPlot独特の多くのオプションに加え,ContourPlotに渡され,プロットを作成するのに使われるオプションもいくつかある.

AspectRatioプロットの縦横比を指定する
Axes座標軸を描画するかどうかを指定する
AxesOrigin座標軸の交点を指定する
ClippingStyle曲線または曲面がプロット範囲を超えたときに,何のスタイルを描画するかを指定する
ColorFunction安定領域と非安定領域を彩色するための関数を指定する
FrameTicksフレームの周囲に目盛を置くか,もし置くならどこに置くかを指定する
MaxRecursionいくつの再帰的細分化を実行することができるかを指定する
PlotPoints等高線プロットの作成で使用するサンプル点の数を指定する
PlotRangeプロットで使用するプロット領域を指定する
Ticks軸に沿って目盛を付けるかどうか,付けるとすればどこに付けるかを指定する

多くのグラフィックス関数に共通のオプション

重要な点はorder starが仮想軸と交差するかどうかである.また,実軸について対称に図示してみるのも面白いだろう.このような比較を簡単にするために,OrderStarPlotは基点を通る軸を描画するグラフィックスオプションを使っている.

OrderStarPlotContourPlotMaxRecursionオプションおよびPlotPointsオプションを利用して,細部の問題を解決する.

デフォルトのプロット領域はorder starの基本機能で決定される.しかし,このデフォルト値は標準のGraphicsオプションでオーバーライドすることができる.

これで関数を定義する.
In[10]:=
Click for copyable input
次は関数のパデ近似を構築する.
In[11]:=
Click for copyable input
標準オプションを使ってデフォルトのプロット領域とプロット密度を変更することができる.しかし,そうすることでプロットがギザギザになる恐れがある.
In[16]:=
Click for copyable input
Out[16]=
MaxRecursionPlotPointsを大きくすることで,この問題は解決できる.
In[18]:=
Click for copyable input
Out[18]=

第2種order starは線形多段階メソッドの安定性の研究に有益である.

第2種のorder starをプロットすることもできる.
In[19]:=
Click for copyable input
Out[19]=

Order starはA安定性のような多くの重要な機能を一目で決定する手段を与える.さらに,関数との相対比較を考慮して,数値スキームの確度の次数を安定性の領域に符号化する.

の有限の零点や極点はないので,以下は近似の零点と極点だけである.この近似は左半分の平面に極点を持たないので,この近似に対応する数値メソッドはA安定性である.さらに,原点に隣接する領域の数により,近似次数が(隣接領域より1少ない)5であることが分かる.
In[12]:=
Click for copyable input
Out[12]=

数値メソッドの安定性理論の説明は[1][2][3]にある.このパッケージの応用方法の詳細は[4]に記載されている.

New to Mathematica? Find your learning path »
Have a question? Ask support »