Solving First Order Differential Equations with Julia

Using DifferentialEquations.jl to solve Linear 1st Order DEs

math
julia
computational-math
Author

Ritobrata Ghosh

Published

March 3, 2025

Solving First Order Differential Equations with Julia

Image by Pete Linforth from Pixabay

Introduction

In this tutorial, we will learn about solving Differential Equations with Julia. Differential Equations are everywhere in Science and Engineering problems. And, being able to solve Differential Equations using computers is very convenient. We will use the DifferentialEquations.jl package for solving Differential Equations.

Who is this for

This is for people who are already familiar with Differential Equations from Mathematics, and who can code, preferably in Julia. Only the very basics of Julia language is required. If you are unfamiliar with Julia, and yet you have non-trivial experience with any of Python, C/C++, Ruby, JavaScript, golang, etc., you are fine, and you can follow along.

What we will do

In this tutorial, we will only focus on using the high level API of the DifferentialEquations.jl package, and we will learn, from scratch, how to translate linear first order Differential Equations to Julia, and how to solve them using the library. We will not learn about and code up numerical methods of solving Differential Equations such as Euler’s method, Runge-Kutta methods, etc. Neither we will talk about testing for accuracy, testing convergence of our algorithms. We will simply and plainly have some DEs in our familiar mathematical notation, and then translate them to a form Julia and the library can understand, and then we call high level APIs of the DifferentialEquations.jl library to solve our DEs. That’s it.

And for this post, we are limiting ourselves to First Order Differential Equations.

What are Differential Equations: a very quick recap

“Normal” Equations

Let’s take a step back and discuss normal equations first. Normal equations are equations involving one or more unknown roots. Like this:

\[x^2 - 3x - 18 = 0\]

And you solve this equation for two values of \(x\), since this equation is quadratic.

How do you create this equation? You formulate an equation from a real-world scenario. Here’s an example, from which you will get the above equation:

A farmer is designing a rectangular vegetable garden. The length of the garden is 3 meters longer than its width. The total area of the garden is 18 square meters. How long are the sides of this garden?

This is how the equation is created: \(x(x-3) = 18\), where \(x\) is the length of the longer side of the garden.

Differential Equations

In normal equations, we solve the equation to find out values of variables previously unknown to us. In case of Differential Equations, we solve to find out functions previously unknown to us.

For normal equations, we find out a relation of the unknown variable with something known to us, concretely.

In DEs, we know how an unknown functions changes with time, or another variable. We form our equation with this derivative of an unknown function, and it’s derivatives- one or more.

It might look like this:

\[\dfrac{dy}{dx} = -y\]

Differential Equations are used to model real-world happenings, and we aim to solve for that function.

Forming a Differential Equation is only possible when there is a function that changes with another variable, let’s say, time. It is not enough for the output of a function to differ with time, it is also necessary that the function itself changes with time.

Here are some examples:

  • An ideal candle burns at a constant rate, i.e. the rate of burning of wax is the same when the candle is in its full length, as well as when it is nearly finished. But, the growth rate of population is not the same when the population is low as well as when it is hight. Not only population itself increases as population increases, but the rate of increase of population also increases with growing population. (Here I am talking about human population in modern age- with the absence of a natural predator, exceptionally good medical facilities, etc.)

  • Another example I can think of is the water emission rate of a hole that is used to empty a water tank. When the water tank is nearly full, the rate of emission of water from the hole is very high, but when the tank is almost empty, the rate of water emission is much lower.

  • When you have just made coffee, and say it’s temperature is \(70 ^\circ C\), and your surrounding temperature is \(8 ^\circ C\), the time taken by the coffee to reach \(30 ^\circ C\) is much lower than the time it takes to reach \(8 ^\circ C\) from \(12 ^\circ C\). The process of cooling is slower when the temperature difference is lower compared to when the temperature difference is larger. We will deal with this situation in our second example.

These are the situations that you model using Differential Equations.

Numerical Solutions of Differential Equations

There are some situations where finding the exact solution to a Differential Equation is not possible, or is impractical. This is why we try to find approximate solution. In many situations, this is good enough. When an analytical solution is impossible, numerical solutions is all we have.

