📂Julia

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