Numerical solution of ordinary differential equations with external forces (Julia DifferentialEquations Package)
Explanation
DifferentialEquations.jl is a Julia package for numerical solutions of differential equations. In this article, we will introduce how to find a numerical solution of an ordinary differential equation with a forcing term using DifferentialEquations.jl. Reading the tutorial might be helpful.
First-order ODE with Forcing Term
Consider the following simple problem.
$$ \begin{align*} \dfrac{du(t)}{dt} &= u(t) + \sin(t) \qquad t \in [0, 1] \\ u(0) &= 0.0 \end{align*} $$
Solving this using the integrating factor method, we obtain the exact solution as follows.
$$ u(t) = -\dfrac{1}{2} \left( \sin(t) + \cos(t) - e^{-t} \right) $$
Now let’s find a numerical solution to the above equation using the package. First, rewrite the equation in the format of $\dfrac{du}{dt} = f(u, p, t)$. Then it becomes $f(u, p, t) = u(t) + \sin(t)$ and is defined in code as follows. Note that the solution $u(t)$ should be written as u
, and the parameter $p(t)$ should be written as p(t)
.
using DifferentialEquations
using Plots
f(u, p, t) = u + p(t)
For performance, you can define it in an in-place manner as follows.
function f!(du, u, p, t)
du[1] = u[1] + p(t)
end
Now, define the initial value, domain, and forcing function. If you defined f
in an in-place manner, be aware that the initial value should also be defined as a vector, like u0 = [0.0]
.
u0 = 0.0
tspan = (0.0, 1.0)
p(t) = sin(t)
Finally, define the problem using ODEProblem
and solve it with solve
.
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)
Checking the graph for verification, it appears as follows.
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")
Complete Code
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")
Second-order ODE with Forcing Term
Consider the following forced harmonic oscillator problem.
$$ \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*} $$
Since this is a second-order linear nonhomogeneous equation, the solution can be obtained as the sum of the complementary solution $\phi(t)$ and the particular solution $X(t)$. Upon calculation, the complementary solution is $\phi(t) = 0.656 e^{(-1+\sqrt{0.5})t}-0.532e^{(-1-\sqrt{0.5})t}$, and the particular solution is $X(t) = -\frac{14}{113}\cos(2t) + \frac{16}{113}\sin(2t)$. Therefore, the exact solution is as follows.
$$ \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*} $$
Now let’s solve this problem using the package. Letting $x_{1} = x$ and $x_{2} = \dot{x}$, the above problem can be represented as the following system of ODEs.
$$ \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*} $$
The above expressions can be coded as follows.
function forced_harmonic_oscillation!(dx, x, p, t)
dx[1] = x[2]
dx[2] = -2x[2] -0.5x[1] + p(t)
end
Of course, you can simply set p(t)
as cos(2t)
and solve. In that case, don’t include p
as an argument to ODEProblem
. Define the initial value, domain, and forcing function as in the previous example, and solve for the numerical solution using the following code.
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)
Plotting to check for correspondence with the actual solution gives the following.
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")
Complete Code
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")
Environment
- OS: Windows11
- Version: Julia 1.10.0, DifferentialEquations v7.15.0, Plots v1.40.3