What we expect here

For the Differential Equations that we will implement here, and will solve numerically, we will see a graph that will predict the behavior of the unknown function.

First Example: Radioactive Decay

Radioactivity is the phenomena in which an unstable nucleus loses energy by radiation.

In radioactive decay, if a sample of radioactive nucleus is kept, through several kinds of radiation, the original sample reduces in mass. If there were \(n\) number of nuclei originally, after a certain amount of time, there would be \(p\) number of nuclei, where \(p < n\).

We know that the amount of decay is directly proportional to the number of nuclei present in a sample.

\[ \begin{align} -\dfrac{dN}{dt} &\propto N \\ \implies -\dfrac{dN}{dt} &= \lambda N \end{align} \]

The sign is negative, because, the rate decreases with decreasing number of nuclei.

We will reframe this equation so that we can implement the equation in Julia, and solve it using the solver.

\[ \begin{align} &-\dfrac{dN}{dt} \propto N \\ &\implies -\dfrac{dN}{dt} = \lambda N \\ &\implies \dfrac{dN}{N} = -\lambda dt \end{align} \]

Julia package DifferentialEquations.jl expects functions to be in this format-

f(du, u, p, t) = p * u

Where du is the first derivative, u is the function, t is the time-span for which we want the function’s output, and p is the parameter, or a set of parameters.

From here, we can frame du to be du[1] = -lambda * u[1].

Importing Packages

import Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")

