맥시마(Maxima)로 미분방정식을 풀 수 있다. 이때 사용되는 대표적인 함수로서 ode2() 함수와 desolve() 함수가 있는데 이 중에서 ode2()함수를 이용하여 아래의 2계 선형 상미방을 풀어보자.

이것을 손으로 풀면 A4용지 한 페이지가 꽉 찰 정도로 풀이 과정이 꽤 복잡한 문제이다. ode2()함수의 첫번째 인자는 미분방정식을 입력해야 된다. wxMaxima에서는

y'(x) 은 'diff(y,x)
y''(x) 는 'diff(y,x,2)

로 입력한다. 그리고 ode2()의 두 번째 인자는 함수(이 경우는 y), 세 번째 인자는 독립변수(이 예제에서는 x이다.)를 넣어주면 된다.


위에서 %k1, %k2는 적분상수를 표시한다. 초기 조건(initial condition)을 입력하여 적분상수를 구하려면 다음과 같이 ic2()함수를 이용하면 된다.


이것으로부터 미방의 해는 y=sin(x) + 3x/25 - 8/25 + 3.1e-0.5x 임을 알 수 있다.


 이 해의 그래프를 그려보고 싶다면 wxplot2d()함수를 이용하면 된다.

위 에서 wxplot2d()의 첫 번째 인자로 도시하려는 함수를 입력해야 하는데 rhs(%o2)라고 되어 있다. 이것은 결과가 저장된 변수 %o2 의 우변(right-hand side)을 추출하는 함수이다. 그리고 두 번째 인자는 [x,0,20] 인데 x축의 범위를 설정하는 것이다.



Posted by 살레시오
,

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. 기본 문법


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


y, infodic = scipy.integrate.odeint (
   func,# callable(y, t, ...) : 시간 t 에서 y(t)의 미분값을 구할 수 있는 함수
   y0,  # array_like : y 의 초기 조건 (벡터도 된다.)
   t,   # array_like : y(t)값을 구할 t 점들. 초기값의 시간이 이 배열의 첫 요소여야 함.
   args=()  # 튜플 : func() 에 넘길 추가적인 인수들
   Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None,
   tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
   mxords=5, printmessg=0
)

 

여기서 필수적인 인수는 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()함수의 초기값도 튜플로 주었다.


Posted by 살레시오
,