Return to 単振り子のエネルギー保存則から周期を求める準備

Maxima で単振り子のエネルギー保存則から周期を求める

単振り子の運動方程式から得られるエネルギー保存則から,Maxima を使って数値積分によって周期を求める。導出については以下を参照:

数値積分 quad_qags()

振幅 $\theta_0$ の単振り子の周期 $T_p$ は,エネルギー保存則の積分から以下のようになる。

\begin{eqnarray}
T_p(\theta_0) &=&
\frac{2}{\pi} \int_{0}^{\theta_0} \frac{1}{\sqrt{2(\cos\theta-\cos\theta_0)}}d\theta \\
&=&
\frac{\sqrt{2}}{\pi} \int_{0}^{\theta_0} \frac{1}{\sqrt{\cos\theta-\cos\theta_0}}d\theta
\end{eqnarray}

積分部分は数値積分 quad_qags() を使う。(romberg() ではなく quad_qags() を使う理由はすぐわかります。)

この積分が一見やっかいなのは,積分範囲の上限 $\theta = \theta_0$ で被積分関数の分母がゼロになることです。$\theta = \theta_0$ で被積分関数が発散してしまうので,何かエラーでも起こるのかと思いきや,quad_qags() は大変優秀で何事もなかったかのように数値積分してしまいます!

ためしに $\theta_0 = 60^{\circ}$ の場合の積分値を数値積分 quad_qags() で求めてみます。

In [1]:
/* 度をラジアンに変換する関数を準備しておく */
radians(deg):= deg * float(%pi)/180$

th: 60$
theta0: radians(th)$

quad_qags(1/sqrt(cos(theta)-cos(theta0)), theta, 0, theta0);
Out[1]:
\[\tag{${\it \%o}_{4}$}\left[ 2.384011014402899 , 3.274180926382542 \times 10^{-10} , 315 , 0 \right] \]

quad_qags() は4つの要素のリストを返します:

積分の近似値
見積もられた近似の絶対誤差
非積分関数の評価数
エラーコード

積分を使って定義された周期を Tp1(th) と定義する。

In [2]:
Tp1(th):= block(
  /* 度で入力した th をラジアンになおす */
  theta0: radians(th), 
  
  /* quad_gags() で数値積分し,tmp に格納する */
  tmp: quad_qags(1/sqrt(cos(theta)-cos(theta0)), theta, 0, theta0), 
  
  /* 積分の近似値を tmp[1] として取り出し,係数をかけて周期をかえす */
  return(sqrt(2.0)/float(%pi) * tmp[1])
)$
In [3]:
/* 振幅と周期のグラフを描くために x, y に追加していく */
x: []$ /* 振幅 */
y: []$ /* 周期 */

/* 参考として,ほぼ単振動となるであろう θ_0 = 1°の場合 */
th: 1$
x: append(x, [th])$
y: append(y, [Tp1(th)])$

/* リストに append したばかりの最後の値を表示 */
printf(true, "θ_0 = ~2d° のとき  ", x[length(x)])$
printf(true, "周期 T = ~8,6f~%", y[length(x)])$

/* θ_0 を 10° から 90° まで 10° ごとにかえて計算 */
for th:10 thru 90 step 10 do(
  x: append(x, [th]),
  y: append(y, [Tp1(th)]),
  printf(true, "θ_0 = ~2d° のとき  ", x[length(x)]), /* 改行しない */
  printf(true, "周期 T = ~8,6f~%", y[length(x)])     /* ~% で改行  */
)$
θ_0 =  1° のとき  周期 T = 1.000019
θ_0 = 10° のとき  周期 T = 1.001907
θ_0 = 20° のとき  周期 T = 1.007669
θ_0 = 30° のとき  周期 T = 1.017409
θ_0 = 40° のとき  周期 T = 1.031341
θ_0 = 50° のとき  周期 T = 1.049783
θ_0 = 60° のとき  周期 T = 1.073182
θ_0 = 70° のとき  周期 T = 1.102145
θ_0 = 80° のとき  周期 T = 1.137493
θ_0 = 90° のとき  周期 T = 1.180341
In [4]:
/* 振幅と周期のグラフを描く */
draw2d(
  font = "Arial", font_size = 16, 
  title = "単振り子の振幅と周期の関係",
  /* 表示範囲 */
  xrange = [0, 92], yrange = [0.95, 1.2], 
  /* グリッド  */
  grid = true, 
  /* 座標軸のラベル */
  xlabel = "振幅 θ (°)", ylabel = "規格化された周期 T", 
  
  /* 点を塗りつぶした丸に */
  point_type = 7,
  /* 点の大きさ */
  point_size = 0.5,
  points(x, y)
)$

第1種完全楕円積分 elliptic_kc(m)

ありがたいことに,Maxima では以下のように定義された第一種完全楕円積分が elliptic_kc(m) として使えます。

$$\mbox{elliptic_kc(m)} = K(m) \equiv \int_0^{\frac{\pi}{2}}\frac{1}{\sqrt{1 – m \sin^2 x}} dx$$

第1種完全楕円積分 elliptic_kc(m) を使って,振幅 $\theta_0$ の単振り子の周期 $T_p$ を以下のように定義する。($m = k^2$ に注意。)

In [5]:
Tp(th):= block(
  theta0: radians(th), 
  m: sin(theta0/2)**2, 
  return(2/float(%pi) * elliptic_kc(m))
)$
In [6]:
/* 参考として,ほぼ単振動となるであろう θ_0 = 1°の場合を表示 */
th: 1$
printf(true, "θ_0 = ~2d° のとき  ", th)$
printf(true, "周期 T = ~14,12f~%", Tp(th))$

/* θ_0 を 10° から 90° まで 10° ごとにかえて表示 */
for th:10 thru 90 step 10 do(
  printf(true, "θ_0 = ~2d° のとき  ", th),
  printf(true, "周期 T = ~14,12f~%", Tp(th))
)$
θ_0 =  1° のとき  周期 T = 1.000019038921
θ_0 = 10° のとき  周期 T = 1.001907188143
θ_0 = 20° のとき  周期 T = 1.007669025792
θ_0 = 30° のとき  周期 T = 1.017408797596
θ_0 = 40° のとき  周期 T = 1.031340519130
θ_0 = 50° のとき  周期 T = 1.049782960623
θ_0 = 60° のとき  周期 T = 1.073182007149
θ_0 = 70° のとき  周期 T = 1.102144909639
θ_0 = 80° のとき  周期 T = 1.137492559924
θ_0 = 90° のとき  周期 T = 1.180340599016
In [7]:
/* 振幅と周期のグラフを描く */
draw2d(
  font = "Arial", font_size = 16, 
  title = "単振り子の振幅と周期の関係",
  /* 表示範囲 */
  xrange = [0, 90], yrange = [0.95, 1.2], 
  /* グリッド  */
  grid = true, 
  /* 座標軸のラベル */
  xlabel = "振幅 θ (°)", ylabel = "規格化された周期 T", 
  
  /* 陽関数 y = f(x) としてグラフにする */
  explicit(Tp(t), t, 1, 90), file_name = mfuriE02A, terminal = 'svg
)$