상미분방정식

상미분방정식 (Ordinary Differential Equations)

lsode는 상미분방정식(ODE)을 푸는 Octave의 대표 함수예요. 형태와 옵션, 그리고 실제 예시까지 하나씩 살펴볼게요.

출처: 문서

본문

lsode 함수는 다음 형태의 ODE를 풀 수 있어요.

dx
-- = f (x, t)
dt

내부적으로는 Hindmarsh의 ODE 솔버 LSODE를 사용해요.

[x, istate, msg] = lsode (fcn, x_0, t)
[x, istate, msg] = lsode (fcn, x_0, t, t_crit)

lsode

상미분방정식(ODE) 솔버예요. 푸는 미분방정식 집합은

dx
-- = f (x, t)
dt

이며 초기 조건은 x(t_0) = x_0이에요.

해는 행렬 x로 반환되는데, 각 행이 벡터 t의 한 원소에 대응해요. t의 첫 원소는 t_0이어야 하고 시스템의 초기 상태 x_0에 대응해야 해서, 출력의 첫 행은 x_0가 돼요.

첫 번째 인자 fcn은 방정식 집합의 오른쪽 변(right hand side) 벡터를 계산하는 함수 f를 가리키는 문자열, inline, 또는 함수 핸들이에요. 이 함수는 다음 형태여야 해요.

xdot = f (x, t)

여기서 xdot과 x는 벡터이고 t는 스칼라예요.

fcn이 두 원소짜리 문자열 배열 또는 두 원소짜리 셀 배열(문자열·inline·함수 핸들)이라면, 첫 번째 원소는 위에서 설명한 함수 f를, 두 번째 원소는 f의 야코비안(Jacobian)을 계산하는 함수를 가리켜요. 야코비안 함수는 다음 형태여야 해요.

jac = j (x, t)

여기서 jac은 편미분 행렬

             | df_1  df_1       df_1 |
             | ----  ----  ...  ---- |
             | dx_1  dx_2       dx_N |
             |                       |
             | df_2  df_2       df_2 |
             | ----  ----  ...  ---- |
      df_i   | dx_1  dx_2       dx_N |
jac = ---- = |                       |
      dx_j   |  .    .     .    .    |
             |  .    .      .   .    |
             |  .    .       .  .    |
             |                       |
             | df_M  df_M       df_M |
             | ----  ----  ...  ---- |
             | dx_1  dx_2       dx_N |

예요.

두 번째 인자는 시스템의 초기 상태 x_0를 지정해요. 세 번째 인자는 해를 구하려는 시간 값을 담은 벡터 t예요.

네 번째 인자는 선택 사항으로, ODE 솔버가 적분을 넘어가지 말아야 할 시간 집합을 지정해요. 특이점(singularity)이나 도함수에서 불연속이 있는 지점의 어려움을 피하는 데 유용해요.

계산이 성공적으로 끝나면 istate 값은 2가 돼요(Fortran 버전의 LSODE와 일치해요). 계산이 실패하면 istate는 2가 아닌 다른 값이 되고 msg에 추가 정보가 담겨요.

lsode_options 함수로 lsode의 선택 매개변수를 설정할 수 있어요.

lsode의 내부 동작에 대한 더 자세한 내용은 Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE Solvers", Scientific Computing, R. S. Stepleman, editor, 1983. 또는 https://computing.llnl.gov/projects/odepack을 참고하세요.

예시: 반 데르 폴(van der Pol) 방정식 풀기

fvdp = @(y,t) [y(2); (1 - y(1)^2) * y(2) - y(1)];
t = linspace (0, 20, 100);
y = lsode (fvdp, [2; 0], t);

See also: daspk, dassl, dasrt.

lsode_options

lsode 함수의 옵션을 조회하거나 설정해요.

lsode_options ()
val = lsode_options (opt)
lsode_options (opt, val)

인자 없이 호출하면 사용 가능한 모든 옵션의 이름과 현재 값이 표시돼요. 인자를 하나 주면 해당 옵션 opt의 값을 반환하고, 두 인자를 주면 lsode_options는 옵션 opt를 값 val로 설정해요.

