最速降下曲線と数値計算

はじめに

斜面を最速で転がり切ることができる曲線は、サイクロイド曲線である。これは、変分法により求めることができる。変分法のイメージとしては、斜面の曲線を少し変化させた場合、転がり切るまでにかかる時間がどの程度変化するのかを考えるといった感じだ。そこで、今回はそのようなアルゴリズムをプログラムで再現することを目的とする。そして、数値計算で得られた斜面の形状とサイクロイド曲線を比較することでそのプログラムの妥当性を検証する。

問題設定

$xy$座標平面上で考える。$(a,0),(0,a)$を結ぶなだらかな曲線を考える。これを斜面として、$(0,a)$からボールを$(a,0)$まで転がす。ただし、$-y$方向に重力がかかっているものとする。
このとき、曲線上任意の点$(x,y)$にあるボールの速さは以下のように表すことができる。

v=\sqrt{2g(a-y)}

ここで、$0\le x\le a$までの範囲を$n$等分することを考える。

区間$x_k\le x\le x_{k+1}$について、ボールの平均の速さは以下の式で与えられる。

v_k=\sqrt{2g\{a-(\frac{y_k+y_{k+1}}{2})\}}

一方でその範囲の斜面の長さ$\Delta s$は、$n$が十分に大きいとき直線にみなせるため、

\Delta s_k=\sqrt{(\Delta x_k)^2+(\Delta y_k)^2}

したがって、斜面を転がり切るまでにかかる時間$t$は以下のように与えられる。

t=\lim_{n\to \infty} \frac{\Delta s_k}{v_k}

したがって、$t$を最小にするための$y$の形状について最急降下法を用いて推定する。

プログラム

そこで、以下のようなプログラムを作成した。

python t_min.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#微小変化量
h=1e-5
#学習率
alpha=0.01
#斜面の始点及び終点の座標
a=10
b=a

#斜面の点群の数
n=20
#重力加速度
g=0.01
x_ary=np.linspace(0,a,n)
y_ary=-(b/a)*(x_ary-a)

def t_gra(y_ary):
    t = 0
    for i in range(1,n):
        dx=x_ary[i]-x_ary[i-1]
        dy=y_ary[i]-y_ary[i-1]
        # 始点からの落下高さを使用（y座標が減少する方向が下）
        h_drop = y_ary[0] - (y_ary[i]+y_ary[i-1])/2
    
        v=(2*g*h_drop)**0.5
        ds=(dx**2+dy**2)**0.5
        dt=ds/v
        t=t+dt
    return t

m_ary=[]
y_ary_t=[]
t_ary=[]
ims = []
m=10000
for p in range(1,m):
    # 始点(k=0)と終点(k=n-1)を固定、内部点のみ更新
    for k in range(1, n-1):
        t_down=t_gra(y_ary)
        y_ary[k]=y_ary[k]+h
        t_up=t_gra(y_ary)
        d_y=(t_up-t_down)/h
        y_ary[k]=y_ary[k]-alpha*d_y
        y_ary[k]=y_ary[k]-h
    y_ary_t.append(np.copy(y_ary))
    t_ary.append(t_gra(y_ary))
    m_ary.append(p)
#試行回数と時間
plt.plot(m_ary,t_ary)
plt.xlabel("試行回数")
plt.ylabel("斜面を転がりきるまでにかかる時間")
plt.savefig("time_min.png")
plt.show()
#斜面の比較
plt.scatter(x_ary,y_ary,color="blue",label="最適計算後の斜面の点群")
y_ary=-(b/a)*(x_ary-a)
plt.plot(x_ary,y_ary,color="black",label="初期直線")
theta_lim=math.pi
theta_ary=np.linspace(0,theta_lim,100)
xx=a*(theta_ary-np.sin(theta_ary))/theta_lim
yy=-b*(1-np.cos(theta_ary))/2+b
plt.plot(xx,yy,color="red",label="サイクロイド曲線")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.savefig("最適化後の斜面.png")
plt.show()

結果

n=20

n=40

n=100

結果として、転がり始めの誤差が大きかった。理由としては最初は斜面の傾きが極めて大きくなるため加速しやすいからだと考えられる。

まとめ

今回は、変分法でおなじみの最速降下曲線を再急降下法を用いた数値計算で検証してみた。結果、アルゴリズムの精度もしくは計算量のせいで上手くサイクロイド曲線を描写することができなかった。ただし、加速度や傾きが小さい区間ではかなり高い精度で推定することができた。