Quick start guide for Duffing Equation


Duffing方程式

例として,Duffing方程式,すなわち,


eq_img01.gif


にみられる周期解の周期倍分岐集合を求めてみましょう.

次の3つのプログラムを作成し,順番に用いることによって計算を行います.

  1. PP
  2. FIX
  3. BF

環境設定

  • まずインストールを参考に BUNKI をインストールし,環境設定を行って下さい.
  • 補足.1
    • Symbolic Math Toolbox (MATLAB 追加パッケージ)を利用出来ない方は,デモ用として以下のサンプルを用意してあります.このモデルの分岐解析に必要な解析ツール,それら設定ファイルなどが 入った Project ディレクトリをひとまとめにしたファイルです.適当な場所に展開して下さい.
    • 注意:BUNKI コマンドで展開したディレクトリを読込むと SE が起動します.ここで system create を行うと,Symbolic Math Toolbox が無い場合はエラーを返されます.Symbolic Math Toolbox が無い場合は,プログラムの作成 の項は飛ばして,PP に進んでください.
  • 補足.2
    • プログラム BUNKI を展開したディレクトリ内にサンプル集を集めた sample ディレクトリーがあります.
      /BUNKI/Sample/Duffing
      このディレクトリ下に置かれた種々のファイルは,以下の説明で使用される設定ファイルです.ファイルを読み込むことで,各フィールドに書き込む手間が省けます.用意された設定ファイルは以下の通り:
      • fileDuffing.sed : モデル方程式が記述されたファイル:SE 用の設定ファイル
      • filePd_point.ppd : 安定周期解の例:PP 用の設定ファイル
      • file2Pd_point.ppd : 安定2周期解の例:PP 用の設定ファイル
      • file4Pd_point.ppd : 安定4周期解の例:PP 用の設定ファイル
      • file8Pd_point.ppd : 安定8周期解の例:PP 用の設定ファイル
      • filechaos.ppd : カオス振動解の例:PP 用の設定ファイル
      • filePd_point.fixd : 周期倍分岐を検出:FIX 用の設定ファイル
      • file2Pd_point.fixd : 2周期解の周期倍分岐を検出:FIX 用の設定ファイル
      • file4Pd_point.fixd : 4周期解の周期倍分岐を検出:FIX 用の設定ファイル
      • filePd_point_PD.bfd : 周期解の周期倍分岐点用設定ファイル:BF 用の設定ファイル

プログラムの作成

これからの作業はサブ・ディレクトリで行うことにしましょう. ディレクトリの名前はなんでもよいです.まず MATLAB command window 内で以下を実行します.

>> BUNKI

実行後,Select BUNKI Project Directory という window が開きます.

project.png

ここで解析を行う Project を作成または選択します. この段階では,まだ何も Project を生成していないので, Name 欄に識別用の Project 名を命名して OK ボタンを押します.ここでは,Duffing 方程式の解析なので, Project 名を Duffing として話を進めます*1

project2.png

新規に Project を生成する場合,ディレクトリ(or フォルダ)を生成しても良いかを尋ねられます. Yes を選択して下さい. SE(system editor) が起動します.

startSE.png

さて,いよいよ解析プログラムの作成です*2

  • システムの記述
    • まず解析する方程式系がどういった種類の系であるかを選択します. Duffing 方程式は,2変数, 周期的外力をもった非自律系方程式 ですので,Dimension は , type は non-autonomous をチェックします.

      startSE2.png

    • Equations 欄に Duffing 方程式を記述します.その際,次のルールに従って記述して下さい.
      • 状態変数:各状態変数は x[i], i = 1,2,...,n として下さい.すなわち,x -> x[1]y -> x[2]と読み替えて式を記述してください.
      • パラメータ: パラメータであることを識別するために, $ 記号でパラメータを挟む形で記述して下さい.例えば $para$,$Vm$,$alpha$,etc,.$ 記号で挟まれたものはパラメータとして認識されます.そのパラメータ名の命名は自由です.
      • 正弦波外力: ここで利用できる周期的外力は,単一周波数をもつ正弦波に限ります.すなわち cos(t) または, sin(t) です.各方程式は,角周波数 omega で正規化されます*3

        startSE3.png

  • 文法チェック&プログラム生成
    • 以下が入力終了後の様子です.方程式の記述に矛盾が無いかをチェックするために, check ボタンを押して下さい.文法チェックが行われます.

      startSE4.png

    • 文法チェックを行った後,エラーが無ければ,いよいよプログラムを生成します.どのプログラムを生成するかを右下の check ボタンで制御します.分岐曲線を求める為には,PPFIXBF 全てにチェックを入れて下さい.Create ボタンを押して下さい.プログラム生成を開始します.

      startSE5.png

      (a): システム生成開始時に status bar が表示される.

      startSE6.png

      (b): システム生成完了後,log window に正常生成のメッセージが出力される.

