MATLABで少し遊んでみた

この記事は KCS Advent Calendar 2023 の15日目の記事です。

14日目←　→16日目

自己紹介

　現在学部一年のストーンチョコです。石というよりかはチョコの人です(?)。初めての作った記事なので細かい所はお許しください(o*。_。)o

MATLABとは

　MATLAB は独自の言語（MATLAB 言語）を持つエンジニアや科学者向けに設計されたプログラミングプラットフォームです。特にベクトル演算を中心とした数値計算がはやく、深層学習や制御システムなど学術界と産業界の様々なアプリケーションに使われているようです。

MATLABを用意しよう

　そんな MATLAB には Campus-Wide License というすべての学生、教職員、および研究者が、キャンパスの内外で、MATLABとそのツールボックスを制限なく使用することができるようなライセンスが用意されており、日本では慶應義塾大学を含む100校以上の大学がこれを採用しています。塾生は次のリンクから手順に沿ってダウンロードすることができます。みなさんもぜひ使ってみてください。LinuxでのC++の環境構築よりはとても簡単。

家族、親族と共有するPCは非対象であるなどのルールがあるので、利用条件や注意事項をしっかり読んでからインストールするようにしましょう。

さっそく遊んでみる

　最近いろいろと忙しくあまり時間をかけずに作れるものはないかな～と考えて実装しようと思ったのが二重振り子です。

二重振り子とはなんぞ

　一応二重振り子について簡単に説明しておきますと、二重振り子とは振り子にまた別の振り子を連結させたものです。腕に糸を用いた振り子もありますが、この記事では剛体を用いて考えることとします（びょいびょいさせないってこと）。二重振り子は初期値（どれくらいの高さから球を落とすか等）をほんの少しでも変えると数秒後になって軌道が全く異なるものとなってしまい、そのためよくカオス理論の入門書とかに出てきます（バタフライ効果っていう用語自分の好きな理系用語TOP3に入ってるかも）。

実装

　それでは実装に入っていきます。丁度よくMATLABの公式のサイトで二重振り子運動の実装の仕方が紹介されていたのでそれを参考にしながら自分で次のコードを書きました。

　初期値が互いに少しだけ異なる（連結させた振り子の球の初期位置を１度(゜)だけずらす）二重振り子を２つ作成し、重力に任せて運動させたときの下の球の軌道をそれぞれ可視化させて比較するGIF画像を出力するものです。

double_pendulu.m

% 手順１：二重振子の質量の変位、速度、加速度の定義

syms theta_1(t) theta_2(t) L_1 L_2 m_1 m_2 g

x_1 = L_1*sin(theta_1);
y_1 = -L_1*cos(theta_1);
x_2 = x_1 + L_2*sin(theta_2);
y_2 = y_1 - L_2*cos(theta_2);

vx_1 = diff(x_1);
vy_1 = diff(y_1);
vx_2 = diff(x_2);
vy_2 = diff(y_2);

ax_1 = diff(vx_1);
ay_1 = diff(vy_1);
ax_2 = diff(vx_2);
ay_2 = diff(vy_2);


% 手順２：運動方程式の定義

syms T_1 T_2

eqx_1 = m_1*ax_1(t) == -T_1*sin(theta_1(t)) + T_2*sin(theta_2(t));
eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g; 

eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t));
eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g;


% 手順３：力の評価とシステム方程式の簡約

Tension = solve([eqx_1 eqy_1], [T_1 T_2]);   % T_1, T_2 について求める
eqRed_1 = subs(eqx_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);   % T_1, T_2 を eqx_2 に代入
eqRed_2 = subs(eqy_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);


% 手順４：システム方程式の求解

L_1 = 1;
L_2 = 1.2;
m_1 = 1.5;
m_2 = 2;
g = 9.8;

eqn_1 = subs(eqRed_1)
eqn_2 = subs(eqRed_2)
[V,S] = odeToVectorField(eqn_1,eqn_2);   % 一回微分方程式に変換

M = matlabFunction(V,'vars',{'t', 'Y'});
% 初期値の設定
initCond1 = [pi/2 0 pi/2 0]
initCond2 = [pi/2-pi/180 0 pi/2 0]

sols1 = ode45(M,[0 10],initCond1);
sols2 = ode45(M,[0 10],initCond2);

% 手順５：振動する二重振子のアニメーションの作成

x_1 = @(t) L_1*sin(deval(sols1,t,3));
y_1 = @(t) -L_1*cos(deval(sols1,t,3));
x_2 = @(t) L_1*sin(deval(sols1,t,3))+L_2*sin(deval(sols1,t,1));
y_2 = @(t) -L_1*cos(deval(sols1,t,3))-L_2*cos(deval(sols1,t,1));

x_3 = @(t) L_1*sin(deval(sols2,t,3));
y_3 = @(t) -L_1*cos(deval(sols2,t,3));
x_4 = @(t) L_1*sin(deval(sols2,t,3))+L_2*sin(deval(sols2,t,1));
y_4 = @(t) -L_1*cos(deval(sols2,t,3))-L_2*cos(deval(sols2,t,1));

fanimator(@(t) plot(x_1(t),y_1(t),'ro','MarkerSize',m_1*10,'MarkerFaceColor','r'));
fanimator(@(t) plot(x_3(t),y_3(t),'ro','MarkerSize',m_1*10,'MarkerFaceColor','r'));
axis equal;

hold on;
fanimator(@(t) plot([0 x_1(t)],[0 y_1(t)],'r-'));
fanimator(@(t) plot(x_2(t),y_2(t),'go','MarkerSize',m_2*10,'MarkerFaceColor','g'));
fanimator(@(t) plot([x_1(t) x_2(t)],[y_1(t) y_2(t)],'g-'));
fanimator(@(t) plot(x_2(0:0.01:t), y_2(0:0.01:t), 'k', 'Color', [0, 0.4, 0.6]));   % 軌道をプロット(青)

fanimator(@(t) plot([0 x_3(t)],[0 y_3(t)],'r-'));
fanimator(@(t) plot(x_4(t),y_4(t),'go','MarkerSize',m_2*10,'MarkerFaceColor','g'));
fanimator(@(t) plot([x_3(t) x_4(t)],[y_3(t) y_4(t)],'g-'));
fanimator(@(t) plot(x_4(0:0.01:t), y_4(0:0.01:t), 'k', 'Color', [0.6, 0, 0]));   % 軌道をプロット(赤)

fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)));
hold off;

writeAnimation(['double_pendulum1.gif'])   % 結果を出力

symsを用いるには "Symbolic Math Toolbox" をインストールする必要があります。Toolbox というものはライブラリと似たようなものであり、便利な機能を追加することができます。これをインストールせずに syms を用いようとすると、これをインストールするように勧めてくるエラーが出るのでそこで与えられたリンクから飛んでインストールすることもできます。

　初期位置を床に水平にした方の二重振り子の軌道を青色、若干角度をずらした方の二重振り子の軌道を赤色として、下のGIF画像が出力されました。

　時間が経つにつれて軌道がずれていっていることを確認することができます。

まとめ

　文法が普段のプログラミング言語と比べると相当異なるので慣れるには時間がかかりそう、、だけど方程式解いてくれるのがすごいし、きれいにプロットしてくれるところに感動したのでそこは MATLOVE(?)。

参考