logo

줄리아 미분방정식 패키지 DiffetentialEquations 튜토리얼 📂줄리아

줄리아 미분방정식 패키지 DiffetentialEquations 튜토리얼

설명

DifferentialEquations.jl는 SciML 그룹에 속하는 패키지 중 하나로, 미분 방정식의 수치적 풀이를 위해 개발되었다. 이 패키지로 풀 수 있는 방정식은 다음과 같다.

  • Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)
  • 상미분방정식(ODE)
  • Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)
  • 확률미분방정식(SDE)
  • Stochastic differential-algebraic equations (SDAEs)
  • Random differential equations (RODEs or RDEs)
  • Differential algebraic equations (DAEs)
  • Delay differential equations (DDEs)
  • Neutral, retarded, and algebraic delay differential equations (NDDEs, RDDEs, and DDAEs)
  • Stochastic delay differential equations (SDDEs)
  • Experimental support for stochastic neutral, retarded, and algebraic delay differential equations (SNDDEs, SRDDEs, and SDDAEs)
  • Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)
  • (Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)

튜토리얼 너머의 사용법은 아래를 참고하자.

설치

줄리아 REPL에서 다음과 같이 입력하여 설치한다.

julia> using Pkg
julia> Pkg.add("DifferentialEquations")

혹은 ]를 눌러 패키지모드에 진입하여 다음과 같이 입력하여 설치한다.

(@v1.10) pkg> add DifferentialEquations

ODEProblem

DifferentialEquations.jl에서 상미분방정식의 풀이는 다음과 같이 간단히 정리할 수 있다.

problem = ODEProblem(f, u0, tspan, p)
solution = solve(prob)

이 과정은 아래에서 구체적으로 다루도록하고, 우선은 ODEProblem에 대해 자세히 알아보자. 함수 ODEProblem은 이름 그대로 상미분방정식 문제를 정의하는 함수이다. 기본적으로 네 개의 인수를 입력받을 수 있는데, 하나씩 살펴보자.

  • f: SciMLBase.ODEFunction 타입의 함수이다. 상미분 방정식을 아래와 같이 정리했을 때, f(u,p,t)f(u, p, t)에 해당하는 부분을 정의하면 된다. dudt=f(u,p,t) \dfrac{du}{dt} = f(u, p, t) 이때 메모리 효율을 위해서 in-place 방식으로 정의할 수 있는데, 이 때는 f!(du, u, p, t)와 같이 정의하면 된다.
  • u0: 초기 조건이다. uu가 스칼라 함수이면 u0 = 1.0과 같이, 벡터함수이면 u0 = [1.0, 1.0, 1.0]과 같이 두면 된다.
  • tspan: 문제에서 정의된 uu의 도메인을 튜플로 입력한다. tspan = (t_start, t_end)
  • p: 솔루션 uu를 제외한 파라매터(외력 등)를 입력한다.

이제 아래와 같은 상미분방정식이 주어졌다고 하자.

du(t)dt=u(t)+sin⁡(t)t∈[0,1]u(0)=0.0 \begin{align*} \dfrac{du(t)}{dt} &= u(t) + \sin(t) \qquad t \in [0, 1] \\ u(0) &= 0.0 \end{align*}

그러면 다음과 같은 코드로 문제를 정의한다.

using DifferentialEquations

f(u, p, t) = u + p(t)
u0 = 0.0
tspan = (0.0, 1.0)
p(t) = sin(t)
prob = ODEProblem(f, u0, tspan, p)

prob를 solve함수의 인자로 넣으면 방정식의 솔루션을 얻을 수 있다. 이제 아래의 예시를 통해 더 구체적으로 알아보자.

상미분방정식의 풀이1

가장 쉬운 상미분 방정식의 예로 다음과 같은 지수 성장 방정식을 고려하자.

dudt=αu \dfrac{du}{dt} = \alpha u

이 방정식의 해는 어떤 상수 AA에 대해서 AeαtA e^{\alpha t}와 같음을 이미 알고 있지만, 예제를 위해 가장 간단한 형태의 미분방정식을 가져왔다. 구체적으로는 다음과 같은 초기값 문제를 고려하자.

du(t)dt=4u(t)t∈[0,1]u(0)=0.5 \begin{align*} \dfrac{du(t)}{dt} &= 4 u(t) \qquad t \in [0, 1] \\ u(0) &= 0.5 \end{align*}

