Next: , Previous:   [Contents][Index]

50 dynamics


50.1 Introduction to dynamics

追加パッケージdynamicsには、 離散力学系とフラクタルの様々なグラフィックス表現を生成するためのいくつかの関数と、 微分方程式系を数値的に解くための4次Runge-Kutta法の実装が含まれています。

このパッケージの関数を使うには、 最初にload("dynamics")でパッケージをロードしなければいけません。

Maxima 5.12で導入された変更

Maxima 5.12以来、現在まで、dynamicsパッケージは グラフを処理するのに関数plot2dを使います。 (juliamandelbrotを除いて) グラフィックスを生成するコマンドはplot2dのいかなるオプションも受け付けます。 オプションには、プロットスタイルや色を使ったり、 一方の、または両方の軸を対数スケールで表したり、 様々なグラフィカルインターフェイス全体に渡って変更を加えるものを含みます。 古いオプションdomain, pointsize, xcenter, xradius, ycenter, yradius, xaxislabel, yaxislabelは この新しいバージョンでは受け付けられません。

現在,すべてのプログラムは、 以前のバージョンのようにxyだけでなく、 任意の変数名を受け付けます。 2つの要求パラメータが2つのプログラムで変えられました: evolution2dは、現在、 2つの独立変数を陽に指定するリストを要求し、 orbitsの水平範囲はもはやステップサイズを要求しません; 範囲は、変数名と最小値、最大値だけを指定しなければいけません; ステップ数は、現在、オプションnticksで変えることができます。


50.2 Functions and Variables for dynamics

関数: chaosgame ([[x1, y1]...[xm, ym]], [x0, y0], b, n, ..., options, ...);

いわゆるカオスゲームを実装します: 初期点(x0, y0)がプロットされ、 m個の点[x1, y1]...[xm, ym]の1つがランダムに選択されます。 プロットされる次の点はプロットされた以前の点からランダムに選ばれた点までの線分上で、 ランダム点からその線分の長さにbを掛けた距離にあります。 手続きは n回繰り返されます。

関数: evolution (F, y0, n, ..., options, ...);

2次元グラフに、 点の水平座標が整数0, 1, 2, ..., n、 垂直座標が再帰関係

        y(n+1) = F(y(n))

で定義された数列の対応する値y(n)である n+1個の点を描きます。 ここで、初期値y(0)y0に等しいです。 Fは1変数にだけ依存する式でなければいけません。 (例の中では、yに依存していますが、他のいかなる変数を使うことができます), y0は実数でなければいけなく、nは正の整数でなければいけません。

関数: evolution2d ([F, G], [u, v], [u0, y0], n, ..., options, ...);

再帰関係を伴う二次元離散力学系によって定義された点列の中の最初のn+1点を 二次元プロットで、 表示します。 初期値u0v0を持つ

        u(n+1) = F(u(n), v(n))    v(n+1) = G(u(n), v(n))

FGは2変数uvのみに依存する式2つでなければいけません。 変数はリストの中に明示的に指名されなければいけません。

関数: ifs ([r1, ..., rm], [A1, ..., Am], [[x1, y1], ..., [xm, ym]], [x0, y0], n, ..., options, ...);

反復関数系法を実装します。 この方法は関数chaosgameで記述した方法に似ていますが、 現在点からランダムに選ばれた点にセグメントを縮める代わりに、 そのセグメントの2成分はランダムに選ばれた点に対応する2行2列行列Aiを乗算されます。

m個のアトラクティブな点の1つのランダムな選択を 重みr1,...,rmで定義された非均一な確率分布で作ることができます。 それらの重みは累積形で与えられます; 例えば、 もし確率0.2, 0.5, 0.3を持つ3点があるなら、 重みr1, r2r3は2, 7, 10とできます。

関数: julia (x, y, ...options...)

複素数(x + i y)に関するJulia集合の表現のグラフィックスファイルを生成します。 パラメータxyは実数でなければいけません。 ファイルは、XPMグラフィックスフォーマットでカレントディレクトリかユーザーディレクトリに生成されます。 プログラムは走らせるのに数秒かかるかもしれません。終了後、生成したファイル名と共にメッセージが印字されます。

Julia集合に属さない点には、異なる色が割り当てられます。 繰り返し回数に従って、その点から始まり、半径2の収束円から動く数列を取ります。 繰り返しの最大回数はオプションlevelsで設定されます; その回数の繰り返しの後、もし数列がまだ収束円内なら、 点はオプションcolorで定義された色で塗られます。

