Langage Julia

Cours ENPC - Pratique du calcul scientifique

Organisation du cours (1/2)

Organisation de chaque séance

Chaque séance comprendra

  • une partie théorique avec diapositives,

  • une partie pratique sur machine.

Modalités d’évaluation

  • 3 devoirs (3 x 10%) sous forme de Jupyter notebooks:

    • Devoir 1 remis aujourd’hui,

    • Devoir 2 remis en 3ème séance (29 avril),

    • Devoir 3 remis en 5ème séance (13 mai).

  • Examen final sur machine, sans IA (70%)

Organisation du cours (2/2)

Programme détaillé

Date Sujet
Mardi 8 avril Introduction à Julia (devoir 1 \(\rightarrow\) lundi 14 avril)
Mardi 15 avril Interpolation et intégration
Mardi 29 avril Systèmes linéaires (devoir 2 \(\rightarrow\) lundi 12 mai)
Mardi 6 mai Problèmes aux valeurs propres
Mardi 13 mai Systèmes non-linéaires (devoir 3 \(\rightarrow\) lundi 26 mai)
Mardi 20 mai Équations différentielles ordinaires
Mardi 27 mai Examen final

Notes importantes

  • Le contenu du premier notebook est essentiel pour la suite du cours.

  • Le cours avance rapidement et demande un travail régulier.

  • La présence en séance est obligatoire (-2/20 par absence non justifiée).

  • Langage de programmation créé en 2009 au MIT

  • Par Jeff Bezanson, Stefan Karpinski, Viral B. Shah, Alan Edelman

  • Version 1.11.4 depuis le 10/03/2025

  • Documentation générale de Julia : https://docs.julialang.org/en/v1/

Installation de Julia → deux options

  1. Gestionnaire de versions juliaup
  1. Versions individuelles depuis https://julialang.org/

Editeur de développement (IDE)

Le REPL et quelques commandes importantes à connaître

  • REPL = read–eval–print loop → console de Julia

  • Le REPL est accessible

    • de manière autonome en tant qu’exécutable ou en tapant julia dans une console (si PATH à jour)

    • intégré dans VSCode en tapant Maj+Ctrl+P puis en cherchant Julia: Start REPL

      Dans VSCode, on peut taper ses lignes de code dans un fichier *.jl et les lancer une par une dans la console intégrée par Maj+Entrée.

  • Depuis le REPL, taper

    • ] donne accès au gestionnaire de bibliothèques

    • ; donne accès au shell mode

    • ? donne accès à l’aide en ligne

    • ⟵ retourne au mode normal

    • Tab pour la complétion automatique

    • les flèches haut et bas pour naviguer dans l’historique (éventuellement filtré par le début de ligne tapé)

Le gestionnaire de bibliothèques

Pour appeler une bibliothèque déjà installée

  • avec import
import LinearAlgebra
LinearAlgebra.inv([4 5 ; 3 4])
2×2 Matrix{Float64}:
  4.0  -5.0
 -3.0   4.0
  • ou avec using (plus besoin de préfixer le nom de la bibliothèque ensuite)
using LinearAlgebra
inv([4 5 ; 3 4])
2×2 Matrix{Float64}:
  4.0  -5.0
 -3.0   4.0

Pour installer une nouvelle bibliothèque : 2 options

  1. dans le gestionnaire accessible depuis le REPL : taper ] puis
pkg> add Plots
  1. ou directement dans le code
import Pkg
Pkg.add("Plots")
add Plots

Important

L’installation d’une bibliothèque dans un environnement ne se fait qu’une seule fois.

Toute tentative de réinstallation est non seulement inutile mais surtout chronophage.

Pourquoi Julia ?

Python Matlab Julia C
Libre ? ✔️ ❌ ✔️ ✔️
Rapide ? ❔ ❔ ✔️ ✔️
Facile à apprendre ? ✔️ ✔️ ✔️ ❌
Adapté au calcul scientifique ? ❔ ✔️ ✔️ ❌
Richesse de l’écosystème ? ✔️ ❔ ❔ ✔️

  • Julia est facile à apprendre, à lire et à maintenir (comparable à Python)

  • Julia est rapide (comparable à C/C++, Fortran)

  • Julia est efficace (dispatch multiple, macros,…)

  • Julia est bien adapté au calcul scientifique