우선 우변을 함수 f(u,p,t)f(u, p, t)로 정의한다. 여기서 pp는 파라매터parameter이고 보통은 각 항의 계수, 이 방정식에서는 α=2\alpha = 2를 의미한다. 그리고 초기값과 시간의 범위를 설정한다.

더 정확하게 설명하자면 uu를 업데이트 하는 규칙을 정해주는 것이다. dudu가 어떻게 표현되는지를 함수로 정의하고, 이를 ODEProblem의 인자로 넣는 것이다. dtdt는 입력하지 않으면 알아서 설정한다. 즉 미분방정식을 du=f(u,p,t)dtdu = f(u,p,t)dt 꼴로 고친 뒤 f(u,p,t)f(u,p,t) 부분을 정의해주면 된다는 것이다. 이와 같이 파라미터가 상수인 경우에는 pp를 명시적으로 함수 내부에서 사용해야할 필요는 없다.

using DifferentialEquations

f(u, p, t) =  4u
u0 = 0.5
tspan = (0.0, 1.0)

이들을 ODEProblem에 인자로 넣어 풀어야할 문제 problem를 만들 수 있다.

julia> problem = ODEProblem(f, u0, tspan)
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 0.5

이제 problem을 solve에 대입하면 미분방정식의 수치적 해를 얻을 수 있다. 솔버가 솔루션을 7차원 벡터로 구해주었다.

julia> @time solution = solve(problem)
  0.000074 seconds (192 allocations: 28.750 KiB)
retcode: Success
Interpolation: 3rd order Hermite
t: 7-element Vector{Float64}:
 0.0
 0.07581611843444613
 0.21771610504885602
 0.39336627353058856
 0.6101928843972066
 0.8638706427431251
 1.0
u: 7-element Vector{Float64}:
 0.5
 0.581866090044407
 0.7728154714821961
 1.0981042913641448
 1.6942469046842847
 2.8139612733277692
 3.6945246786578037

proertynames를 사용하면 solution의 속성을 확인할 수 있다.

julia> propertynames(solution)
(:u, :u_analytic, :errors, :t, :k, :discretes, :prob, :alg, :interp, :dense, :tslocation, :stats, :alg_choice, :retcode, :resid, :original, :saved_subsystem)

솔루션 uu와 도메인 tt를 각각 얻고 싶다면 다음과 같이 할 수 있다.

julia> solution.u
7-element Vector{Float64}:
 0.5
 0.581866090044407
 0.7728154714821961
 1.0981042913641448
 1.6942469046842847
 2.8139612733277692
 3.6945246786578037

julia> solution.t
7-element Vector{Float64}:
 0.0
 0.07581611843444613
 0.21771610504885602
 0.39336627353058856
 0.6101928843972066
 0.8638706427431251
 1.0

시각화

solution에 대한 레시피가 이미 정의되어있기 때문에, 간단히 plot에 대입하면 그래프를 그릴 수 있다. 잘보면 solution.u와 solution.t를 직접 대입하여 그렸을 때와 달리 자동으로 값을 보간interpolation해주는 것을 알 수 있다. 또한 x축이 tt라는 것도 자동으로 표시된다.

using Plots
p1 = plot(solution)
p2 = plot(solution.t, solution.u)
plot(p1, p2, layout = (2, 1))

실제 해인 u(t)=0.5e4tu(t) = 0.5 e^{4t}와 비교해보자.

plot(solution, label="solver", lw=4)
plot!(LinRange(0,1,100), t->0.5*exp(4t), label="exact", lw=3, ls=:dash, lc=:red, dpi=300)

세부설정

키워드인수를 직접 설정하여 솔버를 선택하거나, 허용오차를 설정할 수 있다. 가령 0.050.05 타입스텝 마다 솔루션을 계산하고 싶다면 saveat = 0.05save at 로 설정하면 된다. 자세한 것은 링크에서 확인할 수 있다.

# 솔버를 Tsit5()로 설정
# 매 타임스텝 0.05마다 솔루션을 저장.
# 허용 절대오차는 1e-8로 설정.

julia> solve(problem, Tsit5(), saveat=0.05, abstol=1e-8)
retcode: Success
Interpolation: 1st order linear
t: 21-element Vector{Float64}:
 0.0
 0.05
 0.1
 0.15
 0.2
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0.55
 0.6
 0.65
 0.7
 0.75
 0.8
 0.85
 0.9
 0.95
 1.0
