フーコーの振子のモデル化

モデルを開く

この例では、フーコーの振子のモデル化方法を示します。フーコーの振子は、フランスの物理学者レオン・フーコーが考えたものです。その目的は、地球の自転を証明することでした。地球は自転しているため、フーコーの振子の振動面は終日回転します。振動面が 1 周するのにかかる時間間隔 T は、地理緯度によって異なります。

フーコーの振子のうち最も有名なのは、パリのパンテオンに取り付けられたものです。28 kg の金属球が 67 m ものワイヤーでつるされました。この例では、パリの地理緯度で 67 m の振子のシミュレーションを実行します。

Simulink® モデル

フーコーの振子の問題を Simulink® で解く最も簡単な方法は、このシステムの連立微分方程式を解くモデルを作成することです。このモデルを図 1 に示します。フーコーの振子を説明する方程式は以下のとおりです。このモデルの物理現象とこれらの方程式の導出の詳細については、解析と物理学を参照してください。

$
\ddot{x} = 2\Omega \sin{\lambda} \dot{y} - \frac{g}{L} x
$

$
\ddot{y} = - 2\Omega\sin{\lambda} \dot{x} - \frac{g}{L} y
$

$x, y = \mbox{ pendulum bob coordinates as seen by an observer on Earth}$

$\Omega = \mbox{ Earth's angular speed of revolution about its axis } (rad/sec)$

$g = \mbox{ acceleration of gravity } (m/sec^2)$

$L = \mbox{ length of the pendulum string } (m)$

$\lambda = \mbox{ geographic latitude } (rad)$

モデルを開く

MATLAB® コマンドウィンドウで「sldemo_foucault」と入力してこのモデルを開きます。このモデルはシミュレーションデータのログを変数 sldemo_foucault_output に作成します。ログが作成された信号には青いインジケーターが付きます。詳細については、ログ記録用の信号のマークを参照してください。

図 1: フーコーの振子モデル

初期条件

このモデルは、定数と初期条件を sldemo_foucault_data.m ファイルから読み込みます。このファイルの内容は、以下の表 1 のとおりです。シミュレーションパラメーターを MATLAB ワークスペースで直接変更できます。振子の初期振幅は、振子の長さよりも小さくなければなりません。微分方程式は小さな振動に対してのみ有効であるためです。

表 1: 初期条件

g = 9.83;          % acceleration of gravity (m/sec^2)
L = 67;            % pendulum length (m)
initial_x = L/100; % initial x coordinate (m)
initial_y = 0;     % initial y coordinate (m)
initial_xdot = 0;  % initial x velocity (m/sec)
initial_ydot = 0;  % initial y velocity (m/sec)
Omega=2*pi/86400;  % Earth's angular velocity of rotation about its axis (rad/sec)
lambda=49/180*pi;  % latitude in (rad)

シミュレーションの実行

モデルウィンドウのツールバーの [再生] ボタンをクリックしてシミュレーションを実行します。このシミュレーションでは、可変ステップスティッフソルバー ode23t が使用されます。フーコーの振子は 3600 秒間シミュレートされます (シミュレーション時間は変更可能)。このモデルでは、既定の相対許容誤差 RelTol = 1e-6 が使用されます。

図 2: フーコーの振子のシミュレーション結果 (シミュレーション時間は 3600 秒)

結果

シミュレーション結果は、上記の図 2 のようになります。振子の x 座標および y 座標と、x 軸方向および y 軸方向の速度成分が計算されます。

振子の振動面が 360°回転するには 24 時間超かかります。振動面が 1 周するのにかかる時間は、地理緯度 lambda の関数です (解析と物理学での導出を参照)。

図 3: 振子の振動面が 1 時間でどれだけ回転するかを示す Animation ブロック

シミュレーションを実行した後、Animation ブロックをダブルクリックして結果をアニメーションで表示します。

メモ: 例の "結果をアニメーションで表示" 部分では、Signal Processing Toolbox™ が必要です。これがインストールされていない場合、Animation ブロックをダブルクリックすると、エラーが表示されます。例のその他の部分は、Signal Processing Toolbox がなくても正しく機能します。

sldemo_foucault_animate.m ファイルでは、振子の振り玉の位置がさまざまな時点でプロットされています。振子の振動面が回転する様子をはっきり確認できます。

メモ: 相対許容誤差を既定値よりも大きくしてシミュレーションを実行している場合、結果は長時間にわたって数値的に不安定になります。必ず、スティッフな可変ステップソルバーを使用してください。スティッフな問題の数値的不安定性とソルバーの性能の詳細については、「スティッフなモデルを使った可変ステップソルバーの調査」の例を参照してください。

モデルを閉じる

モデルを閉じます。生成されたデータをクリアします。

解析と物理学

この節では、フーコーの振子を解析し、その背景にある物理現象を説明します。振子は、長さ L のワイヤーにつるされた質点としてモデル化できます。振子の地理緯度は lambda です。図 4 の基準座標系を使用すると便利です。慣性座標系 I (地球の中心に対する) と非慣性座標系 N (地表上の観測者に対する)回転により、非慣性座標系は加速します。

図 4: この問題の慣性座標系と非慣性座標系

点 O は非慣性座標系 N の原点です。これは、振子のつり下げ点の真下にある地表上の点です。非慣性座標系は、z 軸が地球の中心とは逆の、地表に垂直な方向に向くように選択されます。x 軸は南に向かい、y 軸は西に向かいます。

最初に説明したとおり、フーコーの振子の振動面は回転します。振動面は、以下の式で与えられる時間 Trot で 1 周します。Tday は 1 日 (つまり、地球が 1 回自転するのにかかる時間) です。

$T_{rot}=T_{day} \cdot \sin \lambda$

正弦係数についてはさらに説明が必要です。振子の振動面は地球の中心に対する慣性座標系で固定されていると仮定されることが多いのですが、これは誤りです。それが正しいのは北極と南極のみです。このような混同を除外するために、振子がつり下げられている点 S (図 4 を参照) について考えてみましょう。慣性座標系 I の場合、点 S は円を描きます。振子の振り玉は、一定の長さのワイヤーでつるされています。簡単にするために、空気摩擦は無視します。慣性座標系 I の場合、振り玉にかかる力は、ワイヤーの張力 T と重力 Fg の 2 つのみです。

ベクトル r は、振子の振り玉の位置 B (図 4 を参照) を示します。ニュートンの第 2 法則によれば、物体にかかるすべての力の合計は、その物体の質量に加速度を乗算したものと等しくなります。

$m \ddot{\overrightarrow{r}} = \overrightarrow{T} + \overrightarrow{F_g}$

この証明では、ドットは時間微係数を表し、矢印はベクトルを表し、キャップは単一ベクトル (x、y、z 軸方向に i、j、k) を表します。ベクトル矢印の上のドットは、そのベクトルの時間微係数を示します。ドットの上の矢印は、時間微係数のベクトルを示します。合計加速度と半径方向加速度の違いを以下で確認してください。

合計加速度:

$
\ddot{\overrightarrow{r}} = \frac{d^2 \overrightarrow{r}}{d t^2}=
\frac{d^2}{d t^2} \left( x\mathbf{\hat{i}} + y\mathbf{\hat{j}} + z\mathbf{\hat{k}} \right)
$

半径方向加速度:

$
\overrightarrow{ \ddot{r} } = \overrightarrow{ \left( \frac{d^2 r}{dt^2}\right)}=
\ddot{x} \mathbf{\hat{i}}+
\ddot{y} \mathbf{\hat{j}}+
\ddot{z} \mathbf{\hat{k}}
$

重力加速度は地球の中心に向かいます (z 軸の負の方向)。

$$g=9.83 m/sec^2$$

$\overrightarrow{g} = -g\mathbf{\hat{k}}$

$
m \ddot{\overrightarrow{r}} =
\overrightarrow{T} - mg\mathbf{\hat{k}}
$

加速度の項を以下のように分解します。

$
\overrightarrow{r}=
r_x \mathbf{\hat{i}}+
r_y \mathbf{\hat{j}}+
r_z \mathbf{\hat{k}}
$

$
\dot{\overrightarrow{r}}=
\left(
\dot{r}_x \mathbf{\hat{i}}+
\dot{r}_y \mathbf{\hat{j}}+
\dot{r}_z \mathbf{\hat{k}}
\right) +
r_x \dot{ \mathbf{\hat{i}} } +
r_y \dot{ \mathbf{\hat{j}} } +
r_z \dot{ \mathbf{\hat{k}} }
$

非慣性基準座標系 N が空中を回っているため、単位ベクトルの時間微係数が表示されます。つまり、単一ベクトル i、j、k は空中を回っていることを意味します。時間微係数は以下のように与えられます。Omega は地球の自転の角速度です。スカラー Omega は角速度の値です。ベクトル Omega はベクトル角速度です。その方向は右手の法則で決まります。

$
\dot{\mathbf{\hat{i}}}=\overrightarrow{\Omega} \times \mathbf{\hat{i}}
$

$
\dot{\mathbf{\hat{j}}}=\overrightarrow{\Omega} \times \mathbf{\hat{j}}
$

$
\dot{\mathbf{\hat{k}}}=\overrightarrow{\Omega} \times \mathbf{\hat{k}}
$

Omega に対するベクトル r の時間微係数を書き換えます。

$
\dot{\overrightarrow{r}}=
\left(
\dot{r}_x \mathbf{\hat{i}}+
\dot{r}_y \mathbf{\hat{j}}+
\dot{r}_z \mathbf{\hat{k}}
\right) +
\overrightarrow{\Omega} \times \overrightarrow{r}=
\overrightarrow{\dot{r}}+
\overrightarrow{\Omega} \times \overrightarrow{r}
$

同様に、ベクトル r の 2 つ目の時間微係数を表現します。

$
\ddot{\overrightarrow{r}}=
\overrightarrow{\ddot{r}} +
2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}} +
\overrightarrow{\Omega} \times
\left( \overrightarrow{\Omega} \times \overrightarrow{r} \right)
$