Julia est facile à apprendre, à lire et à maintenir

  • Typage dynamique
x = 1 ; y = π ; z = 1//2 ; t = x + y + z
4.641592653589793
for v ∈ (x,y,z,t) println("$v is a $(typeof(v))") end
1 is a Int64
π is a Irrational{:π}
1//2 is a Rational{Int64}
4.641592653589793 is a Float64
  • Mais possibilité de typage statique
f(x) = 2x
f(x::Float64) = 3x
@show f(1), f(1.) ;
(f(1), f(1.0)) = (2, 3.0)
  • Définition simplifiée de fonctions courtes (équivalent du lambda de python) ou de fonctions anonymes
print( (x -> 2x^2)(5) );
50

  • Bibliothèque de base contient les tableaux et leurs opérations courantes (+, -, *, \) puis LinearAlgebra pour det, eigen…

  • Compréhension de tableau (for dans le tableau)

A = [63/(i+2j) for i ∈ 1:2, j ∈ 1:2]
2×2 Matrix{Float64}:
 21.0   12.6
 15.75  10.5
x = [1.2, 3.4]
b = A*x
A\b
2-element Vector{Float64}:
 1.2000000000000008
 3.399999999999998
@show A\b == x
@show A\b ≈ x ;
A \ b == x = false
A \ b ≈ x = true
  • Numérotation commence à 1 et remplissage par colonne
@show x[1]
@show x[2] ;
x[1] = 1.2
x[2] = 3.4
  • Caractères unicode
𝒢(x,μ=0.,σ=1.) = 1/σ√(2π) * exp(-(x-μ)^2/2σ^2)
𝒢 (generic function with 3 methods)
@show Σaᵢⱼ² = sum(A.^2)
@show √Σaᵢⱼ²
ω = 1  # ⚠️ même les emojis !
f(t) = sin(ω*t)
∇f(t) = ω * cos(ω*t)
f²(t) = f(t) * f(t) ;
Σaᵢⱼ² = sum(A .^ 2) = 958.0725
√Σaᵢⱼ² = 30.952746243265718
  • Nouveaux opérateurs binaires
⊕(α,β) = α + β
5 ⊕ 7
12
⊕(α::Array{T,2}, β::Array{T,2}) where {T} = α + β'
⊕ (generic function with 2 methods)
display(A + A)
display(A ⊕ A)
2×2 Matrix{Float64}:
 42.0  25.2
 31.5  21.0
2×2 Matrix{Float64}:
 42.0   28.35
 28.35  21.0

Julia est rapide

CJuliaLuaJITRustGoFortranJavaJavaScriptMatlabMathematicaPythonROctaveiteration_pi_summatrix_multiplymatrix_statisticsparse_integersprint_to_filerecursion_fibonaccirecursion_quicksortuserfunc_mandelbrotbenchmark100101102103104

⚠ Ces résultats ne tiennent pas compte du temps de compilation.

The vertical axis shows each benchmark time normalized against the C implementation. The benchmark data shown above were computed with Julia v1.0.0, SciLua v1.0.0-b12, Rust 1.27.0, Go 1.9, Java 1.8.0_17, Javascript V8 6.2.414.54, Matlab R2018a, Anaconda Python 3.6.3, R 3.5.0, and Octave 4.2.2. C and Fortran are compiled with gcc 7.3.1, taking the best timing from all optimization levels (-O0 through -O3). C, Fortran, Go, Julia, Lua, Python, and Octave use OpenBLAS v0.2.20 for matrix operations; Mathematica uses Intel MKL. The Python implementations of matrix_statistics and matrix_multiply use NumPy v1.14.0 and OpenBLAS v0.2.20 functions; the rest are pure Python implementations. Raw benchmark numbers in CSV format are available under https://github.com/JuliaLang/Microbenchmarks.

These micro-benchmark results were obtained on a single core (serial execution) on an Intel Core i7-3960X 3.30GHz CPU with 64GB of 1600MHz DDR3 RAM, running openSUSE LEAP 15.0 Linux.

Julia est efficace

Un atout majeur est le multiple dispatch

Exemple tiré de la conférence The Unreasonable Effectiveness of Multiple Dispatch de Stefan Karpinski