Julia集合に属さない点に使われる色すべては、同じsaturationvalueを持ちますが、 hueと(hue + huerange)の間に一様に分布する、違った色相角を持ちます。

optionsはオプションの列です。 受け付けられるオプションのリストは以下の節で与えられます。

関数: mandelbrot (options)

Mandelbrot集合の表現のグラフィックスファイルを生成します。 ファイルは、XPMグラフィックスフォーマットでカレントディレクトリかユーザーディレクトリに生成されます。 プログラムは走らせるのに数秒かかるかもしれません。終了後、生成したファイル名と共にメッセージが印字されます。

Mandelbrot集合に属さない点には、異なる色が割り当てられます。 繰り返し回数に従って、その点から始まり、半径2の収束円から動く数列を取ります。 繰り返しの最大回数はオプションlevelsで設定されます; その回数の繰り返しの後、もし数列がまだ収束円内なら、 点はオプションcolorで定義された色で塗られます。

Mandelbrot集合に属さない点に使われる色すべては、同じsaturationvalueを持ちますが、 hueと(hue + huerange)の間に一様に分布する、違った色相角を持ちます。

optionsはオプションの列です。 受け付けられるオプションのリストは以下の節で与えられます。

関数: orbits (F, y0, n1, n2, [x, x0, xf, xstep], ...options...);

パラメータxを持つ一次元離散力学系の族に関する軌道図を描画します; この種の図は一次元離散系の分岐の研究に使われます。

関数F(y)は、 関数evolutionの場合と同様に値y0で始まる数列を定義しますが、 この場合、その関数は、x0からxfまでの区間内の値を取り、 xstepの増分を持つパラメータxにも依存します。 パラメータxに使われるそれぞれの値は水平軸に示されます。 垂直軸は、 数列にn1回の時間発展させた後得られる数列y(n1+1),..., y(n1+n2+1)n2個の値を示します。

関数: rk (ODE, var, initial, domain)
関数: rk ([ODE1,...,ODEm], [v1,...,vm], [init1,...,initm], domain)

4次のRunge-Kutta法を使って、最初の形式は一階常微分方程式一つを数値的に解き、二番目の形式はそれらm個の方程式系を解きます。 varは従属変数を表します。 ODEは独立変数と従属変数にだけ依存する式でなければいけません。 そして、独立変数に関する従属変数の導関数を定義します。

独立変数はdomainで指定されます。 それは4つの要素のリストでなければいけません。 例えば:

[t, 0, 10, 0.1]

リストの最初の要素は独立変数を特定し、二番目と三番目の要素はその変数の初期値と最終値であり、 最後の要素はその区間内で使用されるべき増加分を設定します。

もしmこの方程式が解かれようとしているなら、 m個の従属変数v1, v2, ..., vmが存在しなければいけません。 それらの変数の初期値はinit1, init2, ..., initmとなります。 以前の場合と同様、domainで定義されたただ1つの独立変数が残っています。 ODE1, ..., ODEmは 独立変数に関する従属変数それぞれの導関数を定義する式です。 それらの式に現れるかもしれない変数は独立変数と任意の従属変数だけです。 従属変数と厳密に同じ順序でリストの中に導関数ODE1, ..., ODEmを与えることが重要です; 例えば、リストの三番目の要素は三番目の従属変数の導関数と解釈されます。

プログラムは 方程式を独立変数の初期値から最終値まで一定の増加分を使って積分しようとします。 もしあるステップで従属変数の1つが大きすぎる絶対値を取ったら、 積分はその点で中断されます。 結果は、なされた繰り返しの回数と同じ数の要素を持つリストです。 結果リストの中のそれぞれの要素は、それ自身m+1個の要素を持つもう一つのリストです: 独立変数の値にその点に対応する従属変数の値が続きます。

関数: staircase (F, y0, n, ...options...);

再帰関係によって定義された数列に関する階段図形を描画します。

        y(n+1) = F(y(n))

入力パラメータの解釈と許される値は 関数evolutionに関するものと同じです。 階段図形は線G(y) = yと共に関数F(y)のプロットから構成されます。 垂直区間は、 その線上の点(y0, y0)から 関数Fと交差する点まで描画されます。 水平区間はその点から線上の点(y1, y1)に届くまで描画されます。 手続きは点(yn, yn)に届くまでn回繰り返されます。

