https://ritog.github.io/posts/1st-order-DE-julia/1st_order_DE_julia.html

Posts by Rito Ghosh

  * TILs

  * About
  * Now
  * CV
  * Publications
  *  
  *  
  *  
  *  
  *  
  *  

 

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

[fractal-20]

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 langauge 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 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 smaller 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 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.

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 Equaion 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 olutions
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
behaviour 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
timespan 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
 
 # timespan
 tspan = (0.0, 5e-6)     # timespan 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
timespan 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,
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.

[human-gen-sticker]

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