PP

ここでは,パラメータと状態の初期値を読み込んで位相面図を描き,安定な周期解をつかまえます.

まず SE から PP(phase portrait) プログラムを起動しましょう. Program -> PP を選択します.

pp.png

描画 window とコントロールパネルが起動します.初回起動時には,描画 window は, (x[1],x[2])平面となっています*4.ここでは,描画 window の変更は必要ないでしょう.

pp2.png

次にコントロールパネルです.まず初期状態点を入力します*5.対応する状態(x[1],x[2])にそれぞれスタートさせたい初期点を記述します.ここでは原点(0,0)からスタートさせてみましょう!

描画の前にパラメータ値を記述することを忘れないで下さい.初回起動時には全てのパラメータに1 が入力されます.それでは,各パラメータを以下のように設定します.

parametervalue
B0.1
B00.1
c10
c31
k0.1
omega1


それでは,Startボタンを押してシミュレーションを開始しましょう.原点からスタートした解軌道は,ある周期解に収束していく様子が描かれたでしょうか?青の実線解軌道を,赤点は,時刻 2pi*周期毎にサンプリングされた点,すなわちPoincare写像点を表しています.

pp3.png

パラメータはリアルタイムに変化することができます.取り敢えずパラメータ B を変えてみましょう.パラメータ B を増加することで2周期解が発生する様子をみることができます.この時,''Clear ボタンを押して過渡応答に相当する軌道を消去することで2周期解を確認出来ます.更に B を増加すれば,4 -> 8 -> 16周期解の発生(すなわち周期倍分岐の連鎖発生)をみることができ,カオス振動発生の典型的なルートを体感できるでしょう.

pp4.png

次に,B を初期設定値 0.1 に戻してシミュレーションを再開します.この時,Plot 欄にある states flow にチェックを入れます.こうすることで,数値積分を行っている現在の状態変数値リストをリストボックス内でチェックできます.リストボックス内には,2π周期毎にサンプリングされた状態変数値が出力されます.この例では,Duffing 方程式系は Poincare 写像の安定固定点へ収束していきますので,リストボックス内には,同じ数値が流れていくのを確認できます.

pp5.png

また,この時,Clear ボタンを押すことで,描画 window 内に表示されていた過渡応答分の軌道が Clear されて定常状態に達した時の様子を見ることができます.Orbit のチェックを外して,解軌道の描画も停止させましょう.十分時間が経過した後,安定固定点へ収束します.そのため,描画 window 上では安定固定点に相当する一点だけが表示されます.

pp6.png

それでは,この周期解(Poincare 写像の固定点)の分岐を調べてみましょう.パラメータを変化させた時の周期解に対する固有値(特性乗数)を評価します.この目的には FIX ツールが利用できます.FIX を動かす為の初期値を PP から与えます.ここで重要な事は,定常状態に至った時点での状態変数値,パラメータ値を FIX に与えなければならないということです.

定常状態へ至ったかどうかは,上述した states flow スイッチを入れて,状態変数値に変化がなくなったことを確認するか,Staus ボタンを数度クリックしてみて,リストボックスに表示される情報に変化が無いことを確認する必要があります.

pp7.png

変化がなくなった事を確認した後,メニュ−バーの Tools から, Export current status を選択します.この操作により,定常状態に至った時点での状態変数値,パラメータ値,その他の設定値が保持され,次のプログラムに引き継ぐことが可能になります(FIX の項目を参照).

pp8.png

(a): export 選択.

pp9.png

(b): 設定値の export 完了.

FIX

ここでは,パラメータを変化させて周期解を追跡し,そのときの固有値を表示させることが目的です.

