【倒立振子part2.5】maxiamaのコード

色々調べてもコードが入力できないなぁと思っていたけど、原因が分かったので書いておく。
はてなブログには編集モードが「見たままモード」「はてな記方」「Markdown記法」の3種類があるっぽい。
はじめは「見たままモード」になっており、このモードだと一回記事を書いてしまうと後から変更できない+コードは張れないことが判明した。

基本設定で「はてな記法」にしたのでコードは無事に張れそうなので下記で使用したmaxiamaのコードを保存しておく。
marurunru.hatenablog.com

基本的には難しいことはせずに式を定義してラグランジュの運動方程式のみをmaximaに解いてもらうようなコードのみを書いた。

depends(θp,t);
depends(θr,t);
/*振子*/
/*重心座標*/
xp:lp*sin(θp)/2;
yp:lp*cos(θp)/2;
/*重心速度*/
dxp:diff(xp,t);
dyp:diff(yp,t);
/*エネルギー*/
Kpt:ratsimp((mp*dxp^2)/2+(mp*dyp^2)/2);
Up:mp*g*lp*cos(θp);
Kpw:(Jp*(diff(θp,t))^2)/2;

/*ロータ*/
/*重心座標*/
xr:lr*sin(θp)/2;
yr:lr*cos(θp)/2;
/*重心速度*/
dxr:diff(xr,t);
dyr:diff(yr,t);
/*エネルギー*/
Krt:ratsimp((mr*dxr^2)/2+(mr*dyr^2)/2);
Ur:mr*g*lr*cos(θp);
Krw:(Jr*diff(θp+θr,t)^2)/2;
/*sin^2+cos^2=1で式を整理*/
Kpt:mp*lp^2*(diff(θp,t))^2/8;
Krt:mr*lr^2*(diff(θp,t))^2/8;

/*ラグラシアン*/
L:(expand((Kpt+Kpw+Krt+Krw)-(Up+Ur)));
/*損失エネルギー*/
D:(up*diff(θp,t)^2+ur*diff(θr,t)^2)/2;

/*ラグランジュの運動方程式*/
tp : ratsimp(expand(diff(diff(L ,diff(θp,t)),t)-diff(L,θp)+diff(D,diff(θp,t))),diff(θp,t,2),diff(θr,t,2),diff(θp,t),diff(θr,t),g);
tr :  ratsimp(expand(diff(diff(L ,diff(θr,t)),t)-diff(L,θr)+diff(D,diff(θr,t))),diff(θp,t,2),diff(θr,t,2),diff(θp,t),diff(θr,t),g);

maximaで使用している関数で忘れた箇所があれば下記のまとめサイトで確認できる。
ja.wikibooks.org