[2012/8/28] Brachistochrone Problem

작성자: 허필원(2012/8/28)

1년도 더 된것 같다. 뜻하기 않게 Brachistochrone 문제를 풀 기회가 (부탁을 받아서) 생겨서 여기에 정리해 본다.

문제 정의: 한 점에서 volume이 없는 particle 이 어떤 path 를 따라 중력의 작용에 의해서 움직인다. 이때, 특정 수평방향 위치까지 가장 빨리 가는 path 를 구하여라.

(아래 그림 참조)

문제설정: 시간을 최소화 하는 문제임. 먼저, 아래와 같이 간단히 좌표설정을 해야 한다.

Solution: Conservation of mechanical energy 에 의해서 다음과 같은 식을 얻을 수 있다.

\[ \begin{align*} mgy&=\frac{1}{2}mv^2 \\ 2gy&=v^2=\dot{x}^2+\dot{y}^2 \\ &=\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2 \\ &=\left[ \left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2 \right] \left( \frac{dt}{dx} \right)^2 \left( \frac{dx}{dt}\right)^2 \\ &=\left[ 1 +\left(y' \right)^2 \right] \left( \frac{dx}{dt}\right)^2 \\ \therefore \sqrt{2gy}&=\sqrt{1 +\left(y' \right)^2} \frac{dx}{dt} \\ \implies dt&=\frac{\sqrt{1 +\left(y' \right)^2}}{\sqrt{2gy}}dx \\ \therefore T&=\int_{0}^{b} \frac{\sqrt{1 +\left(y' \right)^2}}{\sqrt{2gy}}dx \end{align*} \]

이 문제는 T 를 최소화 하는 문제이며, Calculus of Variation 의 전통적인 문제중 하나가 된다. 즉, Euler-Lagrange (EL) 방정식을 아래와 같이 풀면 된다.

