본문으로 바로가기

얼마전에 저는 뭐 아무도 쓸일은 없을것 같았지만 그래도 몇 안되는 저의 취밍이자 흥미있어 하는 것이라 블로깅했던 글이 하나 있는데요. 바로 

[The Robot/Prog.Lang.] - MATLAB에서 4차 Runge Kutta를 이용하여 1차 혹은 2차 미분방정식을 푸는 예제

였습니다. MATLAB의 멀쩡한 Simulink나 ODE 명령이 있음에도 불구하고 과감하게(ㅠㅠ) 글을 올렸죠..ㅎㅎㅎ. 뭐 아무튼 그리고 그 글에 대한 응용 예제로 또 하나 후속글을 올릴려고 합니다. 그 예제로는 유명한 진자를 올릴려고 하죠. 사실 저는 꽤 예전에 단순한 진자(pendulum)를 대상으로 연재도 진했었습니만~~^^

저렇게 생긴 아이가 진자라는 아이죠^^. 저 아이의 동역학을 유도했던것은 [바로가기]부터 시작합니다. 뭐 이 글은 그 과정에는 관심이 없으니.. 살짝 결과만 가져오겠습니다. 그 중에서도 모터는 없이 자유 운동을 한다고 보죠.^^

그러면 수식이 위와 같습니다.^^. 위 수식을 이전에 했던 Runge Kutta로 2차 미분방정식을 풀던 예제[바로가기]에 적용하여 코드를 짜야죠^^

h = 0.01;
t = 0;
x = 30*pi/180;
v = 0;

fm = 0.05;
m = 0.1;
l = 100*0.01;
J = 0.02;
g = 9.8;

func = @(tVal, xVal, vVal) -fm/(m*l^2+J)*vVal-m*g*l/(m*l^2+J)*xVal;

figure;

set(gca, 'Xlim', [-80 80], 'Ylim', [-30,130], 'box', 'off', 'XColor', [0.801,0.801,0.801], 'YColor', [0.801,0.801,0.801], 'Color', [0.801,0.801,0.801])
hold(gca, 'on')

goundLine = line([-100, 100],[100, 100],'Color','k','LineWidth',12);
pendulum = line([0, l*100*sin(x)], [100, 100-l*100*cos(x)], 'Color','b','LineWidth',12);

while(1)
    kx1 = v;
    kv1 = func( t, x, v );
    
    kx2 = v + h*kv1/2;
    kv2 = func( t + h/2, x + h*kx1/2, v + h*kv1/2 );
    
    kx3 = v + h*kv2/2;
    kv3 = func( t + h/2, x + h*kx2/2, v + h*kv2/2 );
    
    kx4 = v + h*kv3;
    kv4 = func( t + h, x + h*kx3, v + h*kv3 );
    
    dx = h*(kx1 + 2*kx2 + 2*kx3 + kx4)/6;
    dv = h*(kv1 + 2*kv2 + 2*kv3 + kv4)/6;
    
    x = x + dx;
    v = v + dv;
    
    t = t + h;
    
    updateX = [0,l*100*sin(x)];
    updateY = [100, 100-l*100*cos(x)];
    
    set(pendulum,'Xdata',updateX,'Ydata',updateY);
    
    drawnow;
end

바로 위의 코드입니다. 너무 간단해서 딱히 설명할 것도 없네요. while문 안에 있는것이 Runge Kutta를 풀고 있는 것이구요. 그로부터 구해진 각도를 업데이트(updatedX/Y)해서 그래프를 그리는데 반영하는 것이지요.

미적감각이 무쟈게 아쉬운 저로서는 그냥 저렇게 선으로만 표현할 뿐이지요^^. 흠흠흠...^^

별것이 없습니다만... 그래도 이전에 다룬 Runge Kutta를 예제를 통해 직접 본다는 것과 MATLAB에서 간편하게 그 결과를 애니메이션을 확인할 때 어떻게 할 것인지에 대해 다루었다고 말하면 좀 억지죠???^^ 뭐 여하튼 전 여기까지 이야기할려구요.. ㅎㅎ^^