abstract type Animal end
struct Chien <: Animal; nom::String end
struct Chat <: Animal; nom::String end

function rencontre(a::Animal,b::Animal)
    verb = agit(a,b)
    println("$(a.nom) rencontre $(b.nom) et $verb")
end

agit(a::Chien,b::Chien) = "le renifle" ; agit(a::Chien,b::Chat) = "le chasse"
agit(a::Chat,b::Chien) = "s'enfuit" ; agit(a::Chat,b::Chat) = "miaule"

medor = Chien("Médor") ; 🐶 = Chien("🐶")
felix = Chat("Félix") ; 🐱 = Chat("🐱")

rencontre(medor,🐶)
rencontre(🐶,felix)
rencontre(felix,🐶)
rencontre(🐱,felix)
Médor rencontre 🐶 et le renifle
🐶 rencontre Félix et le chasse
Félix rencontre 🐶 et s'enfuit
🐱 rencontre Félix et miaule

Exemple adapté du blog de Mosè Giordano

abstract type HandShape end
struct Rock     <: HandShape end
struct Paper    <: HandShape end
struct Scissors <: HandShape end
play(::Type{Paper}, ::Type{Rock}) = println("Paper VS Rock ⇒ Paper wins")
play(::Type{Scissors}, ::Type{Paper}) = println("Scissors VS Paper ⇒ Scissors wins")
play(::Type{Rock}, ::Type{Scissors}) = println("Rock VS Scissors ⇒ Rock wins")
play(::Type{T}, ::Type{T}) where {T<: HandShape} = println("Tie between $(T), try again")
play(a::Type{<:HandShape}, b::Type{<:HandShape}) = play(b, a) # Commutativity
play(Paper, Rock)
play(Scissors, Scissors)
play(Paper, Scissors)
Paper VS Rock ⇒ Paper wins
Tie between Scissors, try again
Scissors VS Paper ⇒ Scissors wins

Il est aisé d’ajouter de nouveaux types a posteriori

struct Well <: HandShape end
play(::Type{Well}, ::Type{Rock})     = "Well VS Rock ⇒ Well wins"
play(::Type{Well}, ::Type{Scissors}) = "Well VS Scissors ⇒ Well wins"
play(::Type{Paper}, ::Type{Well})    = "Paper VS Well ⇒ Paper wins"
play(Scissors, Well)
"Well VS Scissors ⇒ Well wins"

Exemple plus mathématique de multiple dispatch

Exemple inspiré de la conférence The Unreasonable Effectiveness of Multiple Dispatch par Stefan Karpinski

\[\sum_{i=1}^p\mathbf{\boldsymbol{v}}_i^\mathrm{T} \mathsf{A} \mathbf{\boldsymbol{v}}_i, \quad \mathsf{A}∈{\mathbb{{R}}}^{n×n},\;\mathbf{\boldsymbol{v}}_i∈{\mathbb{{R}}}^n\]

  • Calcul standard
using LinearAlgebra

inner(v, A, w) = v'A*w

# function inner_sum(A, vs)
#     t = zero(eltype(A))
#     for v ∈ vs t += inner(v, A, v) end
#     return t
# end

inner_sum(A, vs) = sum(inner(v, A, v) for v ∈ vs)

n = 10_000 ; p = 100
A = rand(n, n)
vs = [rand(n) for _ ∈ 1:p]

@time inner_sum(A, vs)
  7.854520 seconds (57.46 k allocations: 10.482 MiB, 0.41% compilation time)
1.251779779546054e9

Cas particulier de vecteurs \(\mathbf{\boldsymbol{v}}_i=(0,\cdots,0,1,0,\cdots,0)^\mathrm{T}\)

\[\mathbf{\boldsymbol{v}}^\mathrm{T} \mathsf{A} \mathbf{\boldsymbol{w}} = A_{\mathrm{ind}(\mathbf{\boldsymbol{v}}), \mathrm{ind}(\mathbf{\boldsymbol{w}})}\]

  • Construction d’un nouveau type de vecteur (size et getindex doivent être définis pour les calculs matriciels usuels)
import Base: size, getindex
struct OneVector <: AbstractVector{Bool}
  len::Int
  ind::Int
