📂줄리아

줄리아 미분방정식 풀이 후 미분계수 얻는 법

개요

줄리아의 미분방정식 패키지인 DifferentialEquations.jl을 사용하여 상미분방정식을 풀 때, 내장된 솔버가 이미 사용한 미분계수를 별도 계산 없이 얻는 방법을 소개한다.

코드

로렌츠 어트랙터

function lorenz(u, p, t)
    x, y, z = u
    σ, ρ, β = p

    return [
        σ*(y - x),
        x*(ρ - z) - y,
        x*y - β*z
    ]
end

u0 = [1.0, 1.0, 1.0]
p = [10.0, 28.0, 8/3]
tspan = (0.0, 10.0)

prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob, RK4())

아주 유명한 예시로써 누구나 아는 로렌츠 어트랙터를 보자. 참고로 결과를 쉽게 보기 위해 룽게-쿠타 메소드 RK4()를 사용한다. OrdinaryDiffEqLowOrderRK에서 따로 불러오는 등 번거로운 점이 있는 건 어쩔 수 없다.

julia> usol = stack(sol.u)
3×238 Matrix{Float64}:
 1.0  1.0317    1.07399   1.19117   …  -6.83541  -5.84547  -5.14917  -4.97414
 1.0  1.41811   1.6532    2.10273      -3.67059  -3.51281  -3.73655  -3.89891
 1.0  0.977369  0.968804  0.961702     29.0521   27.1696   25.2991   24.6334

해 sol은 정상적으로 계산되었다. sol의 필드 u는 각 시간에서 점의 위치를 가지고 있다.

sol.k

한편 sol의 또 다른 필드인 k는 풀이 중에 사용된 미분계수들을 가지고 있다. 가령 RK4라면 각 스텝마다 시작점과 끝점의 미분계수를 계산하는데, 이에 따라서 두 개의 미분계수가 sol.k에 들어가 있다.

julia> dsol = [sol.k[2][1] stack([solk[2] for solk in sol.k[2:end]])]
3×238 Matrix{Float64}:
  0.0       3.86404   5.79208   …   31.6482   23.3266    14.1261   10.7524
 26.0      26.4613   27.3781        10.862    -1.34155  -10.1708  -12.847
 -1.66667  -1.14325  -0.807953     -52.3823  -51.9181   -48.2241  -46.2954

다시 말해, 미분계수는 솔버의 특징에 따라 얼마든지 다를 수 있다. 위의 복잡한 코드는 두 개의 미분계수 중 두번째 것을 참조하되, 첫번째 미분계수는 sol.k[2][1]로 첫번째 것을 따로 가져와서 합쳐놓은 것이다.

julia> dsol = stack(lorenz(u, p, t) for (u, t) in zip(sol.u, sol.t))
3×238 Matrix{Float64}:
  0.0       3.86404   5.79208   …   31.6482   23.3266    14.1261   10.7524
 26.0      26.4613   27.3781        10.862    -1.34155  -10.1708  -12.847
 -1.66667  -1.14325  -0.807953     -52.3823  -51.9181   -48.2241  -46.2954

실제 각 지점에서 lorenz 함수로 직접 계산한 미분계수를 보면 정확히 일치하는 것을 알 수 있다.

전체 코드

using DifferentialEquations
using OrdinaryDiffEqLowOrderRK: RK4

function lorenz(u, p, t)
    x, y, z = u
    σ, ρ, β = p

    return [
        σ*(y - x),
        x*(ρ - z) - y,
        x*y - β*z
    ]
end

u0 = [1.0, 1.0, 1.0]
p = [10.0, 28.0, 8/3]
tspan = (0.0, 10.0)

prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob, RK4())

usol = stack(sol.u)
dsol = [sol.k[2][1] stack([solk[2] for solk in sol.k[2:end]])]

dsol = stack(lorenz(u, p, t) for (u, t) in zip(sol.u, sol.t))