Amor Natura

물리학특강2 수업을 듣는 학생들에게 도움이 될 자료입니다.

python으로 주어진 미분방정식에서의 vector field와 phase portrait를 그릴 수 있는 코드입니다.

Mathematica의 Streamplot 기능까지는 안되지만 도움이 될겁니다.

6.1.8_variation_uploaded.py

Planetary motion

leap-frog.py

미분방정식을 컴퓨터를 이용해 수치적으로 계산해야 할 때가 있다. 이 때 미분방정식을 푸는 방법이 여러개 있는데 가장 쉬운 방법은 Euler method이다. 하지만 에러가 꽤 커서 실제 연구할 때에는 사용하지 않는 방법이기도 하다. 미분의 정의로부터 euler method를 간단히 유도할 수 있다. 미분의 정의에서 극한을 빼면, 

처럼 쓸 수 있다. 물리적인 상황을 떠올리기 위해 위치를 s(t), 속도를 v(t)로 나타내었다. 시간에 대한 위치의 변화량이 속도가 된다. 이 식을 다시 고쳐보면 다음과 같다.

다음 시간(t+dt)에서의 위치는 지금 시간의 위치에서 지금 시간의 속도만큼의 변화를 주어 예측하는 방법이 된다. 이게 바로 Euler method이다. 사실 dt에 대해 테일러 전개를 해보면, 

에서 최대 (dt)의 제곱항만큼의 오차가 생기게 된다. 

이제 Euler method를 프로그램으로 만들어보자. 그 전에 생각해볼 문제가 있다. 바로 euler method를 어떻게 반복시킬 것인가 하는 문제이다. 지금의 위치와 속도를 통해 다음 시간 스텝에서의 위치를 예측하는 것이 이 method의 목적임을 상기하자. 그럼 다음과 같이 순서를 짜야 한다.

1. 우선 현재 시간 t 에서의 위치 x를 보자. 

2. 그 뒤에 현재 시간 t에서의 속도 v(t)를 계산해야 한다. 

3. 그럼 현재 속도 v(t)와 x를  Euler method에 넣어 t+dt 때의 위치를 어림계산할 수 있다. (물론 오차는....ㄷㄷ) 

4. 이후에 이 때의 위치 x+dx 가 t+dt 때의 시간이라는걸 컴퓨터에게 알려준다.

이걸 반복하면 된다. 검은 화살표는 약간의 연산을 해야 하는 과정이다. (속도구하기, Euler method)

녹색은 그냥 컴퓨터에서 대입을 해주는 부분이다.

프로그래밍을 편하게 하기 위해서는 미리 자신이 원하는 dt(시간 간격)와 최종 시간을 정한 뒤 배열의 크기를 잡아준다. 배열의 크기는 (전체 시간 / 시간 간격) 으로 정할 수 있다. 예를들어 2초 동안을 0.1초로 나누어 보고 싶다면 2/(0.1) = 20개의 조각으로 시간이 나뉘게 된다.

이후 결과를 츨력할 때 편하게 하기 위해 초기값을 배열의 맨 앞에 넣어준다. 그럼 배열의 크기 역시 1개가 더 필요한데, 이를 위해 처음부터 배열의 크기를 (전체시간 / 시간간격) + 1 로 잡아주자. 위의 예시에서 배열의 크기는 21이 된다.

아래는 python으로 Euler method를 이용해 지구에서 자유낙하 운동을 만든 예제이다. 질량은 1이고 위로 30m/s 만큼 던졌을 때의 시간에 대한 모습을 보여준다. (NumPy, matplotlib 설치 필요)

계산을 해보면 7.71 ~ 7.72초 사이일 때 땅에 떨어진다고 하는데 등가속도 운동을 계산해보면 약 7.7105 초에 정확히 지면에 떨어진다. 

다음으로 Euler method의 오차를 생각해볼 수 있는데, 한 step에서 계산 결과는 (dt)^2 만큼 쌓이게 된다. 하지만 우리가 최종결과를 원래 이론값과 비교할 때의 오차는 조금 다르다. step마다 계산할 때 생기는 오차를 local error라고 부르며, 우리가 원하는 시간 T에서의 계산값과 이론값의 차이는 global error라고 부른다. 그래서 global error는 local error보다 더 크기 마련인데, 정확히는 (dt)^1 만큼이다. 왜일까? 한 step에서 (dt)^2 만큼 오차가 발생하는데 이 오차가 T/dt = N번 누적된다. 여기서 N은 결국 1/dt에 비례하고, global error는 local error가 N번 누적된 것이므로, 

