Solving Differential Equations in Julia

Julia provides a powerful and efficient way to solve differential equations using the DifferentialEquations.jl package. This tutorial will walk through solving and plotting a simple ordinary differential equation (ODE).

To begin, install the necessary packages if you haven’t already:

using Pkg
Pkg.add(["DifferentialEquations", "Plots"])

Now, consider the first-order differential equation:

\[ dy/dt = -2y \]

with the initial condition y(0) = 1. The form of this differential equation is:

\[ \frac{dy}{dt} = f(y, p, t) \]

The DifferentialEquations.jl package, expects the ODE to be defined in terms of a Julia function that takes the arguments y, p and t.

The following Julia code sets up and solves this equation:

using DifferentialEquations, Plots

function f!(dy, y, p, t)
    dy[1] = -2y[1]  # dy/dt = -2y
end

y0 = [1.0]
tspan = (0.0, 5.0)

prob = ODEProblem(f!, y0, tspan)
sol = solve(prob)

plot(sol, xlabel="Time (t)", ylabel="y(t)", lw=2)

This will generate a plot of the exponential decay function y(t) = e^(-2t).

Note also the indexing with [1] in dy[1] = -2y[1]. Julia’s DifferentialEquations.jl package requires state variables (y) and their derivatives (dy) to be stored in array-like structures rather than scalars. Even for a single equation, they are treated as vectors, allowing the same framework to handle systems of equations.

For example, if we had a system of two equations:

\[ dy1/dt = -2y1 + y2\\ dy2/dt = -y1 - y2 \]

we would define:

function system!(dy, y, p, t)
    dy[1] = -2y[1] + y[2]
    dy[2] = -y[1] - y[2]
end

Even though we did not use the parameters p and time t in our simple example, they must still be included in the function definition. This is because DifferentialEquations.jl expects functions with the signature:

f!(dy, y, p, t)

Keeping p allows for easy parameterization of equations. For example, if the decay rate were variable:

function f!(dy, y, p, t)
    dy[1] = -p * y[1]
end

you could pass p when defining the problem:

prob = ODEProblem(f!, [1.0], (0.0, 5.0), 2.0)  # p = 2.0

Similarly, t is useful for time-dependent equations, such as:

function f!(dy, y, p, t)
    dy[1] = -2 * y[1] + sin(t)
end

These conventions make it easy to extend the model for more complex scenarios while maintaining compatibility with Julia’s solver ecosystem.