댓글을 달아 주세요

  1. BlogIcon 핑구야 날자 2014.10.22 08:15 신고

    동영앙을 보니 프로그램의 흐름이 조금은 이해가 되는군요...

  2. BlogIcon 복돌이^^ 2014.10.22 09:46 신고

    흐미....어지러워지네요....어질어질..^^
    다녀갑니다. 행복한 하루 되세요

  3. BlogIcon 화이트세상 2014.10.22 11:28 신고

    외계언어를 보고 가네요.. ^^*

  4. BlogIcon 또웃음 2014.10.22 12:54 신고

    어려운 내용이지만 잘 보고 갑니다.

  5. BlogIcon 세상속에서 2014.10.22 17:07 신고

    어렵지만 잘보고 갑니다^^
    좋은 하루 되세요~

  6. BlogIcon 딸기향기 2014.10.22 21:13 신고

    물리 너무 어려워요 ㅠㅠㅠㅠㅠㅠㅠ 시험 어제 보고 온지라 아직도 그 여파갘ㅋㅋㅋ
    전 컴공인데 정말 프로그래밍이 여러곳에 사용되더라고요. 전자과 수업 들었을때도 쓰던데 물리에서도 쓰는군요.

  7. BlogIcon 스칼렛 오하라 2014.10.23 01:06 신고

    문외한은 그저 그려러니 한다는 ㅎㅎ 쏘리 ^^^

  8. BlogIcon 양군! 2014.10.23 06:26 신고

    유용한 정보에 대해서 정리 잘해주셔서 잘 보고 갑니다 ^^

  9. BlogIcon 워크뷰 2014.10.23 06:40 신고

    머리가 점점 기름칠이 되어 가네요^^

  10. 한동원 2014.11.06 16:12 신고

    이매트렙을 저 진자운동 시뮬레이션말고 t와 세타의 그래프를 표현하고 싶은데 코드를 어떻게 수정해야 가능할까요?ㅜ

    • BlogIcon PinkWink 2014.11.07 12:07 신고

      위에 있는 x,와 t 변수를
      plot(t,x)
      hold on
      일단 이렇게만 한번 해보세요. 그러면 좀 보기 안좋지만 뭔가가 나타날겁니다.
      좀 더 깨끗한 그래프는 t와 x를 행렬로 저장해서 한번에 뿌리는 것도 괜찮습니다. 시뮬레이션 시간이 몇 초로 짧으니까요^^

  11. 김진완 2014.11.09 17:30 신고

    ODE 명령은 어떻게 사용해야할지 알수 있을까요?
    처음 사용시 ode 정의되지 않은 함수라 나오네요

  12. BlogIcon 열공합시다 2014.12.04 21:36 신고

    저기 처음부분에서 변수들 중에 h, J, l, m, fm 이 각각 뭘 뜻하는 건지 좀 알려주세요ㅠㅠ

    • BlogIcon PinkWink 2014.12.04 22:59 신고

      h는 이전 runge kutta설명에서 이야기했지만 시간간격을 의미합니니다.
      또 본문에서 언급하는 진자의 동역학을 유도하는 글에서 언급이 있지만 (정리는 안되어 있고..ㅠㅠ)
      J는 관성모멘트
      l은 진자의 길이
      m은 진자의 무게
      fm는 마찰계수입니다.

  13. 치명적인남자 2016.10.24 23:00 신고

    혹시 ode45를 사용하여 x''+sin(x)=0, x(0)=0.5, x'(0)=0.5 를 구현 하려고 하는데요 ode45로 구현하는 방법도 올려주실수 있으신가요??

    • BlogIcon PinkWink 2016.10.25 08:49 신고

      ode는 (오래되어서 기억이 가물거리지만...ㅠㅠ) 선형미분방정식에서만 사용가능합니다. sin함수는 근사화하셔서 표현하시면 됩니다.