만큼의 오차가 생긴다. 이를 확인해보자. dt를 바꾸어가면서 어떤 미분방정식을 풀어본 뒤에 log-scale로 global error 와 dt에 대한 그래프를 그려보면 된다. 분명 기울기가 1이 나온다. 위의 프로그램을 조금 고쳐서 확인해보자. T=8일 때의 위치는 정확히 -13.6m 이다.

결과는 다음과 같다. 잘 보면 기울기가 log-log에서 -1이다. 

Global Error를 찍을 때 가끔 결과가 엉뚱하게 나오는 때가 있다.(기울기가 아무리해도 -1인 경우...심지어는 Runge-kutta 4th를 써도!!) 이번에 알게 되었는데 Nt의 크기에서 +1을 하지 않았을 경우 정확히 T=8 일 때가 아닌 T = 8-dt 때의 결과로 에러를 계산하기 때문이다.

이렇게 Euler method를 살펴보았다. 이해하기도 쉽고 만들기도 쉽지만 에러는 우리가 준 dt만큼밖에 줄어들지 않기 때문에 긴 시간뒤의 결과를 실제와 다르게 내주게 된다. 그렇다고 dt를 아주 작게 쪼갤 경우 loop를 많이 돌게 되어 컴퓨터 파워를 많이 써야 하는 부담이 있다. 이런 이유 때문에 Euler method를 상당히 초기조건에 민감한 비선형 미분방정식을 수치적으로 푸는데 적합하지 않다. 이를 극복하기 위해 Improved Euler, Leap-frog, Runge-Kutta 4th order 등의 여러 알고리듬이 있다. 

6.1.8.py

Strogatz의 Nonlinear Dynamics 교과서에서 6장 연습문제에 보면 van der Pol oscillator가 나온다. 식은 다음과 같이 생겼다. 

이 미분방정식들의 Vector Field를 파이썬으로 그릴 수 있는데 최대한 matplotlib함수를 쓰지 않고(mgrid 등) 그리는 코드를 짰다. 거기에 덧붙여 임의의 점 하나를 잡고 trajectory를 시간에 대해 그리도록 하는 기능을 같이 넣었다. portrait 아래에 나오는 두개의 그래프는 x-t 그리고 y-t 그래프이다. 

trajectory를 따라가보면 시간이 조금 지난뒤에는 진동하는 모습을 보여주며 phase portrait에서도 볼 수 있다. 

Vector plot을 하게되면 trajectory가 뚜렷히 보이지 않는 단점이 있다. vector의 크기가 작을 경우 화살표가 상대적으로 작게 나와 잘 보이지 않기 때문인데, mathematica에서는 StreamPlot이라는 명령어로 trajectory들의 portrait를 살펴보기 쉽다. 하지만 matplotlib에서는 Stream을 plot해주는 함수가 없기 때문에 누군가 만들어 놓은 것이 있다. 아래 링크를 살펴보시면 된다.

http://stackoverflow.com/questions/8296617/how-to-plot-a-streamlines-when-i-know-u-and-v-components-of-velocitynumpy-2d

아래 코드에서간단히 1개의 trajectory를 그렸는데 이를 격자형태로 모두 계산해 그리면 Stream plot이 되는건 뻔한 이야기. -.- (하지만 이건 하고싶지 않아져서 패스)

프로그램을 실행시키면 아래와 같다. (이런걸 왜 넣고 있지.)

참! 예쁘게 나올 필요는 없어서 여백조정은 패스~(사실 이걸 알아보려고 또 20분남짓 소비하는게 싫어서...)

코드는 아래에~!

파이썬으로 벡터장을 그려야 하는 경우가 있다. 심심하면 지도를 놓고 위에 풍향과 풍속을 그려볼 수도 있다. (기상청에서 데이터 받아서 그려보셔도 됩니다...정 심심하시면요...)

메뉴얼과 구글링을 통해서 더 깊이있는 방법들을 배울 수 있지만 가장 간단한 코딩은 아래와 같다.

파이썬에서 import 해올 때 matplotlib와 numpy를 다음과 같이 불러온다. (누군가는.... 이렇게 했다,.)

그런데 설치된 패키지에 pylab은 없다. Numpy와 Matplotlib만 있을 뿐이다. pylab의 설명은 matplotlib 홈페이지에 되어있는데 긁어오면 다음과 같다.