オプション

それぞれのオプションは複数の項目のリストです。 最初の項目はオプション名で、残りはオプションの引数からなります。

関数evolution, evolution2d, staircase, orbits, ifs, chaosgameが受け付けるオプションは はplot2dに関するオプションと同じです。 orbitsは、それらのオプションに加えて、 垂直方向に表される異なる点の最大数を設定する特別なオプションpixelsを受け付けます。

以下のオプションを 関数juliamandelbrotは受け付けます:

以下の数列のグファフィックス表現と階段図形: 2, cos(2), cos(cos(2)),...

(%i1) load("dynamics")$

(%i2) evolution(cos(y), 2, 11);

(%i3) staircase(cos(y), 1, 11, [y, 0, 1.2]);
./figures/dynamics1 ./figures/dynamics2

もしシステムが遅いなら、 以下の例のように繰り返し回数を減らさなければいけないでしょう。 そして、もしモニタ上のドットが小さすぎるなら、 [style,[points,0.8]]のように、 違ったスタイルを試したくなるかもしれません。

パラメータaを持つ二次写像の軌跡図。

        x(n+1) = a + x(n)^2
(%i4) orbits(x^2+a, 0, 50, 200, [a, -2, 0.25], [style, dots]);
./figures/dynamics3

x = -1.25近くのより低い分岐の回りの領域を拡大するには、以下を使ってください:

(%i5) orbits(x^2+a, 0, 100, 400, [a,-1,-1.53], [x,-1.6,-0.8],
             [nticks, 400], [style,dots]);
./figures/dynamics4

フラクタルに導かれる二次元系の発展:

(%i6) f: 0.6*x*(1+2*x)+0.8*y*(x-1)-y^2-0.9$

(%i7) g: 0.1*x*(1-6*x+4*y)+0.1*y*(1+9*y)-0.4$

(%i8) evolution2d([f,g], [x,y], [-0.5,0], 50000, [style,dots]);
./figures/dynamics5

そしてそのフラクタルの小さな領域の拡大:

(%i9) evolution2d([f,g], [x,y], [-0.5,0], 300000, [x,-0.8,-0.6],
                  [y,-0.4,-0.2], [style, dots]);
./figures/dynamics6

カオスゲームで得られるSierpinski三角形のプロット:

(%i9) chaosgame([[0, 0], [1, 0], [0.5, sqrt(3)/2]], [0.1, 0.1], 1/2,
                 30000, [style, dots]);
./figures/dynamics7

反復函数系で得られるBarnsleyのシダ:

(%i10) a1: matrix([0.85,0.04],[-0.04,0.85])$

(%i11) a2: matrix([0.2,-0.26],[0.23,0.22])$

(%i12) a3: matrix([-0.15,0.28],[0.26,0.24])$

(%i13) a4: matrix([0,0],[0,0.16])$

(%i14) p1: [0,1.6]$

(%i15) p2: [0,1.6]$

(%i16) p3: [0,0.44]$

(%i17) p4: [0,0]$

(%i18) w: [85,92,99,100]$

(%i19) ifs(w, [a1,a2,a3,a4], [p1,p2,p3,p4], [5,0], 50000, [style,dots]);
./figures/dynamics8

そして数(-0.55 + i 0.6)に関するJulia集合は以下で得られます:

julia (-0.55, 0.6, [iterations, 36], [x, -0.3, 0.2],
      [y, 0.3, 0.9], [grid, 400, 400], [color_bar_tics, 0, 6, 36])$

x方向に-0.3から0.3まで、y方向に0.3から0.9までの領域を表示します。 青から始まり黄で終わる36色が使われます。

./figures/plotting4

微分方程式

          dx/dt = t - x^2

を、初期値x(t=0) = 1で、0から8までのtの区間で0.1の増分で 数値的に解くには、以下を使ってください:

(%i20) results: rk(t-x^2,x,1,[t,0,8,0.1])$

結果はリストresultsに保存されます。

        dx/dt = 4-x^2-4*y^2     dy/dt = y^2-x^2+1

を、t=0のxとyの初期値がそれぞれ-1.25, 0.75で、0から4までの区間で、 数値的に解くには:

(%i21) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$

Next: , Previous:   [Contents][Index]