1. 개요

 상미분방정식(이하 상미방)의 수치 적분을 수행하기 위해서 scipy.integrate.odeint 함수의 사용법에 대해서 알아보도록 하겠다. scipy.integrate 클래스는 적분을 수행하는 다양한 함수들이 모여 있는데 그 중 연립 상미방의 수치적분을 구하는 ode, odeint, complex_ode 함수들이 마련되어 있다.

• odeint             : 상미방의 수치적분

• ode                : 상미방의 수치적분 (generic interface class)

• complex_ode  : 복소 상미방의 수치적분

이중 odeint 는 내부적으로 FORTRAN 라이브러리인 odepack 의 'lsoda'를 사용하여 상미방을 풀어내는 함수이다. 이 함수는 (연립) 상미분 방정식

• y' = f ( y, t )

의 수치해를 구해주며 stiff ode 와 non-stiff ode 모두에 적용된다.

2. 기본 문법

 일단 레퍼런스의 함수 문법은 다음과 같다.

여기서 필수적인 인수는 func, y0, t 세 개이며 각각 ODE함수, 초기값, 시간(벡터)이다. 나머지 인수들은 모두 선택적으로 지정해 줄 수 있으며 지정하지 않으면 설정값이 사용된다. 그리고 y0와 t는 리스트, 튜플, ndarray (이것들을 묶어서 array_like 라고 한다.) 등이 될 수 있다.

 이 함수를 사용하기 전에 먼저 상미방을 기술하는 함수를 먼저 정의해야 한다. 이 함수는 array_like 를 반환해야 하며 ndarray y와 시간 t 를 받아서 그 미분을 구하는 함수이다.

 odeint() 함수의 출력 y는 배열(ndarray) 이고 shape 이 ( len(t), len(y0) ) 이며 2차 배열이다. 입력으로 넘어온 t 배열의 각 점에서의 y값들을 갖는다. 두 번째 출력인 infodict 는 full_output== True 로 지정되었을 경우 추가적인 정보가 넘어오게 된다.

3. 간단한 1차 ode 예제

 다음과 같이 가장 간단한 1차 미방을 시간 구간 [0, 5] 에서 수치해를 구하는 예제를 해 보도록 하겠다.

• y' =  -2y, y(0)=1

위 예제에서는 pylab 모듈에 numpy가 포함되어 있으므로 따로 numpy는 import하지 않았다.

4. 2차 ode : 조화진동자

  조화 진동자 방정식은 2차 상미방으로서 x''(t) = -x(t) 이다. 2차 방정식은 두 개의 1차 방정식으로 분해할 수 있다. 즉, y(t) = [ x(t) ; x'(t) ] =: [y1(t); y2(t)] 라고 정의하면

              y'(t) = [x'(t); -x(t)] = [y2(t); -y1(t)]

와 같이 기술할 수 있다. 초기값이 [0; 1] 일 경우 시간 [0, 5]에 대해서 이 상미방을 풀면 다음과 같다. 이 예제에서는 ode1()함수의 반환값도 튜플이고 odeint()함수의 초기값도 튜플로 주었다.

 만약 x, y가 심볼로 등록이 되어 있다면 함수 y(x)의 x에 대한 미분은 sympy에서 다음과 같이 표현할 수 있다.

• y(x) : 독립 변수 x에 대한 함수 y

• y(x).diff() 혹은 y(x).diff(x,1) : 1차 미분

• y(x).diff(x,2) : 2차 미분

• y(x).diff(x,3) : 3차 미분

이것을 이용하면 상미분식을 표현할 수 있다. 예를 들면 다음과 같다.

[표 1] 미분식의 심파이 표현식 예

함수를 표현할 때 y(x)와 같이 독립변수도 반드시 표기해야 함을 유의하자.

 dsolve()함수는 미분방정식을 풀어주는 sympy의 기본 함수이다. 첫 번째 인자로 미분식을, 두 번째 인자로 구할 함수를 적어주면 된다. 예를 들어서 다음과 같은 1차 상미방을 푼다고 가정하자.

이 미분방정식의 우변을 0으로 만들면 다음과 같다.

이제 dsolve 함수의 첫 번째 인자로 위 식의 좌변을 넣어주고 두 번째 인자로 y(x)를 인가하면 된다.

여기서 C1은 상수이다.  다른 미분 방정식을 풀어보자.

이 방정식도 우변을 먼저 0으로 만들어야 한다.

이제 이 식의 좌변을 dsolve()함수의 첫 번째 인자로 넣어준다.

2차 미분방정식도 입력할 수 있다. 예를 들어

을 풀려면 다음과 같이 하면 된다.

만약 함수의 미분표현식을 매번 입력하기 번거롭다면 다음과 같이 새로 정의하여 사용하면 된다.

이렇게 정의해 놓으면 이후부터는 y(x).diff() 대신 dy를, y(x).diff(x,2) 와 같이 긴 표현식 대신 ddy를 사용하면된다.

