외력이 있는 상미분방정식의 수치적 풀이 (줄리아 미분방정식 패키지 DifferentialEquations) 📂줄리아

외력이 있는 상미분방정식의 수치적 풀이 (줄리아 미분방정식 패키지 DifferentialEquations)

설명

DifferentialEquations.jl는 미분방정식의 수치적 풀이를 위한 줄리아 패키지이다. 본 글에서는 DifferentialEquations.jl을 이용하여 외력^{forcing term}이 있는 상미분방정식의 수치적 해를 구하는 방법을 소개한다. 튜토리얼을 읽고오면 도움이 될 것이다.

외력이 있는 1계 상미분 방정식

아래와 같은 간단한 문제를 고려하자.

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

이를 적분인자법으로 풀면 다음과 같은 정확한 해를 얻을 수 있다.

$u(t) = -\dfrac{1}{2} \left( \sin(t) + \cos(t) - e^{-t} \right)$

이제 패키지를 사용하여 위 방정식의 수치해를 구해보자. 우선 방정식을 $\dfrac{du}{dt} = f(u, p, t)$ 꼴로 고친다. 그러면 $f(u, p, t) = u(t) + \sin(t)$ 이고 코드로는 다음과 같이 정의한다. 여기서 주의할 점은 솔루션 $u(t)$ 는 u와 같이 작성해야하고, 파라미터인 $p(t)$ 는 p(t)와 같이 작성해야한다는 것이다.

using DifferentialEquations
using Plots

f(u, p, t) = u + p(t)

성능을 위해서 in-place 방식으로 정의하려면 다음과 같이 작성한다.

function f!(du, u, p, t)
    du[1] = u[1] + p(t)
end

이제 초기값, 도메인, 외력 함수를 정의하자. f를 in-place 방식으로 정의했다면, 초기값도 벡터로 u0 = [0.0]과 같이 정의해야함에 주의하자.

u0 = 0.0 
tspan = (0.0, 1.0)
p(t) = sin(t)

이제 ODEProblem에 대입하여 문제를 정의하고, solve로 문제를 풀면 끝이다.

prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)

그래프를 그려서 확인해보면 다음과 같다.

plot(sol, lw = 5, label = "numerical solution")
plot!(LinRange(0, 1, 100), t -> -(1/2)*(sin(t) + cos(t) - exp(t)), lw=3, ls=:dash, lc=:red, label="exact solution")

코드 전문

using DifferentialEquations
using Plots

f(u, p, t) = u + p(t)

function g!(du, u, p, t)
    du[1] = u[1] + p(t)
end

u0 = 0.0
tspan = (0.0, 1.0)
p(t) = sin(t)

prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)

plot(sol, lw = 5, label = "numerical solution", dpi=300)
plot!(LinRange(0, 1, 100), t -> -(1/2)*(sin(t) + cos(t) - exp(t)), lw=3, ls=:dash, lc=:red, label="exact solution")

외력이 있는 2계 상미분 방정식

다음과 같은 강제 조화 진동자 문제를 생각하자.

$\begin{align*} \ddot{x} + 2\dot{x} + 0.5 x &= \cos(2t), \qquad t\in [0, 40] \\ x(0) &= 0 \\ \dot{x}(0) &= 1 \end{align*}$

이는 2계 선형 비동차 방정식이므로 해는 보조해 $\phi(t)$ 와 특수해 $X(t)$ 의 합으로 구할 수 있다. 계산해보면 보조해는 $\phi(t) = 0.656 e^{(-1+\sqrt{0.5})t}-0.532e^{(-1-\sqrt{0.5})t}$ 이고 특수해는 $X(t) = -\frac{14}{113}\cos(2t) + \frac{16}{113}\sin(2t)$ 이다. 따라서 정확한 해는 다음과 같다.

$\begin{align*} x(t) &= \phi(t) + X(t) \\ &= 0.656 e^{(-1+\sqrt{0.5})t}-0.532e^{(-1-\sqrt{0.5})t} -\frac{14}{113}\cos(2t) + \frac{16}{113}\sin(2t) \end{align*}$

이제 패키지를 써서 이 문제를 풀어보자. 여기서 $x_{1} = x$ , $x_{2} = \dot{x}$ 라고 두면 위 문제를 다음과 같은 연립 상미분 방정식으로 나타낼 수 있다.

$\begin{align*} \dot{x}_{1} &= \dot{x} = x_{2} \\ \dot{x}_{2} &= \ddot{x} = - 2\dot{x} - 0.5x + \cos(2t) = - 2x_{2} - 0.5x_{1} + \cos(2t) \end{align*}$

$\implies \begin{align*} \dot{x}_{1} &= x_{2} \\ \dot{x}_{2} &= - 2x_{2} - 0.5x_{1} + \cos(2t) \end{align*}$

위 수식은 다음과 같이 코드로 작성한다.

function forced_harmonic_oscillation!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -2x[2] -0.5x[1] + p(t)
end

물론 p(t)를 그냥 cos(2t)로 두고 풀어도 무방하다. 그 경우엔 ODEProblem의 인수로 p를 넣지 않으면 된다. 앞선 예제에서와 같이 초기값, 도메인, 외력 함수를 정의하고 다음의 코드로 수치해를 구할 수 있다.

x0 = [0.0, 1.0]
tspan = (0.0, 40.0)
p(t) = cos(2t)

prob2 = ODEProblem(forced_harmonic_oscillation!, x0, tspan, p)
sol2 = solve(prob2)

실제 해와 일치하는지 그려보면 다음과 같다.

x(t) = 0.6564172054223187exp((-1+sqrt(0.5))*t) -0.5325234001125843exp((-1-sqrt(0.5))*t)-(14/113)*cos(2t) + (16/113)*sin(2t)

plot(sol2, lw = 5, label = "numerical solution", dpi=300, idxs=1)
plot!(LinRange(0, 40, 100), t -> x(t), lw=3, ls=:dash, lc=:red, label="exact solution")

코드 전문

using DifferentialEquations
using Plots

function forced_harmonic_oscillation!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -2x[2] -0.5x[1] + cos(2t)
end

x0 = [0.0, 1.0]
tspan = (0.0, 40.0)
p(t) = cos(2t)

prob2 = ODEProblem(forced_harmonic_oscillation!, x0, tspan)
sol2 = solve(prob2)

x(t) = 0.6564172054223187exp((-1+sqrt(0.5))*t) -0.5325234001125843exp((-1-sqrt(0.5))*t)-(14/113)*cos(2t) + (16/113)*sin(2t)

plot(sol2, lw = 5, label = "numerical solution", dpi=300, idxs=1)
plot!(LinRange(0, 40, 100), t -> x(t), lw=3, ls=:dash, lc=:red, label="exact solution")

환경

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