大学受験に出る近似を含んだ単振り子問題

運動方程式と理論解

単振り子問題の運動方程式は,
$$ m\frac{d^2(l\theta)}{dt^2}=-mg\sin\theta $$
と表される. この方程式を変形すると,
$$ \frac{d^2\theta}{dt^2}+\frac{g}{l}\sin\theta = 0 $$
となる. 大学受験で単振り子問題には必ず,「$\theta<<1$と仮定する」と記載されている. これにより, $\sin\theta\approx\theta$となり, 運動方程式は,
$$ \frac{d^2\theta}{dt^2}+\frac{g}{l}\theta = 0 $$
と定数係数の線形微分方程式が得られる. これをMATLAB Symbolic Math toolboxを用いると,

matlab

% Make symbols
syms theta(t) l g

% Differential Equation
dtheta = diff(theta,t);
eqn_thry = diff(dtheta,t) + g/l*theta == 0

となる.

今回解く問題と初期条件

今回は, 以下の条件を設定した.
* 糸の長さ: $l = 2$[m]
* 重力加速度: $g = 9.81$[m/s$^2$]
* 初期角度: $\theta = \frac{\pi}{6}$
* 初期角速度: $\frac{d\theta}{dt} = 0$

matlab

% Setting
l = 2; % length [m]
g = 9.81; % Gravity Acceleration [m/s^2]

% Show Equation
eqn_thry = subs(eqn_thry)

% Initial Condition
theta0 = pi/6; % [rad]
dtheta0 = 0; % [rad/s]

% Make Initial Boundary Condition
cond = [theta(0) == theta0,dtheta(0) == dtheta0]

この結果, 解くべき方程式は,

$$　\frac{\partial^2}{\partial t^2}\theta(t)+\frac{981\theta(t)}{200} = 0 $$
$$ \left(\theta(0)=\frac{\pi}{6}~~\left(\left(\left.\frac{\partial}{\partial t}\theta(t)\right)\right|_{t=0}\right)=0\right) $$

となる.

微分方程式の求解

上記の微分方程式をdsolve関数で解くと

matlab

eqn_thry_result(t) = dsolve(eqn_thry,cond)

より,
$$ \theta = \frac{\pi\cos\left(\frac{3\sqrt{218}}{20}t\right)}{6} $$

となる.

近似を除くと..

上の理論解は, 以下の近似により導出されています.
$$ \theta << 1 $$
では, この近似を用いなければどのような解になるのでしょうか?
まず, 運動方程式は,
$$ \frac{d^2\theta}{dt^2}+\frac{g}{l}\sin\theta = 0 $$
である. これを数値的に解くには, MATLABでのode45関数を用いる.

2 階微分方程式の数値的な求解

matlab

% Differential Equation
eqn_nmr = diff(theta,2) + g/l*sin(theta)
[V] = odeToVectorField(eqn_nmr);
M = matlabFunction(V,'vars',{'t','Y'});
sol = ode45(M,[0 stop],[theta0,dtheta0]);
eqn_nmr_result = deval(sol,show_t);

$$ \frac{\partial^2}{\partial t^2}\theta(t)+\frac{981\theta(t)}{200} = 0 $$

を実行する. そして, この結果を可視化して, 近似したときと比べると以下のようになる.

matlab

stop = 90; % Time range from 0 to this [s]
show_t = 0:stop/799:stop;
shown_x = l*sin(eqn_nmr_result(1,:));
shown_y = -l*cos(eqn_nmr_result(1,:));

% Create Animation
figure
hold on
plot(0,0,'ko','MarkerEdgeColor','none','MarkerfaceColor','k','MarkerSize',10);
p1 = plot([0,show_x(1)],[0,show_y(1)],'r-','LineWidth',1.2);
p2 = plot(show_x(1),show_y(1),'ro','MarkerEdgeColor','none','MarkerFaceColor','r','MarkerSize',15);
p3 = plot([0,shown_x(1)],[0,shown_y(1)],'b-','LineWidth',1.2);
p4 = plot(shown_x(1),shown_y(1),'bo','MarkerEdgeColor','none','MarkerFaceColor','b','MarkerSize',15);
hold off
axis equal
grid on
title('Theoretical Solution vs Numerical Solution');
xlabel('x');
ylabel('y');
xlim([-1.1*l 1.1*l]);
ylim([-1.1*l 1.1*l]);
legend([p2,p4],{'Theoretical Solution','Numerical Solution'});

for i=1:length(show_t)
    p1.XData = [0,show_x(i)];
    p1.YData = [0,show_y(i)];
    p2.XData = show_x(i);
    p2.YData = show_y(i);
    p3.XData = [0,shown_x(i)];
    p3.YData = [0,shown_y(i)];
    p4.XData = shown_x(i);
    p4.YData = shown_y(i);
    drawnow
end

MATLABは数値解析ソフトウェアです!!

ということで, 大学入試では必ずある, sinの近似を取り除いた微分方程式の求解に本領を発揮していただきました.

これも,Symbolic Math Toolboxで簡単!!#MATLAB #uec #数値解析 pic.twitter.com/Y409H2LBlX
— ヒビヤ@MATLAB Ambassador (@hibs_MATLAB_Amb) September 2, 2021

おわりに

私は, 高校生の時までプログラミングなど全く経験したことがありませんでした. その後, 大学に入学して2~3つのプログラミング言語に触れましたが, MATLABはその中でも, 数式が自然に書けるという点においてとても秀でており, 練習がとてもしやすかった思い出があります.

皆さんも, 今までやったことのある物理や数学の問題をMATLABで解くことでプログラミングの練習してみるのはいかがですか!?

もっとうまくSymbolic Math Toolboxをうまく使って二重振子運動について考えている例もありました.

MATLAB Documentation - 二重振子運動のアニメーションと解

大学受験で出てくる単振り子問題の近似とは...?

大学受験に出る近似を含んだ単振り子問題

運動方程式と理論解

今回解く問題と初期条件

微分方程式の求解

近似を除くと..

おわりに