空中に投げられたバトンの動きを表す方程式の求解

ライブスクリプトを開く

この例では、空中に投げられたバトンのダイナミクスをモデル化する連立常微分方程式を解きます [1]。バトンは、長さ $L$ のロッドで結合された質量 $m_{1}$ および $m_{2}$ の 2 つの質点としてモデル化されます。バトンは空中に投げられた後、重力の影響を受けながら、垂直な xy 平面内を移動します。ロッドは水平方向と $θ$ の角度をもち、1 つ目の質点の座標は $(x, y)$ です。この定式化では、2 つ目の質点の座標は $(x + L \cos θ, y + L \sin θ)$ です。

このシステムの運動方程式は、3 つの座標 $x$ 、 $y$ 、および $θ$ それぞれに対しラグランジュ方程式を適用することで得られます。

$\begin{array}{l} (m_{1} + m_{2}) \ddot{x} - m_{2} L \ddot{θ} \sin θ - m_{2} L {\dot{θ}}^{2} \cos θ = 0, \\ (m_{1} + m_{2}) \ddot{y} - m_{2} L \ddot{θ} \cos θ - m_{2} L {\dot{θ}}^{2} \sin θ + (m_{1} + m_{2}) g = 0, \\ L^{2} \ddot{θ} - L \ddot{x} \sin θ + L \ddot{y} \cos θ + g L \cos θ = 0 . \end{array}$

MATLAB® でこの ODE 系を解くには、ソルバー ode45 を呼び出す前に方程式を関数へとコード化します。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

方程式のコード作成

ode45 ソルバーでは、 $\dot{q} = f (t, q)$ の形式で方程式を記述する必要があります。ここで、 $\dot{q}$ は各座標の 1 階微分です。この問題では、解ベクトルに $x$ 、 $y$ 、角度 $θ$ 、およびそれらの 1 階微分という 6 つの要素があります。

$q = [\begin{array}{c} q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ q_{5} \\ q_{6} \end{array}] = [\begin{array}{c} x \\ \dot{x} \\ y \\ \dot{y} \\ θ \\ \dot{θ} \end{array}] .$

この表記を使って、運動方程式全体を $q$ の要素によって書き換えることができます。

$\begin{array}{l} (m_{1} + m_{2}) \dot{q_{2}} - m_{2} L \dot{q_{6}} \sin q_{5} = m_{2} L q_{6}^{2} \cos q_{5}, \\ (m_{1} + m_{2}) \dot{q_{4}} - m_{2} L \dot{q_{6}} \cos q_{5} = m_{2} L q_{6}^{2} \sin q_{5} - (m_{1} + m_{2}) g, \\ L^{2} \dot{q_{6}} - L \dot{q_{2}} \sin q_{5} + L \dot{q_{4}} \cos q_{5} = - g L \cos q_{5} . \end{array}$

左辺に 1 階微分の項がいくつか存在するため、残念ながらこの運動方程式はソルバーが必要とする形式 $\dot{q} = f (t, q)$ に適合しません。こうした場合は、質量行列を使用して方程式の左辺を表現しなければなりません。

行列表記では、 $M (t, q) \dot{q} = f (t, q)$ の形式で質量行列を使用して、6 つの方程式からなる方程式系として運動方程式を書き換えることができます。質量行列では、方程式の左辺にある 1 階微分の線形結合が、行列とベクトルの積により表現されます。この形式では、方程式系は次のようになります。

$[\begin{array}{c} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & m_{1} + m_{2} & 0 & 0 & 0 & - m_{2} L \sin q_{5} \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & m_{1} + m_{2} & 0 & m_{2} L \cos q_{5} \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & - L \sin q_{5} & 0 & L \cos q_{5} & 0 & L^{2} \end{array}] [\begin{array}{c} \dot{q_{1}} \\ \dot{q_{2}} \\ \dot{q_{3}} \\ \dot{q_{4}} \\ \dot{q_{5}} \\ \dot{q_{6}} \end{array}] = [\begin{array}{c} q_{2} \\ m_{2} L q_{6}^{2} \cos q_{5} \\ q_{4} \\ m_{2} L q_{6}^{2} \sin q_{5} - (m_{1} + m_{2}) g \\ q_{6} \\ - g L \cos q_{5} \end{array}]$

この式から、質量行列の非ゼロ要素を計算する関数を記述できます。関数は 3 つの入力を取ります。 $t$ と解ベクトル $q$ は必須 (関数本体で使用されない場合でも、これらの入力を指定しなければなりません) であり、 $P$ はオプションの追加入力で、パラメーター値を渡すために使用されます。この問題のパラメーターを関数に渡すには、P をパラメーター値を保持する構造体として作成し、関数のパラメーター化で説明されている手法を使用して、その構造体を追加入力として関数に渡します。

function M = mass(t,q,P)
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Mass matrix function
M = zeros(6,6);
M(1,1) = 1;
M(2,2) = m1 + m2;
M(2,6) = -m2*L*sin(q(5));
M(3,3) = 1;
M(4,4) = m1 + m2;
M(4,6) = m2*L*cos(q(5));
M(5,5) = 1;
M(6,2) = -L*sin(q(5));
M(6,4) = L*cos(q(5));
M(6,6) = L^2;
end