using DifferentialEquations, Plots
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.11/Project.toml`
  No Changes to `~/.julia/environments/v1.11/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.11/Project.toml`
  No Changes to `~/.julia/environments/v1.11/Manifest.toml`

Function Definition

function radio_decay!(du, u, p, t)
    lambda = -p
    du[1] = lambda * u[1]
end
radio_decay! (generic function with 1 method)

Declaring Parameters

N = [60.0e18]           # initial condition

# parameter
lambda = 2.28e6         # in SI unit, s^{-1}
p = lambda

# time-span
tspan = (0.0, 5e-6)     # time-span is very small because the half-life is very small
(0.0, 5.0e-6)

Defining the ODE Problem

prob = ODEProblem(radio_decay!, N, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 5.0e-6)
u0: 1-element Vector{Float64}:
 6.0e19

Solving and Plotting the Solution

sol = solve(prob)
retcode: Success
Interpolation: 3rd order Hermite
t: 17-element Vector{Float64}:
 0.0
 4.3859649122807024e-7
 6.41342034530829e-7
 9.724902283236061e-7
 1.2462761785169967e-6
 1.5654877939700708e-6
 1.8726226804630439e-6
 2.196788921416817e-6
 2.5203638024729064e-6
 2.8511887164148543e-6
 3.183703238312926e-6
 3.519781092460714e-6
 3.857529202669002e-6
 4.19723875341719e-6
 4.538177668884617e-6
 4.880273689241893e-6
 5.0e-6
u: 17-element Vector{Vector{Float64}}:
 [6.0e19]
 [2.2085932679488193e19]
 [1.3911043910467977e19]
 [6.538685467576463e18]
 [3.5026563834049787e18]
 [1.6917664046442598e18]
 [8.399157113217204e17]
 [4.011209837867301e17]
 [1.9182281652161018e17]
 [9.022998274997682e16]
 [4.227947963649966e16]
 [1.9650908828375188e16]
 [9.098778059095644e15]
 [4.194139573993325e15]
 [1.927908419090074e15]
 [8.838633947171844e14]
 [6.727182427233061e14]
plot(sol, title="Radioactive Decay", xlabel="Time", ylabel="N")

Second Example: Newton’s Law of Cooling

As we have seen in an example before, we know that something cools much faster when the difference between the temperature of that object and the surrounding temperature is comparatively higher. When the difference is lower, it takes a higher amount of time to cool down.

The problem was quantified by Isaac Newton himself. Let’s understand it.

The temperature of an object, \(T\) varies with time, \(t\), and it depends on the difference of the temperature of the object and its surroundings- the ambient temperature \(T_a\).

\[ \begin{align} \dfrac{dT}{dt} &\propto (T - T_a) \\ \implies \dfrac{dT}{dt} &= k(T - T_a) \end{align} \]

where \(k\) is a constant of proportionality, the value of which depends on the material.

Now we will implement this DE in Julia, and solve it using DifferentialEquations.jl library.

Function Definition

function newtons_cooling!(du, u, p, t)
    k, T_a = p
    du[1] = k * (u[1] - T_a)
end
newtons_cooling! (generic function with 1 method)

Expected format-

f(du, u, p, t) = p * u

Where du is the first derivative, u is the function, t is the time-span for which we want the function’s output, and p is the parameter, or a set of parameters.

Here, we are using in-place updating to define du because it is memory efficient, and will be a crucial where there is a system of Differential Equations present.

Here we have two parameters- k and T_a, and they are related as \(T_1 = k \times (T_0 - T_a)\). That is what we implement on line 3 of the function definition.

Declaring Parameters

T = [80.0]                # Temperature of object
k = -0.06                 # k, the cooling constant, negative, because cooling
T_a = 20.0                # the ambient temperature
tspan = (0.0, 300.0)      # time-period for which we want to tract the values of the function, i.e. temperature of the object
p = (k, T_a)              # defining a tuple of parameters as expected by the library
(-0.06, 20.0)

Defining the ODE Problem

Here we define the ODEProblem as required by the solver package.

prob = ODEProblem(newtons_cooling!, T, tspan, p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 300.0)
u0: 1-element Vector{Float64}:
 80.0

Solving the ODE

sol = solve(prob)
retcode: Success
Interpolation: 3rd order Hermite
t: 19-element Vector{Float64}:
   0.0
   0.18593390178623953
   2.045272919648635
   6.383848011370185
  11.959325821307548
  18.73361986647917
  26.74893491369098
  35.99530039828558
  46.545554755847775
  58.51808266239741
  72.13894261696798
  87.73678833989489
 105.79306069340441
 126.99454945496298
 152.3564523878484
 183.42977158769716
 222.68158609043894
 273.9980369813687
 300.0
u: 19-element Vector{Vector{Float64}}:
 [80.0]
 [79.33435782064105]
 [73.07086984267411]
 [60.9075122329383]
 [49.276503335091164]
 [39.49832238666797]
 [32.054188527309826]
 [26.921535861341653]
 [23.675325765485024]
 [21.792012228502433]
 [20.79153270072668]
 [20.31058542718253]
 [20.10523457704901]
 [20.029607729780903]
 [20.006573373974568]
 [20.0011132329579]
 [20.00018200375497]
 [20.000082850038083]
 [20.00001776796253]
plot(sol,
    xlabel="Time, t",
    ylabel="Temperature, T",
    title="Temperature vs. Time")

We see very clearly how the value of the function \(T\) changes with time.

Conclusion

We have seen the very basics and intuitive explanation of First Order Differential Equations, and learned about how DifferentialEquations.jl solves DEs, and how to frame DEs according to the format expected by that library.

Acknowledgements

I learned a lot about the content of this tutorial from- some amazing teachers that I had in college, and some amazing textbooks like Mathematical Methods for Physics and Engineering by Riley, Hobson, Bence and Mathematical Methods for Physicists by Arfken, Weber, and Harris. I learned about the DE solver package from many content from the creator: Chris Rackauckas, and some example videos from this YouTube channel- doggo dot jl.

Discuss

If you have read this post, and found it interesting or edifying, please let me know. I would like that very much. If you have any criticism, suggestion, or want to tell me anything, just add a comment or let me know privately. Discuss this post on the Fediverse, Hacker News, or Twitter/X.

Changelog

This is an Open Source blog. Feel free to inspect diffs in the GitHub repo.

Cite this Article

@ONLINE {,
    author = "Ritobrata Ghosh",
    title  = "Solving First Order Differential Equations with Julia",
    month  = "mar",
    year   = "2025",
    url    = "https://ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html"
}

Reuse

CC-BY-NC-SA