PythonでDuffing振動子を解析してアトラクターを見る

  • このエントリーをはてなブックマークに追加

Duffing振動子は非線形の2階常微分方程式です。この方程式をある特定のパラメータ下で解析すると特徴的な模様のアトラクターを観察できます。ここではDuffing振動子の概要と、Pythonによる計算方法を紹介し、自分の手でアトラクターを描くまでを目標とします。

こんにちは。wat(@watlablog)です。ここではDuffing振動子の解析から特徴的なアトラクターを可視化してみます!

Duffing振動子の概要

Duffing振動子の方程式

Duffing振動子Duffing oscillator)は式(1)で表される2階常微分方程式です。

\[ \ddot{x} + \delta \dot{x} + \alpha x + \beta x^{3} = \gamma \cos(\omega t + \phi) \tag{1} \]

ここで\(\delta\)は\(\dot{x}\)(速度)にかかり振動子の減衰を表し、\(\alpha\)や\(\beta\)は\(x\)にかかり剛性を表すパラメータです。
特に\(\beta\)がかかっている\(x^{3}\)項は線形でないため、\(\beta\)の数値によって方程式は非線形特性を持つようになります。

方程式の右辺は強制振動における外力を示し、\(\gamma\)が外力の振幅、\(\omega\)が角振動数、\(\phi\)が位相を意味し時間\(t\)で振動します。

Duffing振動子については以下のScholarpediaのページが一番わかりやすいと思いました。今回はこのページの最初に載っている動画をPythonで作ってみます。

http://www.scholarpedia.org/article/Duffing_oscillator

アトラクターの周期変化を観察するコード

Advertisements

Duffing方程式を解析するコード

先ほどのリンク先によると、式(1)のパラメータを以下にセットして\(\frac{2\pi}{\omega}\)の周期で出力を観察することで、特徴的なアトラクターを見ることができるそうです。

  • \(\alpha\) = -1.0
  • \(\beta\) = 1.0
  • \(\delta\) = 0.2
  • \(\gamma\) = 0.3
  • \(\omega\) = 1.0

今回はSciPyodeintを使ってみました。これまで当ブログでは自作のRunge-Kutta法を使っていましたが、やはり計算速度が段違いであることに気付きました。非線形振動を計算する場合計算速度は重要です。ここは素直に巨人の肩に乗ってみましょう。

公式ドキュメントを以下に示します。

SciPy odeint:https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

それでは、基本の解析コードを紹介します。

関数f()でDuffing振動子を定義しています。メイン処理の部分でパラメータをセットしodeintを使っていますが、このコードはそれほど長くないのでコード内コメントを見て頂ければ内容がわかると思います。

コードを実行すると以下の結果を得ます。時間波形として変位と速度をプロットしていますが、変位と速度(変位の微分)で構成される位相平面Phase plane)もプロットしています。この位相平面はポアンカレ切断面Poincaré section)、またはポアンカレ写像図Poincaré map)とも呼ばれます。

ダフィング振動子の解析結果例

特徴的なプロットが確認されました!

ポアンカレ写像図のアニメーションを作成するコード

次にDuffing振動子解析時の位相\(\phi\)を少しずつずらしながら\(2\pi\)まで計算し、結果をアニメーションにしてみます。これが先ほど紹介したScholarpediaのページに載っていた動画となります。

以下の動画が実行結果として得られるgifファイルです。
\(\frac{2\pi}{\omega}\)の周期で抽出した時間波形の見た目は完全なランダムなように見えますが、位相平面を見ると全体的にはある範囲にまとまり、さらにその形状はなんとも不思議な形が変形を伴って周期的な挙動をしています。

Duffing振動子のアトラクター動画

不思議ですね…!計算結果だけではなく、これを見る僕もアトラクターの影響を受けているようです。

ジャパニーズ・アトラクターも見てみる

ジャパニーズ・アトラクターJapanese Attractor)というのは、日本人の上田睆亮氏が発見したアトラクターで、有名なローレンツ・アトラクターとほぼ同時期に発見された\(^{[1]}\)とのことです。

上田の名前をとってウエダ・アトラクターと呼ばれることもあります。このウエダ・アトラクターは非線形インダクタンスを持つ直列共振回路の回路方程式から導出されたもので、式(1)のDuffing振動子を以下のパラメータで解析すると可視化可能です。

  • \(\alpha\) = 0.0
  • \(\beta\) = 1.0
  • \(\delta\) = 0.1
  • \(\gamma\) = 12.0
  • \(\omega\) = 1.0

このアトラクターは位相により広範囲に可動するので、初期値[0, 0]では発散してしまう場合があります。そのため位相0時のPythonコードは以下のようにパラメータと共に初期値を設定しました。

静止画として可視化した結果がこちら。この結果もすごく不思議な形をしていますね。

ジャパニーズ・アトラクターの静止画

動画を作成するコードはこちら。

そして結果がこちら。ジャパニーズ・アトラクターも位相を進めていくと時計回りに回転しているような挙動をしていますが、バウンドする回数に違いがあるみたいです。これが何を意味するのかはさっぱりわかりませんが、見ていて飽きない不思議な結果ですね…。

ジャパニーズ・アトラクターの動画

まとめ

この記事はここまで。
Duffing振動子は3乗項が非線形ばねの特性を持ち工学的に面白い周波数応答関数を持つのですが、それは次の記事で扱います。

今回はDuffing振動子の魅力として、アトラクターの可視化のみを扱いました。ここで紹介した形以外もたくさん発見されているようなので、是非皆さんのお手持ちのPCでカオスを体験してみてください。

参考文献

[1]:逢沢明, カオスの縁から<複雑適応系>を探検する偏 複雑系はいつも複雑, 現代書館, 1997, pp174

アトラクターの魅力にすっかりはまりました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*