여기서 C1과 C2는 상수이다. 이것이 dsolve()함수를 이용한 미분방정식을 푸는 가장 기본적인 방법이다.

 심파이에서 행렬객체에 수행할 수 있는 선형대수 관련 연산을 다음 표에 정리하엿다.

[표 1] 선형대수 연산

숫자로만 이루어진 행렬뿐만 아니라 대수 기호가 포함된 행렬에 대한 연산도 수행 가능하다.

만약 대수 기호 대신에 숫자를 입력하고 싶다면 subs() 메쏘드를 이용한다.

 행렬 간 산술 연산은 +, -, *, ** 연산자로 수행할 수 있다. *은 행렬 간 곱셈을, **은 거듭제곱을 수행한다.

[표 1] 행렬의 기본 연산

덧셈과 곱셈은 두 행렬의 크기가 같아야 하고 행렬 간 곱셈은 차수 조건에 맞아야 한다. 즉, 첫 번째 행렬의 열 수와 두 번재 행렬의 행 수가 같아야 곱셈이 성립한다.

 파이썬에서 **가 거듭제곱 연산자이므로 심파이에서도 행렬의 거듭제곱은 ** 연산자로 수행한다. 거듭제곱을 수행하려면 행렬이 정방행렬이어야 한다.

  행렬 객체에 대해서 파이썬 인덱싱과 슬라이싱도 사용할 수 있다. 주의할 점은 수학의 행렬은 1행, 1열부터 시작하지만 Matrix객체는 인덱스가 0으로 부터 시작한다는 점이다. 따라서 첫 번째 행(열)의 인덱스는 0이고 두 번째 행(열)의 인덱스는 1이다.

특정한 행이나 열을 뽑아내고 싶다면 row(n), col(n) 멤버 함수를 이용한다.

행렬 M의 복사본을 생성하고 싶다면 다음과 같이 하면 된다.

이와 같이 M2를 변경해도 M에는 영향을 미치지 않는다.

행렬의 특정 행이나 열을 삭제하고 싶을 경우는 row_del(), col_del() 멤버함수를 이용하면 된다.

반대로 특정한 행이나 열을 끼워 넣고 싶다면 row_insert(n), col_insert(n) 을 이용한다. 이 경우 기존 행렬의 요소들을 밀어내고 그 위치에 입력한 행이나 열을 끼워 넣게 된다.

 그리고 A.row_join(B) 는 A행렬의 오른편에 B행렬을 병합시키고 A.col_join(B)는 A행렬의 하단에 B행렬을 병합시킨다.

또한 행렬객체의 메쏘드 중에 extract()라는 함수가 있다. 이것을 이용해서 원하는 부분을 뽑아낼 수 있다.

예를 들면 다음과 같다.

이 예를 보면 [0행, 1행, 3행] 과 [0열, 1열]이 겹치는 부부만 뽑아져 나온 것을 알 수 있다. 해당 행과 열에 선을 그어서 선이 겹치는 부분을 확인해 보면 쉽게 이해할 수 있을 것이다.

 심파이를 가장 간단하게 실행해 볼 수 있는 방법은 live.sympy.org 에 접속하는 것이다. 이 페이지에 접속하면 심파이의 최신 버전이 임포트된 파이썬 쉘을 사용할 수 있다.

[그림 1] live.sympy.org의 접속 화면

맨 처음 접속했을 때 상단에 보면 다음과 같이 미리 실행되는 명령들이 있다.

 미리 실행되는 이  명령들에 의해서 심파이 모듈이 글로벌 환경에 임포트되어 있고 또한  변수 x, y, z, t 는 심파이의 대수 기호로 사용할 수 있다. 그리고 k, m, n은 정수 대수 기호 그리고 f, g, h 는 함수 기호로 바로 사용할 수 있다.

 만약 윈도를 사용하고 있고 로컬PC에서 심파이를 실행하고 싶다면 winpython을 사용하면 된다. 이 패키지는 파이썬과 필요한 모듈들이 대부분 포함되어 있으며 설치방식이 아니라 단순히 압축을 푼 후 실행하는 방식 (포터블)이라 사용하기 용이하다. 파이썬 2.x버전, 3.x 버전 각각 별도로 그리고 윈도 32bit, 64bit 버전별로 배포본이 마련되어 있다. 자기 환경에 맞는 것을 선택하여 다운로드하면 된다.

필자는 winpython 3.4.3.2 버전을 다운로드 받아서 실행하였다. 그러면 그 실행 파일의 하위 폴더에 다음과 같은 파일들이 생성되어 있다. 그 후 Spyder.exe를 실행하면 된다.

[그림 2] winpython의 설치 폴더

Spyder 가 실행된 후 Tools>Preferences 메뉴에서 다음 그림과 같은 항목을 체크한 후 ipython쉘을 실행시키면 바로 sympy를 사용할 수 있다.

[그림 3] winpython에서 sympy를 기본으로 실행시키기

