맥시마(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 살레시오

댓글을 달아 주세요