Differential equations with Julia

# Who am I ?

* My name is *Pierre Navaro*
  * **Fortran 77 + PVM** : during my PhD 1998-2002 (Université du Havre)
  * **Fortran 90-2003 + OpenMP-MPI** : Engineer in Strasbourg (2003-2015) at IRMA
  * **Numpy + Cython, R + Rcpp** : Engineer in Rennes (2015-now) at IRMAR
  * **Julia v1.0** since July 2018

## Instructions to open the notebook and slides

* https://github.com/pnavaro/JuliaSMAI2021
  * https://pnavaro.github.io/JuliaSMAI2021

---

# Why Julia?

* Started in 2009 and first version was released in 2012.
  * High-level languages like Python and R let one explore and experiment rapidly, but can run slow.
  * Low-level languages like Fortran/C++ tend to take longer to develop, but run fast.
  * This is sometimes called the "two language problem" and is something the Julia developers set out to eliminate.
  * Julia's promise is to provide a "best of both worlds" experience for programmers who need to develop novel algorithms and bring them into production environments with minimal effort.

**Julia: A Fresh Approach to Numerical Computing**

*Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah*

SIAM Rev., 59(1), 65–98. (34 pages) 2012

**SciML Scientific Machine Learning Software** https://sciml.ai/citing/

**Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia,**

*Christopher Rackauckas and Qing Nie*,

Journal of Open Research Software,, volume 5, number 1, 2017.

---

# Implement your own numerical methods to solve

$$ y'(t) = 1 - y(t),  t \in [0,5],  y(0) = 0. $$

## Explicit Euler

```julia
euler(f, t, y, h) = t + h, y + h * f(t, y)
```

```
euler (generic function with 1 method)
```

## Runge-Kutta 2nd order

```julia
rk2(f, t, y, h) = begin
    ỹ = y + h / 2 * f(t, y)
    t + h, y + h * f(t + h / 2, ỹ)
end
```

```
rk2 (generic function with 1 method)
```

---

## Runge-Kutta 4th order

```julia
function rk4(f, t, y, dt)

y₁ = dt * f(t, y)
    y₂ = dt * f(t + dt / 2, y + y₁ / 2)
    y₃ = dt * f(t + dt / 2, y + y₂ / 2)
    y₄ = dt * f(t + dt, y + y₃)

t + dt, y + (y₁ + 2 * y₂ + 2 * y₃ + y₄) / 6

end
```

```
rk4 (generic function with 1 method)
```

---

## Solve function

```julia
function dsolve(f, method, t₀, y₀, h, nsteps)

t = zeros(Float64, nsteps)
    y = similar(t)

t[1] = t₀
    y[1] = y₀

for i = 2:nsteps
        t[i], y[i] = method(f, t[i-1], y[i-1], h)
    end

t, y

end
```

```
dsolve (generic function with 1 method)
```

---

## Plot solutions

```julia
using Plots

nsteps, tfinal = 7, 5.0
t₀, x₀ = 0.0, 0.0
dt = tfinal / (nsteps - 1)
f(t, x) = 1 - x

t, y_euler = dsolve(f, euler, t₀, x₀, dt, nsteps)

t, y_rk2 = dsolve(f, rk2, t₀, x₀, dt, nsteps)

t, y_rk4 = dsolve(f, rk4, t₀, x₀, dt, nsteps)
```

```
([0.0, 0.8333333333333334, 1.6666666666666667, 2.5, 3.3333333333333335, 4.166666666666667, 5.0], [0.0, 0.5624678497942387, 0.8085656175363232, 0.9162413030173554, 0.9633528772107507, 0.9839657055671652, 0.9929844806797695])
```

---

```julia
plot(t, y_euler; marker = :o, label = "Euler")
plot!(t, y_rk2; marker = :d, label = "RK2")
plot!(t, y_rk4; marker = :p, label = "RK4")
plot!(t -> 1 - exp(-t); line = 3, label = "true solution")
```