FIX 動作の設定値はユーザ側でそれぞれ調整する必要があります.ここでは簡単のために,そのまま動作させてみましょう.FIX プログラムは,PP プログラムで保持された値 (Export Current status) をデータを読み込んでニュートン法への初期値として実行するのが通常の方法です*6

まず PP コントロールパネルからFIX プログラムを起動しましょう.Program -> FIX を選択します.

fix.png

初回起動時には,各フィールドにはデフォルト値が入力された状態で起動します.

fix2.png

PPで求められた周期解の固定点に関する情報を初期設定値としてインポートします.メインパネル Tools -> Import initial point を選択します.

fix3.png

選択後,PP を利用して得たパラメータ値,固定点の座標値などが各項目に反映されます.

fix4.png

PP でパラメータ B を大きくすることで, 系のダイナミクスに変化が起きることが判ったので,変化させるパラメータを B とし、その増分ステップサイズを設定します.

fix5.png

Start ボタンを押して計算を開始します.計算開始と同時に,別ウィンドウが開きます.これは,Gauss 平面上のどの位置に固有値が位置しているかを示す図です.周期倍分岐の発生条件は,固有値のひとつが "-1" になるときです.パラメータの変化とともに,ウィンドウ上の固有値の一つが "-1" へと近づいていく様子を観察できます.

パラメータを連続的に変化させていき,分岐点が検出されるとプログラムは停止します*7.上の例では,固有値の一つが -1 となったことで周期倍分岐の発生を検出した様子を示しています.

fix6.png

停止した時のパラメータ値,状態変数値は分岐曲線を追跡するプログラム BF の初期点として利用できます.停止した際のパラメータ値などを BF 用に Export します.ここで Export BF point を選択すると,状態・パラメータ値,その他必要な設定値が保持されます.この保持された値を次のプログラムに引き継ぎます(BF の項目を参照).

fix7.png

(a): Export 選択.

fix8.png

(b): リストボックス内に Export した情報が表示される.

BF

ここでは,パラメータを変化させて分岐曲線を求めます.

BF 動作の設定値はユーザ側でそれぞれ調整する必要があります.ここでは簡単のために,そのまま動作させてみましょう.BF プログラムは,FIX プログラムで検出された分岐点での状態変数やパラメータ値を Newton 法の初期値として読み込んで実行するのが通常の方法です*8

まず FIX メインパネルから BF プログラムを起動しましょう. Program -> BF を選択します.

bf.png

初回起動時には,各フィールドにはデフォルト値が入力された状態で起動します.

bf2.png

FIX で求めた分岐点に関する各種の情報を初期値として Import します.メインパネル Tools -> Import BF point を選択します.

bf3.png

選択後,分岐の種類, 固有値の偏角情報などの付属情報が各項目に入力されます.(B,B0) 平面上の分岐曲線を追跡することを考え, X パラメータとして B を, Y パラメータとして B0 にチェックを入れます.また step size,end などを設定します.計算結果をデータファイルとして残す場合は, Output に出力ファイル名を書いておきます.

bf4.png

ここまで設定が終われば,後は Start ボタンを押して分岐曲線の追跡を開始します.そのとき別ウィンドウとして,Eigenvalues window と Bifurcation diagram window が起動します.Eigenvalues window は,分岐曲線を計算している際,Neimark-Sacker 分岐に対応する固有値だけではなく,それ以外の固有値の配置変化を視覚化する目的で表示されます.例えば,余次元2の分岐点が発生したことなどを視覚的に理解するために使用します.Bifurcation diagram window は,周期倍分岐曲線を追跡していく様子をリアルタイムに確認するための window です.

bf5.png

(a):固有値配置図とリストボックスの様子.

bf6.png

(b):(B,B0) 平面の分岐図.リアルタイムに計算結果を見ることができる.




