Amor Natura

http://nosee.tistory.com/294

python3 부터는 ANSI 기준으로 작성된 파일만 읽을 수 있습니다. 그래서 UTF-8로 작성된 파일은 보통 방법으로 읽을 때 에러가 납니다. 다음과 같이 말이죠.

UnicodeDecodeError: 'cp949' codec can't decode bytes in position ~~~~

그래서 우리가 읽으려는 파일의 인코딩 형식을 변환해주어야 합니다.

1. 리눅스에서 변경하기

블로그에 iconv 에 대한 설명이 있습니다. (http://amornatura.tistory.com/42) 

그냥 UTF-8에서 cp949(ANSI)로 파일을 변환하면 되기 때문에 다음과 같이 입력해줍니다.

iconv -f UTF-8 -t CP949 -o rename.txt original.txt

간단합니다. 

2. 파이썬3 에서 그냥 읽어오기

사실 윈도우에서 작업하게 되면 다른 프로그램들을 써서 변환해야 하기 때문에 매우 번거롭습니다. 리눅스에서도 많은 양의 파일을 변환하려면 일이죠. (뭐...1분이면 되지만요...ㅎ)

파이썬3 에서 직접 UTF-8 문서를 읽어올 수 있습니다. 자세한 내용은 다음의 링크를 따라가보세요.

(http://www.evanjones.ca/python-utf8.html)

위의 사이트에서 설명해주는 내용을 따라하면, 

import codecs

fileObj = codecs.open( "read.txt", "r", "utf-8" )

u = fileObj.readlines()

for i in u :

print(i)

print("\n")

codecs 라는 라이브러리를 불러와서 파일의 형식을 지정해준 후 읽어옵니다. 읽을 때 read()로 하면 한글자씩, readlines()하면 단위줄로 읽어오죠. 

쉽죠? 

물리학특강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 등의 여러 알고리듬이 있다. 

덴마크의 유명한 물리학자이자 양자역학의 확립에 큰 공을 세운 닐스 보어(Niels Bohr)의 기압계 이야기가 있다. 이걸 처음 접한 것은 네이버 웹툰이었는데 아래 링크를 통해 볼 수 있다.

http://comic.naver.com/webtoon/detail.nhn?titleId=163295&no=5

여기서 등장하는 닐스 보어의 기압계 이야기가 무엇인지 궁금하실거다. 위의 웹툰을 보면 금방 알겠지만 간단히 말해서 기압계로 건물의 높이를 재는 문제를 푸는 이야기다. 높이에 따라 기압이 달라지므로 그 차이를 이용해 건물의 높이를 재는 것인데, 닐스 보어는 그냥 기압계에 줄을 달고 떨어뜨린 뒤 줄의 길이를 재면 된다고 했다고 한다. 물리학적인 답을 요구하자 자유낙하 식을 이용해 구해도 된다고 말했고 이후에 여러 답을 더 제시했다는 이야기다.

그럴듯하지 않은가? 사실 양자역학을 처음 접했을 1900년대 초에는 고전역학과 너무 다른 (심지어 말도 안되는!) 가정과 생각들을 통해 기존의 물리학에서 자리잡고 있던 개념들을 일부 포기해야 하는 상황이었다. 닐스보어는 주어진 고정관념에서 벗어나 양자역학의 기초를 다진 사람들 중 한명이었고, 창의적인 생각의 일화가 있었어도 이상하지 않았을 것이다.

그래서 기압계 이야기는 너무나도 유명했다.

하지만... 이건 모두다 뻥이었다. -.-

구글에서 겁색어를 입력하다 보면 Niels Bohr barometer 뒤에 true가 자동완성된다. 응??? 왜???

사실 이 이야기는 1958년에 <리더스 다이제스트>에 처음 소개되었고 Dr.Alexander calandra가 만들었다고 전해진다. 1961년에는 그가 쓴 교과서에 이 이야기가 실려있다. 정말 그럴듯한 이야기들로 만들어져 있지만 실제 이야기도 아니고 닐스 보어가 그런 것도 아니었다. (심지어 원래 이야기에 닐스 보어가 언급되는지도 확실치 않다...) 단지 교육이 학생들에게 단순한 공식과 기계적인 풀이를 연습시키는 것이 아니라 스스로 답을 찾아가는 과정을 연습하는 것이어야 한다는 교훈을 주기 위한 우화였던 것이다.  

중요한건 이야기이지 사람이 아니었다. ㅎㅎㅎ

자세한 이야기는 아래 링크 참고~

http://www.snopes.com/college/exam/barometer.asp

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분남짓 소비하는게 싫어서...)

코드는 아래에~!

물빨래만 있던 시절이 있었다. 19세기까지는 그랬다. 기름이 물에 녹지 않기에 기름때를 지우는 일은 쉽지 않았을것이다. 하지만 우연한 발견으로 물 대신 다른 것으로 빨래를 하는 방법이 만들어졌다. 19세기 중반, 프랑스의 염색공장 사장인 jean baptiste jolly는 식탁보에 엎질러진 가솔린을 보게 된다. 그의 하녀가 실수로 가솔린을 쏟은 것인데 신기하게도 염색된 식탁보에서 가솔린이 묻은 부분만 하얀색으로 돌아왔다. 가솔린이 염색된 성분을 녹인 후 증발한 것이다. 보통 사장님들이라면 하녀를 불러 혼을 냈겠지만, jean baptiste jolly는 이 방법이 때를 녹이는 방법이 될 수 있겠다는 생각을 했다. 그 후 이 방법으로 가솔린이나 등유로 사람들의 옷을 세탁해주는 사업을 했고 큰 돈을 벌게 되었단다.

우연히 일어나고 지가나는 일들에서 우리는 생각보다 많은 진실들을 놓치고 있을지도 모른다. 하지만 jean baptiste jolly는 우연히 쏟아진 가솔린 때문에 망가진 식탁보를 관찰하고 이유를 추적하면서 생각의 전환을 해냈다. 물로만 세탁하리란 법이 어디있나. 사실 가솔린이 묻어 망가진 직물들은 많이 있었을텐데도 이전까지 '드라이 크리닝'을 생각해내지 못했다. 우연은 항상 준비된 사람에게 기회가 되는 법이다. 집요한 관찰과 상상력 그리고 고정관념에서 벗어날 수 있는 힘이 우리가 우연이라는 박을 썰어 안의 보석을 발견하도록 만들어 줄 것이다.

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

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

파이썬에서 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 )