次に、方程式系 $M (t, q) \dot{q} = f (t, q)$ における各方程式の右辺の関数を記述できます。質量行列関数と同様、この関数は $t$ と $q$ という 2 つの必須入力と、パラメーター値 $P$ という 1 つのオプション入力を取ります。

function dydt = f(t,q,P)
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Equation to solve
dydt = [q(2)
        m2*L*q(6)^2*cos(q(5))
        q(4)
        m2*L*q(6)^2*sin(q(5))-(m1+m2)*g
        q(6)
        -g*L*cos(q(5))];
end

方程式を解く

まず、適切に命名されたフィールドを構造体に設定することで、 $m_{1}$ 、 $m_{2}$ 、 $g$ 、および $L$ のパラメーター値からなる構造体 P を作成します。構造体 P は追加入力として、質量行列と ODE 関数に渡されます。

P.m1 = 0.1;
P.m2 = 0.1;
P.L = 1;
P.g = 9.81

P = struct with fields:
    m1: 0.1000
    m2: 0.1000
     L: 1
     g: 9.8100

積分の時間範囲として、0 ～ 4 の 25 点からなるベクトルを作成します。時間のベクトルを指定すると、ode45 は要求された時点での内挿された解を返します。

tspan = linspace(0,4,25);

問題の初期条件を設定します。バトンはある角度で上向きに投げられるため、初期速度には非ゼロ値を $\dot{x_{0}} = 4$ 、 $\dot{y_{0}} = 20$ 、 $\dot{θ_{0}} = 2$ のように使用します。バトンは直立位置からの開始となるため、初期位置は $θ_{0} = - π / 2$ 、 $x_{0} = 0$ 、 $y_{0} = L$ となります。

y0 = [0; 4; P.L; 20; -pi/2; 2];

odeset を使用して質量行列関数を参照するオプション構造体を作成します。ソルバーでは、質量行列関数が入力を 2 つのみもつものと想定されているため、無名関数を使用して、ワークスペースからパラメーター P を渡します。

opts = odeset('Mass',@(t,q) mass(t,q,P));

最後に、次の入力を指定し、ode45 を使用して方程式系を解きます。

方程式 @(t,q) f(t,q,P) の無名関数
解が要求される時間のベクトル tspan
初期条件のベクトル y0
質量行列を参照するオプション構造体 opts

[t,q] = ode45(@(t,q) f(t,q,P),tspan,y0,opts);

結果のプロット

ode45 からの出力には、要求された各タイムステップでの運動方程式の解が含まれます。結果を確認するために、時間経過に沿ったバトンの運動をプロットします。

解の各行をループし、タイムステップごとにバトンの位置をプロットします。時間経過に沿った回転が確認できるよう、バトンのそれぞれの端に色を付けます。

figure
title('Motion of a Thrown Baton, Solved by ODE45');
axis([0 22 0 25])
hold on
for j = 1:length(t)
   theta = q(j,5);
   X = q(j,1);
   Y = q(j,3);
   xvals = [X X+P.L*cos(theta)];
   yvals = [Y Y+P.L*sin(theta)];
   plot(xvals,yvals,xvals(1),yvals(1),'ro',xvals(2),yvals(2),'go')
end
hold off

Figure contains an axes object. The axes object with title Motion of a Thrown Baton, Solved by ODE45 contains 75 objects of type line.

参照

[1] Wells, Dare A. Schaum’s Outline of Theory and Problems of Lagrangian Dynamics: With a Treatment of Euler’s Equations of Motion, Hamilton’s Equations and Hamilton’s Principle. Schaum's Outline Series. New York: Schaum Pub. Co, 1967.

ローカル関数

ここでは、ODE ソルバーが解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。

function dydt = f(t,q,P) % Equations to solve
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% RHS of system of equations
dydt = [q(2)
        m2*L*q(6)^2*cos(q(5))
        q(4)
        m2*L*q(6)^2*sin(q(5))-(m1+m2)*g
        q(6)
        -g*L*cos(q(5))];
end
%------------------------------------------------
function M = mass(t,q,P) % Mass matrix function
% Extract parameters
m1 = P.m1;
m2 = P.m2;
L = P.L;
g = P.g;

% Set nonzero elements in mass matrix
M = zeros(6,6);
M(1,1) = 1;
M(2,2) = m1 + m2;
M(2,6) = -m2*L*sin(q(5));
M(3,3) = 1;
M(4,4) = m1 + m2;
M(4,6) = m2*L*cos(q(5));
M(5,5) = 1;
M(6,2) = -L*sin(q(5));
M(6,4) = L*cos(q(5));
M(6,6) = L^2;
end

参考

ode45

空中に投げられたバトンの動きを表す方程式の求解

方程式のコード作成

方程式を解く

結果のプロット

参照

ローカル関数

参考

関連するトピック