*1 Project 名は自由に命名出来ます.ユーザーが識別しやすい名前にすれば良いです
*2 サンプルディレクトリの下に fileDuffing.sed という SE 用の設定ファイルがあります.このファイルを読み込む事で,あとはプロクラムの生成をするだけの状態にすることが可能です.すなわち,これから行う諸々の設定を省略することができます.ファイルの読み込み方法は,SE設定ファイルを読み込むには?を参照して下さい.
*3 自動的にパラメータ omega が付加されます
*4 この設定は,コントロールパネル Setting から変更することが出来ます
*5 サンプルディレクトリの下に filePd_point.ppd という PP 用の設定ファイルがあります.このファイルは以下の説明にある周期解をみるためのパラメータ設定や初期条件が保存された PP 設定ファイルです.あとは読み込んだ後にスタートボタンを押すだけでシミュレーションが可能です.ファイルの読み込み方法は,PPの設定ファイルを読込むには?を参照して下さい.なお,file2Pd_point.ppdfile4Pd_point.ppdfile8Pd_point.ppdfilechaos.ppd といった初期設定ファイルが併せて置いてあります.それぞれ,2周期点初期設定ファイル,4周期点初期設定ファイル,8周期点初期設定ファイル,カオス初期設定ファイルです.興味のある設定ファイルを PP で読込んでシミュレーションを楽しんで下さい.
*6 一方,パラメータ,状態変数値が記述さたファイルを読み込んで実行する方法もあります.サンプルディレクトリの下に filePd_point.fixd という FIX 用の設定ファイルがあります.このファイルは以下で説明される,周期解のパラメータ設定や初期条件が保存された FIX 用の設定ファイルです.読み込んだ後にスタートボタンを押すだけで数値計算を開始することが可能です.ファイルの読み込み方法は,FIX の設定ファイルを読込むには?を参照して下さい.なお,file2Pd_point.fixdfile4Pd_point.fixd とした2周期点用,4周期点用設定ファイルも併せて置いてあります.どの設定ファイルを読込んで FIX を動作させても,対応する周期点の周期倍分岐の発生をみることが出来ます.興味のある設定ファイルを FIX で読込んで分岐点の探索を楽しんで下さい.
*7 この動作の詳細は Setting から option パネルを動作させてコントロールすることが可能です.詳細はマニュアルまたは,個別対策マニュアルを参考にして下さい.
*8 一方,パラメータ,状態変数値が予め記述さたファイルを読み込んで実行する方法もあります.サンプルディレクトリの下に filePd_point_PD.bfd という BF 用の設定ファイルがあります.このファイルは以下で説明される固定点の周期倍分岐曲線を追跡するためのパラメータ設定や初期条件が保存された BF 設定ファイルです.読み込んだ後にスタートボタンを押すだけで分岐曲線の追跡が可能です.ファイルの読み込み方法は,BFの設定ファイルを読込むには?を参照して下さい.

Attach file: fileeq_img02.gif 318 download [Information] fileeq_img01.gif 493 download [Information] fileDuffing-sample.zip 457 download [Information] filechaos.ppd 463 download [Information] file4Pd_point.fixd 478 download [Information] file4Pd_point.ppd 452 download [Information] file2Pd_point.fixd 503 download [Information] file2Pd_point.ppd 474 download [Information] filePd_point_PD.bfd 501 download [Information] filePd_point.fixd 489 download [Information] filePd_point.ppd 455 download [Information] filebf6.png 558 download [Information] filebf2.png 541 download [Information] fileproject2.png 514 download [Information] fileDuffing.sed 462 download [Information] filefix7.png 508 download [Information] filefix3.png 496 download [Information] filestartSE5.png 501 download [Information] filefix6.png 564 download [Information] filepp8.png 498 download [Information] filefix.png 571 download [Information] filepp4.png 586 download [Information] filestartSE.png 569 download [Information] filestartSE4.png 472 download [Information] filefix2.png 530 download [Information] filepp9.png 497 download [Information] filebf3.png 560 download [Information] filepp5.png 498 download [Information] filefix8.png 487 download [Information] filestartSE6.png 501 download [Information] filestartSE2.png 515 download [Information] filefix4.png 551 download [Information] filebf5.png 544 download [Information] filepp7.png 496 download [Information] filepp3.png 496 download [Information] filepp.png 538 download [Information] fileproject.png 479 download [Information] filebf4.png 520 download [Information] file8Pd_point.ppd 479 download [Information] filebf.png 530 download [Information] filestartSE3.png 484 download [Information] filepp2.png 522 download [Information] filefix5.png 548 download [Information] filepp6.png 555 download [Information]

Front page   Edit Freeze Diff Backup Upload Copy Rename Reload   New List of pages Search Recent changes   Help   RSS of recent changes
Last-modified: 2009-07-23 (Thu) 20:18:53 (3309d)