MATLABで運動方程式の導出(Symbolic Math Toolbox使用)

windows10
MATLAB 2018b (with Symbolic Math Toolbox)

今回は題材として以下の単振り子を取り上げます．簡単なので…

まずはシンボリック変数の定義を行います．

syms theta(t) L m g

theta(t)はtを従属変数とするシンボリック変数thetaを表しています．
天井と糸の付け根の部分を原点とします．
重りの位置座標を $L$ と $\theta$ で表すと，
$x(t)=L\sin\theta$
$y(t)=L\cos\theta$
と表されるので，

x =L*sin(theta)
y = L*cos(theta)

となります．
ここで，xとyには(t)がありませんが，ｔが従属変数として扱われてます．
xもx(t)も同じなようです．
運動方程式の導出には速度が必要なので，これの微分は，

vx = diff(x,t)
vy = diff(y,t)

でできます．

さて，運動エネルギーは
$T = \frac{1}{2}*m*(v_x^2+v_y^2)$
なので，

T = m*(vx^2 + vy^2)/2
T = simplify(T)

でできます．
二行目のsimplifyは数式を整形してくれる関数です．
位置エネルギーは

U = m*g*(1-y)

となります．
位置エネルギーの基準を最下点にしています．

L = T-U

となり，
運動方程式はthetaについて解くので，

eqn = functionalDerivative(Lag,theta)==0

でできます．
式整形したいときはisolate関数を使用して，

isolate(eqn,diff(theta,t,t))

とすれば良いです．isolateは'isolate(方程式，左辺にしたい変数)'と使います．
今回はeqnにたいしてthetaの二階微分を左辺に移行したかったのでこのような表記になっています（diff(theta,t,t)が二階微分）

ここまできたら一般解まで求めてしまいましょう．
ただし $\theta$ は微小という条件で．
まず，mとLは絶対に正なので，

assume(m,'positive')
assume(L,'positive')

とします．
そして， $\theta$ が微小なら $\sin(\theta) \approx \theta$ なので，

eqn = subs(eqn,sin(theta),theta)

とします．subsはeqn中のsin(theta)をthetaに置き換えています．
解はdsolve関数を用いて，

Sol = dsolve(eqn,theta(0) == 1)

とすることで解けます．
t=0でtheta=1[rad]という初期条件です（微小か？というツッコミはなしで）

運動方程式の導出～解までできました．

Tagebuch von Spargel