```
"/private/var/folders/56/lhgl67bx4qj6tvjx41b1wgdr0000gn/T/jl_BAvc0r/build/dsolve.png"
```

![dsolve](dsolve.png)

---

## DifferentialEquations.jl

```julia
using DifferentialEquations

f(y, p, t) = 1.0 - y
y₀, t = 0.0, (0.0, 5.0)

prob = ODEProblem(f, y₀, t)

sol_euler = solve(prob, Euler(), dt = 1.0)
sol = solve(prob)
```

```
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 15-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0010999999999999998
 0.011099999999999997
 0.076742091466188
 0.2256219976613778
 0.4455760022489606
 0.7272505097568545
 1.0899643857570624
 1.5331652827982714
 2.06972399829104
 2.705750446504391
 3.4562444530999095
 4.3378403240538095
 5.0
u: 15-element Vector{Float64}:
 0.0
 9.999500016666247e-5
 0.001099395221772342
 0.011038622307372232
 0.07387132069817183
 0.20198030129465658
 0.35954472556518147
 0.5167641452122544
 0.6637713825376198
 0.7841482271226826
 0.8737784044472122
 0.9331779821204208
 0.9684489616227439
 0.9869310637824636
 0.9932596937939658
```

---

```julia
plot(sol_euler, label = "Euler")
plot!(sol, label = "default")
plot!(1:0.1:5, t -> 1.0 - exp(-t), lw = 3, ls = :dash, label = "True Solution!")
```

![diffeq](diffeq.png)

---

`sol.t` is the array of time points that the solution was saved at

```julia
sol.t
```

```
15-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0010999999999999998
 0.011099999999999997
 0.076742091466188
 0.2256219976613778
 0.4455760022489606
 0.7272505097568545
 1.0899643857570624
 1.5331652827982714
 2.06972399829104
 2.705750446504391
 3.4562444530999095
 4.3378403240538095
 5.0
```

`sol.u` is the array of solution values

```julia
sol.u
```

```
15-element Vector{Float64}:
 0.0
 9.999500016666247e-5
 0.001099395221772342
 0.011038622307372232
 0.07387132069817183
 0.20198030129465658
 0.35954472556518147
 0.5167641452122544
 0.6637713825376198
 0.7841482271226826
 0.8737784044472122
 0.9331779821204208
 0.9684489616227439
 0.9869310637824636
 0.9932596937939658
```

---

```julia
function lorenz(du, u, p, t)
    du[1] = 10.0 * (u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)

sol = solve(prob)
```

```
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 1263-element Vector{Float64}:
   0.0
   3.5678604836301404e-5
   0.0003924646531993154
   0.0032624077544510573
   0.009058075635317072
   0.01695646895607931
   0.02768995855685593
   0.04185635042021763
   0.06024041165841079
   0.08368541255159562
   ⋮
  99.30760258626904
  99.39665422328268
  99.49536147459878
  99.58822928767293
  99.68983993598462
  99.77864535713971
  99.85744078539504
  99.93773320913628
 100.0
u: 1263-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9996434557625105, 0.0009988049817849058, 1.781434788799208e-8]
 [0.9961045497425811, 0.010965399721242457, 2.146955365838907e-6]
 [0.9693591634199452, 0.08977060667778931, 0.0001438018342266937]
 [0.9242043615038835, 0.24228912482984957, 0.0010461623302512404]
 [0.8800455868998046, 0.43873645009348244, 0.0034242593451028745]
 [0.8483309847495312, 0.6915629321083602, 0.008487624590227805]
 [0.8495036669651213, 1.0145426355349096, 0.01821208962127994]
 [0.9139069574560097, 1.4425599806525806, 0.03669382197085303]
 [1.088863826836895, 2.052326595543049, 0.0740257368585531]
 ⋮
 [4.669609096878053, 3.061564434452441, 25.1424735017959]
 [4.188801916573263, 4.617474401440693, 21.09864175382292]
 [5.559603854699961, 7.905631612648314, 18.79323210016923]
 [8.556629716266505, 12.533041060088328, 20.6623639692711]
 [12.280585075547771, 14.505154761545633, 29.332088452699942]
 [11.736883151600804, 8.279294641640229, 34.68007510231878]
 [8.10973327066804, 3.2495066495235854, 31.97052076740117]
 [4.958629886040755, 2.194919965065022, 26.948439650907677]
 [3.8020065515435855, 2.787021797920187, 23.420567509786622]
```