\[ \begin{align*} L=\frac{\sqrt{1 +\left(y' \right)^2}}{\sqrt{2gy}} \\ \textrm{EL eq:}\quad L_y-\frac{d}{dx}L_{y'}=0 \end{align*} \]

그런데, Lagrange \(L\) 을 잘 살펴보면, \(x\) 에 independent 함을 알 수 있고, 이 경우에는 Hamiltonian 이 \(x\) 에 대해 일정한 값을 갖게 된다. 이를 이용해서 계속 풀어보면 아래와 같다.

\[ \begin{align*} H&=y' L_{y'}-L=C \\ L_{y'}&=\frac{y'}{\sqrt{2gy}\sqrt{1 +\left(y' \right)^2}} \\ \therefore H&=\frac{\left(y'\right)^2}{\sqrt{2gy}\sqrt{1 +\left(y' \right)^2}}-\frac{\sqrt{1 +\left(y' \right)^2}}{\sqrt{2gy}} \\ &=-\frac{1}{\sqrt{2gy}\sqrt{1 +\left(y' \right)^2}} \\ &=C \\ \implies &y\left[ 1+\left(y' \right)^2\right]=C \end{align*} \]

여기서 한가지 생각해 볼 것은, \(y\) 는 \(x\) 의 함수 (\(y=y(x)\)) 이고, 다음과 같은 2개의 조건을 만족시켜야 한다. 즉, \(y(0)=0\) 이고 \(y(b)'=0\) 이다. \(y(0)=0\) 는 문제 설정상 원점에서 시작하기 때문에 만족시켜야 하며, \(y(b)'=0\) 는 calculus of variation 에서 variable end point 문제의 natural boundary condition 에 의해서 반드시 성립되어야 한다. 이를 analytic 하게 풀면 결과는 cycloid 가 된다. 즉,

\[ \begin{align*} x=\frac{1}{2}C(\theta-\sin \theta) \\ y=\frac{1}{2}C(1-\cos \theta) \end{align*} \]

여기에 위의 두 조건을 대입하면 다음과 같은 path 를 얻게 된다.

\[ \begin{align*} x(\pi)&=\frac{1}{2}C \pi = b \\ \therefore C&=\frac{2b}{\pi} \\ \implies x&=\frac{b}{\pi}(\theta-\sin \theta) \\ y&=\frac{b}{\pi}(1-\cos \theta) \\ &\textrm{where } \theta \in [0,\pi] \end{align*} \]

그런데, 이것을 수치적으로 풀려면 어떻게 해야 할까? 즉, 우리는 \(y[1+(y')^2]=C\) 라는 미분방정식 하나만 가지고 있고, solution 이 cycloid 라는 것을 모르는 상황에서 어떻게 수치적으로 그래프를 얻어낼 수 있을까?

먼저, 위의 두가지 조건을 살펴보자. \(y(0)=0\) 이되면 위의 미분방정식 \(y[1+(y')^2]=C\) 를 만족시킬 수 없다. 왜냐하면 어떤 수에 \(0\) 을 곱하면 항상 \(0\) 이 되기 때문이다. 위 방정식을 만족시키기 위한 유일한 방법은 \(y'=\infty\) 이 되면 된다. 그러므로, path 는 origin 에서 무한대의 기울기를 가질거라는 것을 예상할 수 있다 (위의 analytic solution 의 경우와 일치한다). 그리고, \(y'(b)=0\) 이므로, \(C=y(b)\) 임을 알 수 있다. 그렇다면, 비록 아직까지 path 를 알지는 못하지만, \(C\) 와 \(y(b)\) 를 사용하여서, 수치적 방법으로 문제를 풀때 solution 을 얻어내는 하나의 criterion 으로 사용하자. 즉, 여러가지 \(C\) 를 사용해서 미분방정식을 수치적으로 풀고, \(y(b)\) 와 \(C\) 가 가장 가까워지는 경우의 \(C\) 를 찾아내면 문제를 푼것으로 간주할 수 있다.

수치적으로 문제를 풀기위해서 Matlab 의 ode45 함수를 사용할 것이다. \(y\) 의 초기값은 \(0\) 임을 알고 있지만, \(y'\) 의 분모에 \(y\) 가 들어가 있기 때문에 수치적으로 문제를 풀때는 \(y=0\) 을 사용할 수가 없다. 그러므로 \(\texttt{eps}\) 를 \(y\) 의 초기값으로 둔다.

결과 그래프는 아래와 같다. 약갂의 오차가 발생하지만, 비슷한 결과를 얻을 수 있다. 이때 \(C\) 는 \(0.01\) 부터 \(0.001\) 의 간격으로 \(2\) 까지 변화시켰다. \(y(b)-C\) 의 오차값도 함께 plot 하였다.

수치적으로 찾은 \(C\) 는 \(1.301\) 이지만, 해석적으로 찾은 \(C\) 는 \(4/\pi=1.2732\) 이다. \(C=4/\pi\) 를 대입하여 다시 수치적으로 계산한 결과를 보니 그 결과가 거의 정확히 일치하는 것을 알 수 있다. \(0.001\) 의 간격이 정확한 값을 찾아내기에는 충분히 미세하지 않는 것이었나 보다.

function Brachistochrone()
	clear;clc
	global c
	dif=[];
	b=2;
	xspan=[0 b];
	y0=eps;
	cspan=0.01:0.001:2
	n=length(cspan);
	for i=1:n
		c=cspan(i);
		[x,y]=ode45(@myfun,xspan,y0);
		dif(i)=abs(y(end)-c);
	end

	[temp ind]=min(dif);
	c=cspan(ind);
	% c=4/pi;
	[x,y]=ode45(@myfun,xspan,y0);
	plot(cspan,dif);
	grid on;
	xlabel('c');
	ylabel('y(b)-c')
	figure;
	plot(x,-y);
	hold on;
	grid on;

	theta=0:0.01:pi;
	x=b/pi*(theta-sin(theta));
	y=b/pi*(1-cos(theta));
	plot(x,-y,'r')
	xlabel('x (m)');
	ylabel('-y (m)');
	title('Brachistochrone path')
	legend('numerical','analytic')

	function dy=myfun(t,y)
		global c
		dy=zeros(1);
		dy(1)=sqrt(c/y(1)-1);
	end
end