end
size(v::OneVector) = (v.len,)
getindex(v::OneVector, i::Integer) = i == v.ind
print(OneVector(5, rand(1:5)))
Bool[0, 0, 1, 0, 0]
  • Calcul standard
vs = [OneVector(n, rand(1:n)) for _ ∈ 1:p]
@time inner_sum(A, vs)
  0.026785 seconds (96.88 k allocations: 5.091 MiB, 99.90% compilation time)
52.43257200179093
  • Dispatch multiple : spécialisation de inner
inner(v::OneVector, A, w::OneVector) = A[v.ind, w.ind]
@time inner_sum(A, vs)
@time inner_sum(A, vs)
  0.019712 seconds (84.77 k allocations: 4.537 MiB, 99.68% compilation time: 100% of which was recompilation)
  0.000011 seconds (1 allocation: 16 bytes)
52.43257200179093

Julia est bien adapté au calcul scientifique

  • Différentiation automatique
using Zygote
f(x) = log(x)
f'(2), f''(2)
(0.5, -0.25)
  • Calculs dimensionnels
using Unitful
m = u"m" ; cm = u"cm" ;
@show 1.0m + 1.0cm
println()
@show 1.0u"MPa" + 2.0u"bar" + 3.0u"daN/cm^2"
println()
@show uconvert(u"MPa",1u"bar")
;
1.0m + 1.0cm = 1.01 m

1.0 * u"MPa" + 2.0 * u"bar" + 3.0 * u"daN/cm^2" = 1.5e6 kg m^-1 s^-2

uconvert(u"MPa", 1 * u"bar") = 1//10 MPa
  • Equations différentielles
using ModelingToolkit, DifferentialEquations, Plots
@parameters t
∂_∂(x) = y -> expand_derivatives(Differential(x)(y))
d_dt = ∂_∂(t)
q = @variables θ(t) ϕ(t) ψ(t)
θ̇, ϕ̇, ψ̇ = q̇ = d_dt.(q)
q̈ = d_dt.(q̇)
@parameters g R M m₁ m₂ ℓ₁ ℓ₂ L ;

… éq. d’Euler-Lagrange \(\frac{{\mathrm{d}}}{{\mathrm{d}}t}\left(\frac{∂ℒ}{∂\dot{q_i}}\right) -\frac{∂ℒ}{∂q_i}=0\)

K = Jᴵ * θ̇^2 / 2 + m₁ * (ẋ₁^2 + ẏ₁^2) / 2 + m₂ * (ẋ₂^2 + ẏ₂^2) / 2
V = -M * g * yᴳ * cos(θ) + m₁ * g * y₁ + m₂ * g * y₂
ℒ = K - V
EL = [d_dt(∂_∂(v̇)(ℒ)) - ∂_∂(v)(ℒ) for (v̇, v) ∈ zip(q̇, q)] ;

…

  • Calcul tensoriel
using SymPy, LinearAlgebra, TensND
Spherical = coorsys_spherical()
θ, ϕ, r = getcoords(Spherical)
𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical)
@set_coorsys Spherical
𝕀, 𝕁, 𝕂 = ISO() ; 𝟏 = tensId2()
k, μ = symbols("k μ", positive = true)
λ = k -2μ/3
u = SymFunction("u", real = true)
𝛆 = SYMGRAD(u(r) * 𝐞ʳ)
𝛆 |> intrinsic
(u(r)/r)𝐞ᶿ⊗𝐞ᶿ + (u(r)/r)𝐞ᵠ⊗𝐞ᵠ + (Derivative(u(r), r))𝐞ʳ⊗𝐞ʳ
𝛔 = λ * tr(𝛆) * 𝟏 + 2μ * 𝛆
eq = DIV(𝛔) ⋅ 𝐞ʳ
sol = dsolve(eq, u(r))

\(u{\left(r \right)} = \frac{C_{1}}{r^{2}} + C_{2} r\)

T̂ = tfactor(tsimplify(tsubs(𝐞ʳ ⋅ 𝛔 ⋅ 𝐞ʳ, u(r) => sol.rhs())))

\(\frac{- 4 C_{1} μ + 3 C_{2} k r^{3}}{r^{3}}\)

Références pour auto-formation