옵션은 다음과 같아요.

  • "absolute tolerance" — 절대 허용 오차. 벡터 또는 스칼라일 수 있으며, 벡터라면 상태 벡터의 차원과 일치해야 해요.
  • "relative tolerance" — 상대 허용 오차 매개변수. 절대 허용 오차와 달리 스칼라만 가능해요. 각 적분 단계에서 적용되는 국소 오차 검사는 
      abs (local error in x(i)) <= ...
      rtol * abs (y(i)) + atol(i)
    이에요.
  • "integration method" — ODE 시스템을 푸는 데 사용할 적분 방법을 지정하는 문자열.
    • "adams", "non-stiff" — 야코비안을 사용하지 않아요(가능해도요).
    • "bdf", "stiff" — 강성 역차분 공식(BDF) 방법을 사용해요. 야코비안 계산 함수를 제공하지 않으면 lsode가 야코비안 행렬의 유한차분 근사를 계산해요.
  • "initial step size" — 첫 번째 단계에서 시도할 단계 크기(기본값은 자동 결정).
  • "maximum order" — 해법의 최대 차수를 제한해요. Adams 방법을 쓰면 이 옵션은 1과 12 사이여야 해요. 그렇지 않으면 1과 5 사이(포함)여야 해요.
  • "maximum step size" — 최대 단계 크기를 설정해 매우 넓은 영역을 지나쳐 넘어가는 것을 피해요(기본값은 미지정).
  • "minimum step size" — 허용되는 최소 절대 단계 크기(기본값 0).
  • "step limit" — 허용되는 최대 단계 수(기본값 100000).
  • "jacobian type" — 강성 역차분 공식(BDF) 적분 방법에서 사용할 야코비안의 종류를 지정하는 문자열.
    • "full" — 기본값. 모든 편미분을 근사하거나 사용자가 제공한 야코비안 함수에서 가져와요.
    • "banded" — 대각선과, 옵션 "lower jacobian subdiagonals"·"upper jacobian subdiagonals"로 지정한 아래·위 부대각선 수만 근사하거나 사용자 야코비안 함수에서 가져와요. 사용자가 제공한 야코비안 함수는 나머지 편미분을 임의의 값으로 설정해도 돼요.
    • "diagonal" — 사용자가 야코비안 함수를 제공하면 이 설정은 효과가 없어요. lsode가 근사하는 야코비안은 대각선으로 제한되는데, 각 편미분을 상태의 모든 원소에 유한 변화를 적용해 계산해요. 실제 야코비안이 정말 항상 대각선이라면 이는 상태의 각 원소에만 유한 변화를 적용하는 것과 같은 효과를 주면서 더 효율적이에요.
  • "lower jacobian subdiagonals" — 옵션 "jacobian type"이 "banded"일 때 사용하는 아래 부대각선 수. 기본값 0.
  • "upper jacobian subdiagonals" — 옵션 "jacobian type"이 "banded"일 때 사용하는 위 부대각선 수. 기본값 0.

lsode를 이용한 예시

lsode로 세 개의 미분방정식 집합을 푸는 예시를 볼게요. 함수가 다음과 같이 주어졌을 때,

## oregonator differential equation
function xdot = f (x, t)

  xdot = zeros (3,1);

  xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ...
            - 8.375e-06*x(1)^2);
  xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
  xdot(3) = 0.161*(x(1) - x(3));

endfunction

초기 조건이 x0 = [ 4; 1.1; 4 ]라면, 방정식 집합은 다음 명령으로 적분해요.

t = linspace (0, 500, 1000);

y = lsode ("f", x0, t);

이걸 직접 실행해 보면 t = 0과 5 사이, 그리고 t = 305 부근에서 결과 값이 극적으로 변하는 걸 볼 수 있어요. 더 효율적인 출력 지점 집합은 다음과 같아요.

t = [0, logspace(-1, log10(303), 150), ...
        logspace(log10(304), log10(500), 150)];

위에서 사용한 미분방정식의 m-파일은 Octave 배포판의 examples 디렉터리에 oregonator.m이라는 이름으로 포함되어 있어요.

더 알아보기

  • daspk, dassl, dasrt는 미분-대수 방정식(Differential-Algebraic Equations) 항목을 참고하세요.
  • ODE의 수치적 안정성과 강성(stiff) 문제도 함께 보면 좋아요.