Search Results for '2012/08'

ATOM Icon

1 POSTS

  1. 2012/08/29 Brachistochrone Problem by Pilwon Hur

Brachistochrone Problem

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

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

사용자 삽입 이미지

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

Solution: Conservation of mechanical energy 에 의해서 다음과 같은 식을 얻을 수 있다.
사용자 삽입 이미지
이 문제는 T 를 최소화 하는 문제이며, Calculus of Variation 의 전통적인 문제중 하나가 된다. 즉, Euler-Lagrange (EL) 방정식을 아래와 같이 풀면 된다.
사용자 삽입 이미지
그런데, Lagrange L 을 잘 살펴보면, x에 independent 함을 알 수 있고, 이 경우에는 Hamiltonian 이 x 에 대해 일정한 값을 갖게 된다. 이를 이용해서 계속 풀어보면 아래와 같다.
사용자 삽입 이미지
여기서 한가지 생각해 볼 것은, 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 가 된다. 즉,
사용자 삽입 이미지
여기에 위의 두 조건을 대입하면 다음과 같은 path 를 얻게 된다.
사용자 삽입 이미지

그런데, 이것을 수치적으로 풀려면 어떻게 해야 할까? 즉, 우리는y[1+(y’)2]=C 라는 미분방정식 하나만 가지고 있고, solution 이 cycloid 라는 것을 모르는 상황에서 어떻게 수치적으로 그래프를 얻어낼 수 있을까?
먼저, 위의 두가지 조건을 살펴보자. y(0)=0 이되면 위의 미분방정식 y[1+(y’)2]=C 를 만족시킬 수 없다. 왜냐하면 어떤 수에 0을 곱하면 항상 0이 되기 때문이다. 위 방정식을 만족시키기 위한 유일한 방법은 y’ =∞ 이 되면 된다. 그러므로, 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 을 사용할 수가 없다. 그러므로 eps 를 y의 초기값으로 둔다.
결과 그래프는 아래와 같다. 약갂의 오차가 발생하지만, 비슷한 결과를 얻을 수 있다. 이때 C 는 0.01 부터 0.001 의 갂격으로 2까지 변화시켰다. y(b)-C 의 오차값도 함께 plot 하였다.

사용자 삽입 이미지

사용자 삽입 이미지

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

사용자 삽입 이미지


사용자 삽입 이미지

mfile
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);

Posted by Pilwon Hur

2012/08/29 09:03 2012/08/29 09:03
Response
0 Trackbacks , 0 Comments
RSS :
http://pilwonhur.net/tc/pilwonhur/rss/response/64

Trackback URL : 이 글에는 트랙백을 보낼 수 없습니다