이제는 로컬PC에서 ipython과 sympy를 사용할 수 있다.

[그림 4] winpython에서 ipython console을 실행시킨 화면

 다음 표는 산술/논리 연산에 관련된 내장 함수들이다.

[표 1] 산술/논리 연산에 관련된 내장 함수들

hex(), oct(), bin() 함수는 각각 ‘0x’, ‘0o’, ‘0b’ 로 시작하는 ‘문자열’로 결과를 반한한다. 이 문자열을 파이썬 값으로 변환하려면 eval()함수를 이용하면 된다.

 다음은 열거형(sequence)의 정보를 얻는 내장 함수들이다.

[표 1] 열거형의 정보를 얻는 내장 함수들

iter()은 seq형으로부터 이터레이터를 반환하는 함수이다. next()는 이터레이터의 현재 요소를 반환하고 이터레이터의 포인터를 그 다음 요소로 넘기는 함수이다.

enumerate()은 시퀀스로부터 enumarate형을 반환한다. 이것은 인덱스와 해당 값의 튜플을 반환하는 이터레이터를 가진다.

보통은 for 문과 같이 사용되어서 반복문 안에서 이터러블의 요소값 뿐만 아니라 그 인덱스 정보도 필요할 경우 편리하게 사용할 수 있다.

filter(func, iterable)는 함수와 이터러블을 입력으로 받아서 이터러블의 요소가 하나씩 함수에 인수로 전달될 때, 참을 반환시키는 것만을 따로 모아서 이터레이터로 반환하는 함수이다.

 다음 예를 살펴보자.

위의 positive함수는 리스트를 입력값으로 받아서 각각의 요소를 판별해 양수값만 따로 리스트에 모아 그 결과값을 돌려주는 함수이다.

 filter()함수를 이용하면 아래와 같이 똑같은 일을 동일하게 구현할 수 있다.

filter() 함수는 첫 번째 인수로 함수명을, 두 번째 인수로는 그 함수에 차례로 들어갈 이터러블(리스트, 터플, 문자열 등)을 받는다. filter 함수는 두 번째 인수인 이터러블의 각 요소들이 함수에 들어갔을 때 리턴값이 참인 것만을 이터레이터로 묶어서 돌려준다. 위의 예에서는 1, 2, 6 만이 양수로 x > 0 이라는 문장이 참이 되므로 [1, 2, 6]이라는 결과 값을 돌려주게 된다.

익명함수(lambda)를 쓰면 더욱 간편하게 쓸 수 있다.

 filter()와 유사한 map() 함수가 있다. 이것은 함수와 이터러블을 입력으로 받아서 각각의 요소가 함수의 입력으로 들어간 다음 나오는 출력값을 묶어서 이터레이터로 돌려주는 함수이다. 만약 어떠 리스트의 모든 요소의 제곱값을 갖는 새로운 리스트를 만들고 싶다면 다음과 같이 하면 된다.

파이썬 2.7은 filter()와 map()의 결과가 리스트이므로 위 예에서 list함수를 이용하여 리스트로 변환하지 않아도 된다. 파이썬 3.x에서는 이들 함수의 반환값이 리스트가 아닌 반복형(iterable)이므로 이것을 리스트로 만들기 위해서는 list()함수를 명시적으로 사용해야 한다.

  이전 포스트에서 객체를 생성하는 것이 외형상으로 함수를 호출하는 것과 유사하다고 했는데, 사실 객체를 생성할 때 클래스의 __init__() 메소드를 호출하게 된다. 이 메쏘드는 다음과 같이 사용자가 정의할 수 있다.

이 예에서 __init__(self) 메쏘드는 객체를 생성할 때 자동으로 호출되는 특수한 메소드이고 반드시 첫 번째 인자는 self 이어야 한다. 이 메쏘드 내부에 클래스 변수를 생성할 수 있다.

이제

이라고 Robot객체를 생성하면 nLegs 와 nArms 필드가 생성되고 초기화 되었음을 알 수 있다.

객체의 필드를 확인하고 싶다면 내장 함수 vars()를 이용하면 된다.

이 내장 함수는 사실 객체의 __dict__ 내부 필드를 반환한다. 이 내부 필드에 객체의 필드들이 딕셔너리의 요소로 저장되어 있다.

객체를 생성할 때 이름을 사용자가 입력하도록 하려면 다음과 같이 __init__() 메소드에 일반 인자를 self 뒤에 주면 된다.

이제 객체를 생성할 때 이름을 반드시 입력해야 한다.

만약 이름이 주어지지 않았을 때 ‘dummy’ 라는 이름을 주고 싶다면 다음과 같이 기본값 인자를 사용하면 된다.

이제 객체를 생성할 때 이름을 입력하지 않아도 된다.

이와 같이 __init__() 생성자 함수도 파이썬의 일반 함수처럼 가변 개수 인자, 기본값 인자, 키워드 인자 등을 적용할 수 있다.