$\overrightarrow{\dot{r}} = \mbox{ acceleration in the non-inertial frame N (x, y, z components)}$

$2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}} = \mbox{ Coriolis acceleration}$

$\overrightarrow{\Omega} \times
\left( \overrightarrow{\Omega} \times \overrightarrow{r} \right)
= \mbox{ additional term due to rotation of non-inertial frame N}
$

この方程式を簡単にするために、地球の Omega が非常に小さいと仮定します。そうすることで、上記の方程式の 3 番目の項を無視できます。実際、2 つ目の項 (1 つ目の項よりも既にかなり小さい) は 3 つ目の項よりも 4 けた大きくなっています。これにより、方程式は以下のような形になります。

$
\ddot{\overrightarrow{r}} \simeq
\overrightarrow{\ddot{r}} +
2 \overrightarrow{\Omega} \times \overrightarrow{\dot{r}}
$

ニュートンの第 2 法則を書き換え、以下のように x、y、z 成分に分解できます。

$
m \overrightarrow{\ddot{r}} =
\overrightarrow{T} - mg\mathbf{\hat{k}} -
2 m \overrightarrow{\Omega} \times \overrightarrow{\dot{r}}
$

$
m \ddot{x} = T_x + 2m\Omega \dot{y} \sin{\lambda}
$