pylab combines pyplot with numpy into a single namespace. This is convenient for interactive work, but for programming it is recommended that the namespaces be kept separate, e.g.:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 5, 0.1);
y = np.sin(x)
plt.plot(x, y)

사실 pylab은 Numpy와 Matplotlib의 pyplot이라는 모듈을 같이 불러오는 방법이다. 따라서 pylab을 부르면 Numpy를 import하지 않아도 알아서 가져오게 된다. 하지만 Namespace을 같이 쓰는 것이 나중에 문제될 수 있으므로 되도록이면 위의 방법처럼 따로 불러와 쓰라고 한다.

matplotlib는 pylab에서 불러오는 pyplot이외에도 다양한 모듈이 존재한다. 나머지 모듈에서 제공하는 수많은 기능들을 보려면 역시 공식 홈페이지를 방문하시라. ( http://matplotlib.org/py-modindex.html )

함수들의 자세한 arguments는 메뉴얼을 참고.

import matplotlib.pyplot (이외에도 많은 모듇들이 존재함)

1. 라벨붙이기.

pylab.xlabel('aaa')

pylab.ylabel('bbb')

2. 범위 정해주기

pylab.xlim(1, 100)

pylab.ylim(1, 100)

3. grid 보이기

pylab.grid('on' / 'off')

4. title 쓰기

pylab.title('aaaa')

메스메티카에서는 쉽게 할 수 있는데 WolframAlpha에서 Streamplot이나 vectorplot처럼 함수를 부르면 벡터장을 그려주지 않는다. 어쩌다 발견한건데 더 좋은방법이 있는지는 모르겠음.

1. http://www.wolframalpha.com 에 들어가서 vector plot 이라고만 검색한다.

2. 그럼 예제가 하나 나오는데 그 위에 보면 옵션박스들이 있다!

Vector field에 원하는 함수를 넣고 나머지에는 그릴 구간을 정해준다.

그리고 Enter를 치면 다시 그림을 그려준다!

이걸로 그려볼 수 있게 되었음ㅋ

plot.py

꽤 많은 데이터들을 gnuplot으로 plot을 하고 싶을 때가 있다. 사실 python에서 matplotlib와 numpy로 싹 그릴수도 있지만 matplotlib에서 지원해주는 옵션을 자유자재로 다루기가 어려운 경우가 다반사다. 그렇다고 plot을 안할수는 없는 노릇이고... 

그냥 gnuplot의 스크립트 파일(e.g. plot.plt)을 만들어 실행시키는 방법을 써보자.

python이 빠르고 쉽게 문자열을 다룰수 있다는 강점이 있어 이를 사용했다.

(결국 python 관련 글....)

이를 실행시키면 한 디렉토리 안에 있는 dat파일을 모두 읽어 자신이 지정한 방법대로 그림을 그려 png로 저장해준다. 이 프로그램을 그대로 쓰기엔 조심할 점이 있는데, 

원래 데이터파일이 위와같이 tab이 여러번 존재한다고 해보자. gnuplot은 그래도 착하게 우리가 눈으로 보는 것과 같이 행을 구분해준다. using 1:4 와 같은 식으로 옵션을 써주어도 우리가 의도한대로 잘 실행된다.

하지만 위의 프로그램에서는 Column의 개수를 tab의 수로 가정했기 때문에 중복된 tab(위의 예제에선 3과 4 사이)를 없애주어야 한다. 없애는 방법은 아래 링크를 참조!

http://amornatura.tistory.com/62

removetab.py

한 디렉토리 안에 있는 모든 데이터파일을 읽어서 중복된 tab을 없애고 다른 확장자로 저장하는 코드. python을 완전 잘하는 편이 아니어서 그냥 기초적인 방법으로 주욱주욱 만들었음...

Class가 list인지 str인지에 따라 쓰이는 함수(method)가 다르기 때문에 꽤나 헷갈린다. 자주 다루지 않다보니 벌어지는일! 아무튼 이렇게 하면 중복된 tab을 없앨 수 있다. 

주의!!

용량이 큰 데이터파일들의 tab을 지워야 할 경우 이 프로그램은 컴터를 넉다운시킬수 있다. 왜냐면 readlines()로 그냥 한번에 주욱 읽어서 리스트로 저장하기 떄문이다... 메모리관리는 전혀 신경쓰지 않았는데 지금 다루는 데이터파일들 용량이 작은편이어서이다.ㅋㅋㅋ