---

```julia
plot(sol, vars = (1, 2, 3))
```

![lorenz](lorenz.png)

---

```julia
using ParameterizedFunctions

lotka_volterra = @ode_def begin
  d🐁 = α*🐁  - β*🐁*🐈
  d🐈 = -γ*🐈 + δ*🐁*🐈
end α β γ δ

u0 = [1.0, 1.0] # Initial condition

tspan = (0.0, 10.0) # Simulation interval
tsteps = 0.0:0.1:10.0 # intermediary points

p = [1.5, 1.0, 3.0, 1.0] # equation parameters: p = [α, β, δ, γ]

prob = ODEProblem(lotka_volterra, u0, tspan, p)
sol = solve(prob)
```

```
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 34-element Vector{Float64}:
  0.0
  0.0776084743154256
  0.23264513699277584
  0.4291185174543143
  0.6790821987497083
  0.9444046158046306
  1.2674601546021105
  1.6192913303893046
  1.9869754428624007
  2.2640902393538296
  ⋮
  7.584863345264154
  7.978068981329682
  8.48316543760351
  8.719248247740158
  8.949206788834692
  9.200185054623292
  9.438029017301554
  9.711808134779586
 10.0
u: 34-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0454942346944578, 0.8576684823217127]
 [1.1758715885138267, 0.639459570317544]
 [1.4196809607170826, 0.4569962601282084]
 [1.876719395008001, 0.32473342927911314]
 [2.5882500645533466, 0.26336255535952163]
 [3.8607089092207665, 0.2794458098285253]
 [5.750812667710396, 0.5220072537934558]
 [6.814978999130169, 1.9177826328390666]
 [4.3929992925714245, 4.194670792850584]
 ⋮
 [2.614253967788294, 0.26416945387525886]
 [4.241076127191749, 0.3051236762921916]
 [6.791123785297795, 1.1345287797146113]
 [6.265370675764892, 2.74169350754023]
 [3.7807651118880545, 4.431165685863461]
 [1.816420140681761, 4.064056625315978]
 [1.1465021407690728, 2.7911706616216976]
 [0.9557986135403302, 1.6235622951850799]
 [1.0337581256020607, 0.9063703842886133]
```

---

# Type-Dispatch Programming

* Centered around implementing the generic template of the algorithm not around building representations of data.
  * The data type choose how to efficiently implement the algorithm.
  * With this feature share and reuse code is very easy