u: 21-element Vector{Float64}:
  0.5
  0.6107013642059927
  0.7459122224392803
  0.9110596077813341
  1.1127694183818455
  1.3591415739400332
  1.6600546629913426
  2.0275986016325755
  2.4765152918979516
  3.0248023813153906
  3.6945244502370813
  4.512485296082109
  5.511506995093366
  6.731841509374868
  8.22227376044907
 10.042494301852722
 12.266093755920075
 14.981974632546809
 18.298554052093934
 22.34964366466662
 27.298573538274052

코드 전문

using DifferentialEquations

f(u, p, t) =  4u
u0 = 0.5
tspan = (0.0, 1.0)

problem = ODEProblem(f, u0, tspan)
@time solution = solve(problem)

using Plots
p1 = plot(solution, title="plot(solution)")
p2 = plot(solution.t, solution.u, title="plot(solution.t, solution.u)")
plot(p1, p2, layout = (2, 1), dpi=300)

plot(solution, label="solver", lw=4)
plot!(LinRange(0,1,100), t->0.5*exp(4t), label="exact", lw=3, ls=:dash, lc=:red, dpi=300)

solve(problem, Tsit5(), saveat=0.05, abstol=1e-8)

연립 상미분방정식의 풀이2

연립 1계 미분방정식의 예로 다음과 같은 로렌츠 방정식을 고려하자.

dxdt=−σx+σydydt=−xz+ρx−ydzdt=xy−βz \begin{align*} \dfrac{dx}{dt} &= -\sigma x + \sigma y \\ \dfrac{dy}{dt} &= -xz + \rho x - y \\ \dfrac{dz}{dt} &= xy - \beta z \end{align*}

연립 미분방정식을 푸는 것도 미분방정식을 푸는 것과 비슷하다. du=(dx,dy,dz)du = (dx, dy, dz)를 업데이트하는 규칙을 정의하고 이를 ODEProblem의 인자로 넣으면 끝이다. 파라미터를 σ=10\sigma = 10, β=8/3\beta = 8/3, ρ=28\rho = 28로 설정하고, du=f(u,p,t)dtdu = f(u,p,t)dt꼴로 고치자.

dx=(−10x+10y)dtdy=(−xz+28x−y)dtdz=(xy−83z)dt \begin{align*} dx &= (-10 x + 10 y)dt \\ dy &= (-xz + 28 x - y)dt \\ dz &= (xy - \frac{8}{3} z)dt \end{align*}

이제 우변의 식을 dtdt를 제외하고 함수로 정의한다. 이때 lorenz가 아니라 in-place 함수인 lorenz!로 정의하는 것은 성능을 위함이다2. 이 경우에는 함수의 첫번쨰 인자에 출력인 du를 추가해주면 된다. u=(x,y,z)u = (x,y,z), du=(dx,dy,dz)du = (dx,dy,dz)라 두면,

dx=10(y−x)dtdy=[x(28−z)−y]dtdz=(xy−83z)dt  ⟹  du1=10(u2−u1)du2=u1(28−u3)−u2du3=u1u2−83u3 \begin{align*} dx &= 10(y-x)dt \\ dy &= \left[x(28-z)- y\right]dt \\ dz &= (xy - \frac{8}{3} z)dt \end{align*} \implies \begin{align*} du_{1} &= 10(u_{2}-u_{1}) \\ du_{2} &= u_{1}(28-u_{3})- u_{2} \\ du_{3} &= u_{1}u_{2} - \frac{8}{3} u_{3} \end{align*} 

function lorenz!(du, u, p, t)
du[1] = 10(u[2] - u[1])
du[2] = u[1]*(28-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)u[3]
end

초기값을 u(0)=(1,1,1)u(0) = (1,1,1), 도메인을 [0,100][0, 100]으로 설정한 뒤 해를 계산한다. 솔루션을 플랏할 때는 인수 idx=(1,2,3)을 꼭 추가해주어야한다.

u0 = [1.0, 1.0, 1.0]
tspan = (0.0, 100)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob)

plot(sol, idxs=(1,2,3))

환경

  • OS: Windows11
  • Version: Julia 1.10.0, DifferentialEquations v7.15.0, Plots v1.40.3

  1. https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Example-1-:-Solving-Scalar-Equations ↩︎

  2. https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Example-2:-Solving-Systems-of-Equations ↩︎ ↩︎