$
m \ddot{y} = T_y - 2m\Omega \left(\dot{x} \sin{\lambda}+\dot{z}\cos{\lambda}\right)
$

$
m \ddot{z} = T_z - mg + 2m \Omega \dot{y} \cos{\lambda}
$

振動の角度振幅はわずかです。したがって、垂直速度と垂直加速度 (z ドットと z 二重ドット) は無視できます。ストリング張力の成分は小角度近似を使用して表現でき、問題が 2 次元へと大幅に単純化されます (以下を参照)。

$T_z = mg - 2m\Omega \dot{y} \cos{\lambda} \simeq mg$

$T_x= -T\frac{x}{L}$

$T_y= -T\frac{y}{L}$

$T_z= T\frac{L-z}{L}\simeq T$

特性微分方程式

最後に、この問題の物理現象は以下に示す連立方程式系で説明できます。x 座標と y 座標は、地球上の観測者から見た振子の振り玉の位置を指定しています。

$
 \ddot{x} - 2\Omega \sin{\lambda} \dot{y} + \frac{g}{L} x =0
$

$
 \ddot{y} + 2\Omega\sin{\lambda} \dot{x} + \frac{g}{L} y =0
$

解析的な解 (近似)

以下は、フーコーの振子の問題に対する解析的な解です。残念ながら、厳密なものではありません。解析解を微分方程式に代入しようとすると、Omega 2 乗のオーダーの約されない項が残ります。しかし、Omega は非常に小さいため、実用目的では、約されていない項を無視できます。

$\eta = x+ i\cdot y \mbox{ (complex number)}$

$\ddot{\eta}+(2i\Omega \sin{\lambda})\dot{\eta} + \frac{g}{L} \eta = 0$

$\eta = \left( C_1 e^{i\sqrt{g/L}t} + C_2 e^{-i\sqrt{g/L}t}\right)
e^{-i\Omega t \sin{\lambda}}$

$C_1, C_2 \mbox{ are complex integration
constants}$

実際の微分方程式系は非対称

微分において、Omega の自乗に関する項は無視されました。これにより、微分方程式で xy 対称性が生じます。Omega の自乗に関する項を考慮した場合、微分方程式系は非対称になります (以下を参照)。

$
\ddot{x} - 2\Omega \sin{\lambda} \dot{y} + (\frac{g}{L}-\Omega^2 \sin^2{\lambda}) x =0
$

$
\ddot{y} + 2\Omega\sin{\lambda} \dot{x} + (\frac{g}{L} - \Omega^2) y =0
$

非対称の微分方程式が考慮されるように現在のフーコーの振子モデルを変更することは簡単です。g/L を含んでいる該当の Gain ブロックを編集し、必要な式を追加するだけです。この変更による数値結果の全体的な訂正は、ごくわずかです。