[JuliaCon 2019 | The Unreasonable Effectiveness of Multiple Dispatch | Stefan Karpinski](https://youtu.be/kc9HwsxE1OY)

---

Simple gravity pendulum

```julia
using DifferentialEquations, Plots

g = 9.79 # Gravitational constants
L = 1.00 # Length of the pendulum

#Initial Conditions
u₀ = [0, π / 60] # Initial speed and initial angle
tspan = (0.0, 6.3) # time domain

#Define the problem
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6)
```

```
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 51-element Vector{Float64}:
 0.0
 0.031087055826046293
 0.08561410455758335
 0.15017614691650324
 0.22845000693621464
 0.31505859610319265
 0.40952552265148245
 0.5094804665223178
 0.6142606828816117
 0.7233814675983777
 ⋮
 5.276991527121083
 5.427396011335077
 5.568207973942741
 5.706096737512237
 5.841862506521234
 5.97857442065412
 6.118893944150791
 6.267408889178013
 6.3
u: 51-element Vector{Vector{Float64}}:
 [0.0, 0.05235987755982988]
 [0.0016251489879623073, 0.052112381623717514]
 [0.004429323578461451, 0.05049245161124966]
 [0.007577026053634464, 0.04668511432388565]
 [0.010968722238292433, 0.03954356773072573]
 [0.013951484712956963, 0.02891353707715841]
 [0.01603824120236177, 0.014943963379430084]
 [0.01672973991691839, -0.0012207134336849744]
 [0.015713051180260403, -0.018011260476534265]
 [0.012878614479109221, -0.03343350596936019]
 ⋮
 [-0.012041966804622781, -0.03635808653571255]
 [-0.016001764791371276, -0.015321724384266154]
 [-0.016561978903112534, 0.007493947984585336]
 [-0.014042737908113273, 0.028477612500254594]
 [-0.009043594380336785, 0.04405515312673323]
 [-0.0023878772866644496, 0.051823979712627824]
 [0.0048791392846661815, 0.050084788515492334]
 [0.011535290105166612, 0.037932317215479905]
 [0.012709480659187212, 0.034061109155362805]
```

---

Analytic and computed solution

```julia
u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

scatter(sol.t, getindex.(sol.u, 2), label = "Numerical")
plot!(sol.t, u, label = "Analytic")
```

![](pendulum1.svg)

---

[Numbers with Uncertainties](http://tutorials.juliadiffeq.org/html/type_handling/02-uncertainties.html)

```julia
using Measurements

g = 9.79 ± 0.02; # Gravitational constants
L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions
u₀ = [0 ± 0, π / 60 ± 0.01] # Initial speed and initial angle

#Define the problem
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6);
```

---

Analytic solution

```julia
u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")
plot!(sol.t, u, label = "Analytic")
```

![](pendulum2.svg)

---

# Poisson Equation

$$ \frac{\partial^2 u}{\partial x^2} = b  \qquad x \in [0,1] $$

$$ u(0) = u(1) = 0, \qquad b = \sin(2\pi x) $$

```julia
using Plots, SparseArrays

Δx = 0.05
x = Δx:Δx:1-Δx ## Solve only interior points: the endpoints are set to zero.
n = length(x)
B = sin.(2π*x) * Δx^2

P = spdiagm( -1 =>    ones(Float64,n-1),
              0 => -2*ones(Float64,n),
              1 =>    ones(Float64,n-1))

u1 = P \ B
```

```
19-element Vector{Float64}:
 -0.007892189393343805
 -0.015011836300750241
 -0.020662020077425496
 -0.024289661368163382
 -0.025539661368163383
 -0.024289661368163386
 -0.0206620200774255
 -0.015011836300750248
 -0.007892189393343811
 -6.308085367188389e-18
  0.0078921893933438
  0.01501183630075024
  0.020662020077425496
  0.024289661368163382
  0.025539661368163387
  0.024289661368163386
  0.020662020077425496
  0.015011836300750243
  0.007892189393343806
```

---

```julia
plot([0;x;1],[0;u1;0], label="computed")
scatter!([0;x;1],-sin.(2π*[0;x;1])/(4π^2),label="exact")
```

![](poisson1.png)

---

# DiffEqOperators.jl

```julia
using DiffEqOperators

Δx = 0.05
x = Δx:Δx:1-Δx ## Solve only interior points: the endpoints are set to zero.
n = length(x)
b = sin.(2π*x)

# Second order approximation to the second derivative
order = 2
deriv = 2

Δ = CenteredDifference{Float64}(deriv, order, Δx, n)
bc = Dirichlet0BC(Float64)

u2 = (Δ * bc) \ b
```

---

```julia
plot([0;x;1],[0;u2;0], label="computed")
scatter!([0;x;1],-sin.(2π*[0;x;1])/(4π^2),label="exact")
```

![](poisson2.png)

---

# HOODESolver.jl

The objective of this Julia package is to valorize the recent developments carried out within [INRIA team MINGuS](https://team.inria.fr/mingus/) on Uniformly Accurate numerical methods (UA) for highly oscillating problems. We propose to solve the following equation

$$
\frac{d u(t)}{dt} = \frac{1}{\varepsilon} A u(t) + f(t, u(t)), \qquad u(t=t_0)=u_0, \qquad \varepsilon\in ]0, 1], \qquad (1)
$$

with

* $u : t \in [t_0, t_1] \mapsto u(t)\in \mathbb{R}^n, \quad t_0, t_1 \in \mathbb{R}$,
  * $u_0 \in \mathbb{R}^n$,
  * $A\in {\mathcal{M}}_{n,n}(\mathbb{R})$ is such that $\tau \mapsto \exp(\tau A)$ is $2 \pi$-periodic,
  * $f : (t, u) \in  \mathbb{R}\times \mathbb{R}^n \mapsto \mathbb{R}^n$.

https://pnavaro.github.io/HOODESolver.jl/stable/

Philippe Chartier, Nicolas Crouseilles, Mohammed Lemou, Florian Mehats and Xiaofei Zhao.

Package: Yves Mocquard and Pierre Navaro.

---

## Two-scale formulation

First, rewrite equation (1) using the variable change $w(t)=\exp(-(t-t_0)A/\varepsilon) u(t)$ to obtain

$$
\frac{d w(t)}{dt} = F\Big(\frac{t-t_0}{\varepsilon}, w(t) \Big), $$

$$w(t_0) = u_0, \varepsilon \in ]0, 1], $$

where the function $F$ is expressed from the data of the original problem (1)

$$F\Big( \frac{s}{\varepsilon}, w \Big) = \exp(-sA/\varepsilon) \; f( \exp(sA/\varepsilon), \; w).
$$

We then introduce the function $U(t, \tau), \tau\in [0, 2 \pi]$ such that $U(t, \tau=(t-t_0)/\varepsilon) = w(t)$. The two-scale function is then the solution of the following equation.

$$
\frac{\partial U}{\partial t} + \frac{1}{\varepsilon} \frac{\partial U}{\partial \tau} =  F( \tau, U), \;\;\; U(t=t_0, \tau)=\Phi(\tau), \;\; \varepsilon\in ]0, 1], \;\;\;\;\;\;\;\;\;\; (2)
$$

where $\Phi$ is a function checking $\Phi(\tau=0)=u_{0}$ chosen so that the $U$ solution of (2) is smooth.

---

# Hénon-Heiles Example

We consider the system of Hénon-Heiles satisfied by $u(t)=(u_1, u_2, u_3, u_4)(t)$.

$$
\frac{d u }{dt} = \frac{1}{\varepsilon} Au + f(u),  $$
$$ u(t_0)=u_0 \in \mathbb{R}^4,
$$

where $A$ and $f$ are selected as follows

$$
A=
\\begin{pmatrix}
0 & 0 & 1 & 0  \\\\
0 & 0 & 0 & 0  \\\\
-1 & 0 & 0 & 0  \\\\
0 & 0 & 0 & 0
\\end{pmatrix}, \qquad
f(u) = \left(
\begin{array}{cccc}
0 \\\\
u_4\\\\
-2 u_1 u_2\\\\
-u_2-u_1^2+u_2^2
\end{array}
\right).
$$

---

# SplitODEProblem

The `SplitODEProblem` type from package [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/types/split_ode_types/) offers a interface for this kind of problem.

```julia
using Plots, DifferentialEquations

epsilon = 0.002
A = [ 0 0 1 0 ;
      0 0 0 0 ;
     -1 0 0 0 ;
      0 0 0 0 ]

f1 = DiffEqArrayOperator( A ./ epsilon)

function f2(du, u, p, t)
    du[1] = 0
    du[2] = u[4]
    du[3] = 2*u[1]*u[2]
    du[4] = -u[2] - u[1]^2 + u[2]^2
end

tspan = (0.0, 0.1)

u0 = [0.55, 0.12, 0.03, 0.89]

prob1 = SplitODEProblem(f1, f2, u0, tspan);
sol1 = solve(prob1, ETDRK4(), dt=0.001);
```

---

```julia
using HOODESolver, Plots

A = [ 0 0 1 0 ;
      0 0 0 0 ;
     -1 0 0 0 ;
      0 0 0 0 ]

f1 = LinearHOODEOperator( epsilon, A)

prob2 = SplitODEProblem(f1, f2, u0, tspan);

sol2 = solve(prob2, HOODEAB(), dt=0.01);
```

```
solve function prob=HOODESolver.HOODEProblem{Float64}(HOODESolver.HOODEFunction{true, 4}(Main.f2), [0.55, 0.12, 0.03, 0.89], (0.0, 0.1), missing, [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0], 0.002, missing),
 nb_tau=32, order=4, order_prep=6, dense=true,
 nb_t=10, getprecision=true, verbose=100

10/10

12/12
```

---

```julia
plot(sol1, vars=[3], label="EDTRK4")
plot!(sol2, vars=[3], label="HOODEAB")
plot!(sol2.t, getindex.(sol2.u, 3), m=:o, label="points")
```

```
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.
│   caller = ip:0x0
└ @ Core :-1
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.
│   caller = ip:0x0
└ @ Core :-1
```

![](sol1sol2.png)

---

# High precision

```julia
u0 = BigFloat.([90, -44, 83, 13]//100)
t_end = big"1.0"
epsilon = big"0.0017"

prob = HOODEProblem(f2, u0, (big"0.0",t_end), missing, A, epsilon)
```

```
HOODEProblem with uType Vector{BigFloat} and tType BigFloat. In-place: nothing
timespan: (0.0, 1.0)
u0: 4-element Vector{BigFloat}:
  0.8999999999999999999999999999999999999999999999999999999999999999999999999999965
 -0.4400000000000000000000000000000000000000000000000000000000000000000000000000014
  0.830000000000000000000000000000000000000000000000000000000000000000000000000001
  0.1300000000000000000000000000000000000000000000000000000000000000000000000000006
```

The float are coded on 512 bits.

---

## Precision of the result with ε = 0.015

![](assets/error_order.png)

JOSS paper https://joss.theoj.org/papers/10.21105/joss.03077

Much of HOODESolver.jl was implemented by Y. Mocquard while he was supported by Inria through the AdT (Aide au développement technologique) J-Plaff of the center Rennes- Bretagne Atlantique.

---

# Why use Julia language!

* **You develop in the same language in which you optimize.**
  * Packaging system is very efficient (5858 registered packages)
  * PyPi (311,500 projects) R (17739 packages)
  * It is very easy to create a package (easier than R and Python)
  * It is very easy to access to a GPU device.
  * Nice interface for Linear Algebra and Differential Equations
  * Easy access to BLAS, LAPACK and scientific computing librairies.
  * Julia talks to all major Languages - mostly without overhead!

---

# What's bad

* It is still hard to build shared library or executable from Julia code.
  * Compilation times lattency. Using [Revise.jl](https://github.com/timholy/Revise.jl) helps a lot.
  * Plotting takes time (5 seconds for the first plot)
  * OpenMP is better than the Julia multithreading library but it is progressing.
  * Does not work well with vectorized code, you need to do a lot of inplace computation to avoid memory allocations and use explicit views to avoid copy. There are some packages like [LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl).

[What's Bad About Julia by Jeff Bezanson](https://www.youtube.com/watch?v=TPuJsgyu87U)

---

## From zero to Julia!

https://techytok.com/from-zero-to-julia/

## Python-Julia benchmarks by Thierry Dumont

https://github.com/Thierry-Dumont/BenchmarksPythonJuliaAndCo/wiki

## Mailing List

* https://listes.services.cnrs.fr/wws/info/julia

---

*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*