Semaine 1 : Introduction à Julia¶

In [1]:
using Plots
using LaTeXStrings

Comprendre la sortie des cellules Jupyter¶

  • La sortie d'une série de commandes Julia est la dernière expression :
In [2]:
# Jupyter will display the result of `a + b` after running the cell

a = 5
b = 4
c = a + b
Out[2]:
9
  • L'utilisation de ; en fin de commande supprime l'affichage
In [3]:
# Jupyter will not display anything

a = 5
b = 4
a + b;
  • Cela s'applique aussi aux graphiques. Comparer les deux situations suivantes :
In [4]:
plot(cos)
Out[4]:
No description has been provided for this image
In [5]:
plot(sin);
  • On peut aussi utiliser @show ou print :
In [6]:
a = 5
b = 4
@show a
println(b)
a + b

a = 5
4
Out[6]:
9

Manipulation de chaînes de caractères¶

  • Les chaînes de caractères sont définies en Julia entre guillemets doubles ". Contrairement à Python, les guillemets simples ' ne sont utilisés que pour les caractères simples.
In [7]:
s = "Hello world!"
Out[7]:
"Hello world!"
In [8]:
'o' ∈ s # ∈ is an alias for in
Out[8]:
true
  • Contrairement à de nombreux langages qui utilisent +, Julia utilise * pour la concaténation de chaînes. L'idée est de souligner la non-commutativité de cette opération
In [9]:
s1 = "Hello " ; s2 = "world!"
s = s1 * s2
println(s)

Hello world!
In [10]:
print(s1^5)

Hello Hello Hello Hello Hello
  • L'interpolation dans une chaîne s'obtient avec $
In [11]:
col = "red" ; x = 20_000 ; discount = 5
s = "My $col car cost $(x*(1-discount/100)) €"
print(s)

My red car cost 19000.0 €
In [12]:
s = "Hello world!"
ordinal_suffix(i) = i == 1 ? "st" : i == 2 ? "nd" : "th"
for (i, c) ∈ enumerate(s) 
    println("The character at $i$(ordinal_suffix(i)) place is $c") 
end

The character at 1st place is H
The character at 2nd place is e
The character at 3th place is l
The character at 4th place is l
The character at 5th place is o
The character at 6th place is  
The character at 7th place is w
The character at 8th place is o
The character at 9th place is r
The character at 10th place is l
The character at 11th place is d
The character at 12th place is !

Notez que les chaînes LaTeX peuvent être construites grâce au package LaTeXStrings, importé ci-dessus et utilisé plus loin dans ce notebook.

Utilisation et définition de fonctions¶

[Exercice 1]¶

Écrire une fonction qui retourne la somme des n premiers entiers au carré en utilisant une boucle for.

Indices :

  • lire la documentation de function (ouvrir un REPL et taper ? suivi du mot recherché, ici function),
  • un intervalle d'entiers de m à n s'écrit simplement m:n (contrairement à Python, les bornes sont incluses dans l'intervalle).
In [13]:
function sum_of_squares(n)
    ### BEGIN SOLUTION
    result = 0
    for i in 1:n
        result += i*i
    end
    return result
    ### END SOLUTION
end;
In [14]:
@assert sum_of_squares(1) == 1 "❌"
@assert sum_of_squares(2) == 5 "❌"
@assert sum_of_squares(5) == 55 "❌"
@assert sum_of_squares(10) == 385 "❌"

Remarques :

  • La même fonction peut s'écrire sous forme compacte
In [15]:
short_sum_of_squares(n) = (1:n)'*(1:n)

@show short_sum_of_squares(2)
@show short_sum_of_squares(4);

short_sum_of_squares(2) = 5
short_sum_of_squares(4) = 30
  • Lire la documentation de la fonction sum pour comprendre pourquoi le code suivant fonctionne aussi.
In [16]:
other_sum_of_squares(n) = sum(x -> x^2, 1:n)

@show other_sum_of_squares(2)
@show other_sum_of_squares(4);

other_sum_of_squares(2) = 5
other_sum_of_squares(4) = 30
  • La fonction peut être appliquée élément par élément à un tableau en utilisant l'opérateur de diffusion .
In [17]:
@show short_sum_of_squares.([2, 4, 5])

# The above command is syntactic sugar for
@show broadcast(short_sum_of_squares, [2, 4, 5]);

short_sum_of_squares.([2, 4, 5]) = [5, 30, 55]
broadcast(short_sum_of_squares, [2, 4, 5]) = [5, 30, 55]
  • La diffusion fonctionne même pour des fonctions à plusieurs arguments
In [18]:
geometric_mean(a, b) = √(a*b)
@show geometric_mean.(2, [8, 18]);

geometric_mean.(2, [8, 18]) = [4.0, 6.0]
  • Il est parfois commode de définir des fonctions anonymes avec la notation -> :
In [19]:
# Plot the sum of cubes for different n
scatter(n -> sum(x -> x^3, 1:n), 1:5, label="sum of cubes")
Out[19]:
No description has been provided for this image

[Exercice 2]¶

La méthode babylonienne pour calculer la racine carrée de $a > 0$ est donnée par l'itération $$ \forall k\in\mathbb{N},\qquad x_{k+1}=\frac{1}{2}\left(x_k+\frac{a}{x_k}\right) $$

Écrire une fonction my_sqrt(a, x₀) qui retourne une approximation de la racine carrée de a par 10 itérations de cette méthode, en partant de x₀.

In [20]:
function my_sqrt(a, x₀)
    ### BEGIN SOLUTION
    x = x₀
    for i in 1:10
        x = (x + a/x)/2
    end
    return x
    ### END SOLUTION
end;
In [21]:
@assert my_sqrt(4, 1) ≈ 2 "❌"
@assert my_sqrt(144, 1) ≈ 12 "❌"

Tracé de courbes et création d'animations¶

  • La fonction plot peut être appelée de différentes façons. Une première approche consiste à lui passer une fonction :
In [22]:
f(x) = cos(x) * exp(-abs(x))
plot(f, xlims=(-5, 5))
Out[22]:
No description has been provided for this image
  • Une autre option, souvent plus flexible, est de lui passer des tableaux $x$ et $y$ :
In [23]:
f(x) = cos(x) * exp(-abs(x))
x = -5:.01:5
plot(x, f.(x))
Out[23]:
No description has been provided for this image
  • Utiliser les commandes suivies de ! pour modifier un graphique existant.
In [24]:
plot(cos, label = "f1")
plot!(sin, label = "f2")
xlims!(-5, 5)
title!("title")
xlabel!("x")
ylabel!("y")
Out[24]:
No description has been provided for this image
  • Les options peuvent être regroupées dans un appel à plot! :
In [25]:
plot(cos, label = "f1")
plot!(sin, label = "f2")
plot!(xlims = (-5, 5), title = "title", xlabel = "x", ylabel = "y")
Out[25]:
No description has been provided for this image
  • Les animations peuvent être créées de la façon suivante :
In [26]:
anim = Animation()
for θ in LinRange(0, 2π, 50)
    plot(x -> cos(x - θ), label=nothing)
    title!("θ = $(round(θ, digits=3))")
    frame(anim)
end
gif(anim, fps=10)

[ Info: Saved animation to /home/urbain/work/teaching/casc_private/notebooks/2026/source/week1/tmp.gif
Out[26]:
No description has been provided for this image
  • La bibliothèque Plots définit la macro @animate pour simplifier la création d'animations :
In [27]:
anim = @animate for θ in LinRange(0, 2π, 50)
    plot(x -> cos(x - θ), label=nothing)
    title!("θ = $(round(θ, digits=3))")
end
gif(anim, fps=10)

[ Info: Saved animation to /home/urbain/work/teaching/casc_private/notebooks/2026/source/week1/tmp.gif
Out[27]:
No description has been provided for this image

[Exercice 3]¶

Créer une animation illustrant le mouvement d'une balle de golf dont la position dans le plan $x,y$ est donnée par $$ \begin{aligned} x(t) &= u_0 t, \\ y(t) &= v_0 t - g \frac{t^2}{2}, \end{aligned} $$ Utiliser les paramètres fournis ainsi que la fonction scatter, et rendre l'animation aussi soignée que possible.

In [28]:
u₀ = 50
v₀ = 50
g = 9.81
### BEGIN SOLUTION
X(t) = u₀*t
Y(t) = v₀*t - g*t^2/2
tf = 2v₀/g
nframes = 100

anim = @animate for t in LinRange(0, tf, nframes)
    scatter([X(t)], [Y(t)], label=nothing, markersize=10)
    title!("t = $(round(t, digits=3))")
    xlims!(0, X(tf))
    ylims!(0, Y(tf/2) * 1.1)
end
gif(anim, fps=nframes/tf)
### END SOLUTION

[ Info: Saved animation to /home/urbain/work/teaching/casc_private/notebooks/2026/source/week1/tmp.gif
Out[28]:
No description has been provided for this image

Formats à virgule flottante¶

Seul un nombre fini de réels peut être représenté exactement sur un ordinateur. Par conséquent, certains réels sont représentables exactement et d'autres non. La norme IEEE 754, à laquelle se conforment pratiquement tous les langages de programmation, spécifie un petit nombre de formats à virgule flottante. L'ensemble des réels représentables dans chaque format est de la forme $$ \begin{align*} F(p, E_{\min}, E_{\max}) = \Bigl\{ & (-1)^s 2^E (b_0. b_1 b_2 \dots b_{p-1})_2 \colon \\ & \qquad s \in \{0, 1\}, b_i \in \{0, 1\} \, \text{et} \, E_{\min} \leq E \leq E_{\max} \Bigr\}. \end{align*} $$ Ici $(b_0.b_1 b_2 \dots b_{p-1})_2$ est un nombre en représentation binaire, c'est-à-dire $$ (b_0.b_1 b_2 \dots b_{p-1})_2 = b_0 + \frac{1}{2} b_1 + \dotsc + \frac{1}{2^{p-1}} b_{p-1}. $$ Pour un élément $x \in F(p, E_{\min}, E_{\max})$, $(b_0. b_1 b_2 \dots b_{p-1})_2$ est appelée la mantisse, $s$ est le signe et $E$ est l'exposant. Trois paramètres apparaissent dans la définition de l'ensemble ci-dessus :

  • $p$ est le nombre de bits significatifs (appelé aussi la précision),
  • $E_{\min}$ est l'exposant minimum,
  • $E_{\max}$ est l'exposant maximum.

À partir de la précision, l'epsilon machine associé au format est défini comme $$ \varepsilon = 2^{-(p-1)} $$

Les paramètres des formats les plus couramment utilisés spécifiés par la norme sont les suivants :

Paramètre Demi-précision (Float16) Simple précision (Float32) Double précision (Float64)
$p$ 11 24 53
$E_{\min}$ -14 -126 -1022
$E_{\max}$ 15 127 1023

En Julia, les formats à virgule flottante demi, simple et double précision sont donnés par Float16, Float32 et Float64, le format par défaut étant Float64.

In [29]:
a = √2
b = Float32(√2)
c = Float16(√2)

@show a
@show b
@show c

@show typeof(a)
@show typeof(b)
@show typeof(c);

a = 1.4142135623730951
b = 1.4142135f0
c = Float16(1.414)
typeof(a) = Float64
typeof(b) = Float32
typeof(c) = Float16

L'epsilon machine doit être compris comme une mesure de l'espacement relatif entre un nombre à virgule flottante et le suivant, au sein d'un format. En effet, le nombre 1 est représentable exactement dans tous les formats, et le prochain nombre représentable est précisément $1 + \varepsilon$

In [30]:
a = Float16(1.0)
@show nextfloat(a) - a
@show eps(Float16);

nextfloat(a) - a = Float16(0.000977)
eps(Float16) = Float16(0.000977)

De même, le prochain nombre représentable après $2$ est $2 + 2\varepsilon$ et le prochain nombre représentable après $4$ est $4 + 4\varepsilon$ :

In [31]:
ε = eps(Float64)
@show nextfloat(2.0) - 2.0
@show 2ε;
@show nextfloat(4.0) - 4.0
@show 4ε;

nextfloat(2.0) - 2.0 = 4.440892098500626e-16
2ε = 4.440892098500626e-16
nextfloat(4.0) - 4.0 = 8.881784197001252e-16
4ε = 8.881784197001252e-16

Comme le suggèrent leurs noms, un nombre Float16 est stocké sur 16 bits, un Float32 sur 32 bits, et un nombre Float64 sur 64 bits. Parmi ces bits, un est utilisé pour stocker le signe, et les autres pour stocker la mantisse et l'exposant. De plus, quelques combinaisons de bits sont réservées pour représenter les valeurs spéciales Inf, -Inf et NaN, ce dernier étant l'acronyme de « not a number » (pas un nombre).

Un octet, abrégé par un B majuscule, est une unité composée de 8 bits. C'est l'unité généralement utilisée pour décrire la quantité de mémoire disponible dans les produits commerciaux ; par exemple, un ordinateur portable dispose typiquement de 8 Go ou 16 Go de mémoire vive. Il est souvent utile d'avoir une idée approximative de l'espace occupé en mémoire par les objets mathématiques. Par exemple, une matrice 1000 × 1000 de nombres Float64 nécessite $8 \times 10^6$ octets pour être stockée, et une matrice 10000 × 10000 de nombres Float64 requiert près de 800 Mo de mémoire vive.

In [32]:
A = zeros(10_000, 10_000)
@show Base.summarysize(A);

Base.summarysize(A) = 800000048
























Pour approfondir notre compréhension des formats à virgule flottante,
traçons sur un même graphique l'espacement entre
les nombres successifs Float16, Float32 et Float64
dans l'intervalle $[1, 10^4]$. On constate que la densité de
nombres à virgule flottante représentables autour de $x$ diminue
avec $x$.


















In [33]:






n = 1000
x_16 = (^).(10, LinRange{Float16}(0, 4, n))
x_32 = (^).(10, LinRange{Float32}(0, 4, n))
x_64 = (^).(10, LinRange{Float64}(0, 4, n))
spacing16 = nextfloat.(x_16) - x_16
spacing32 = nextfloat.(x_32) - x_32
spacing64 = nextfloat.(x_64) - x_64
plot(x_16, spacing16, label="Float16", )
plot!(x_32, spacing32, label="Float32")
plot!(x_64, spacing64, label="Float64")
plot!(xscale=:log10, yscale=:log10,
      xlabel="x", ylabel="Δx", legend=:bottomright)





















Out[33]:



No description has been provided for this image























[Exercice 4]¶
Pour approfondir votre compréhension,


    

  • Lire la documentation de la fonction nextfloat.
    • 

  • En utilisant cette fonction,
tracer l'espacement relatif, normalisé par l'epsilon machine du format Float32,
entre les nombres Float32 successifs dans l'intervalle $[10^{-2}, 10^2]$.
Plus précisément, tracer la fonction
$$
x \mapsto \frac{\text{nextfloat}(x) - x}{x \varepsilon}
$$
    • 



Utiliser une échelle logarithmique pour l'axe $x$ uniquement.
Il peut être utile d'utiliser LinRange{Type}(a, b, n)
pour créer un vecteur de $n$ nombres de type Type équidistants entre $a$ et $b$.


















In [34]:






### BEGIN SOLUTION
x = (^).(10, LinRange{Float32}(-2, 2, 10000))
spacing = (nextfloat.(x) - x) ./ (x*eps(Float32))
plot(x, spacing, label="Relative spacing")
plot!(xscale=:log2, xlabel="x", ylabel="Δx/x (in ε units)")
### END SOLUTION





















Out[34]:



No description has been provided for this image























Abordons maintenant une question naturelle découlant de cette discussion :
que se passe-t-il si le résultat d'une opération mathématique n'est pas exactement représentable dans le format à virgule flottante ?
La réponse est simple : le résultat est arrondi au nombre représentable le plus proche.
Cela explique les résultats surprenants suivants :


















In [35]:






@show 1 + eps()/3      # the result is 1
@show 1 + 2*eps()/3    # the result is 1 + ε
@show .1 + .2 == .3    # this is False

























1 + eps() / 3 = 1.0
1 + (2 * eps()) / 3 = 1.0000000000000002
0.1 + 0.2 == 0.3 = false









Out[35]:




false
























Les erreurs causées par ce processus d'arrondi sont appelées erreurs d'arrondi.
Elles sont particulièrement importantes pour le calcul des dérivées.
En effet, supposons que l'on souhaite calculer,
pour une fonction différentiable $f\colon \mathbb R \to \mathbb R$, une approximation de $f'(1)$.
L'approche la plus naturelle est de calculer pour un petit $\delta > 0$
$$
d(\delta) = \frac{f(1+ \delta) - f(1)}{\delta}.
$$
En arithmétique exacte, cela fournirait une approximation de $f'(1)$ avec une erreur en $\mathcal O(\delta)$,
en supposant $f''(1) \neq 0$.
En arithmétique d'ordinateur, cependant, une erreur supplémentaire apparaît en raison des arrondis,
dont on peut montrer qu'elle est de l'ordre de $\mathcal O(\varepsilon/\delta)$.
Le meilleur choix est de prendre $\delta \approx \sqrt{\varepsilon}$,
de sorte que les erreurs soient approximativement du même ordre de grandeur.


















In [36]:






f(x) = exp(x)
d(δ) = (f(1+δ) - f(1))/δ
δs = 10 .^(0:-.1:-17)
err = abs.(d.(δs) .- exp(1))
plot(δs, err, xscale=:log10, yscale=:log10, label=L"|d(\delta) - f'(1)|", size=(900,450))
plot!([eps()], seriestype=:vline, label = L"\varepsilon")
plot!([sqrt(eps())], seriestype=:vline, label = L"\sqrt{\varepsilon}")





















Out[36]:



No description has been provided for this image























Pour des $\delta$ très petits de l'ordre de l'epsilon machine $\varepsilon$,
l'approximation obtenue est complètement fausse.
Comprenez-vous les observations suivantes ?


















In [37]:






δ = eps()/4; @show (exp(δ) - exp(0)) / δ
δ = eps()/2; @show (exp(δ) - exp(0)) / δ
δ = eps(); @show (exp(δ) - exp(0)) / δ;

























(exp(δ) - exp(0)) / δ = 0.0
(exp(δ) - exp(0)) / δ = 2.0
(exp(δ) - exp(0)) / δ = 1.0

























Remarque. Dans de nombreux contextes, il est utile d'utiliser des nombres à virgule flottante avec plus de précision que les types standard,
afin de réduire les erreurs d'arrondi et d'observer l'erreur intrinsèque,
liée au caractère inexact de nombreux algorithmes même en arithmétique exacte.
Le type BigFloat en Julia est un format à précision arbitraire.
Le nombre de bits significatifs correspondant peut être défini à l'aide de la fonction setprecision.


















In [38]:






@show Float64(π)    

setprecision(100, base=10)  # We want 100 digits in base 10
BigFloat(π)

























Float64(π) = 3.141592653589793









Out[38]:




3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706807
























Manipulation des tableaux¶






















    

  • Définition explicite simple
    • 



















In [39]:






u = [1, 2, 3]  # defines a column Vector
display(u)
@show typeof(u)

























3-element Vector{Int64}:
 1
 2
 3












typeof(u) = Vector{Int64}




















Out[39]:




Vector{Int64} (alias for Array{Int64, 1})
























    

  • Remarquer la différence avec
    • 



















In [40]:






u = [1 2 3] # defines a Matrix of 1 line
display(u)
@show typeof(u)

























1×3 Matrix{Int64}:
 1  2  3












typeof(u) = Matrix{Int64}




















Out[40]:




Matrix{Int64} (alias for Array{Int64, 2})
























    

  • Définition explicite d'une matrice
    • 



















In [41]:






A = [1 2 3
     4 5 6
     7 8 9]





















Out[41]:




3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9




















In [42]:






A = [1 2 3; 4 5 6; 7 8 9]





















Out[42]:




3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9
























    

  • Extraction d'une colonne ou d'une ligne
    • 





Note : en Julia, les indices commencent à 1 comme en Fortran (et non à 0 comme en Python ou C)




















In [43]:






A[1,:]





















Out[43]:




3-element Vector{Int64}:
 1
 2
 3




















In [44]:






A[:,1]





















Out[44]:




3-element Vector{Int64}:
 1
 4
 7
























    

  • : est un raccourci pour begin:end ou 1:end
    • 



















In [45]:






A[2,1:end]





















Out[45]:




3-element Vector{Int64}:
 4
 5
 6




















In [46]:






A[begin:end,2]





















Out[46]:




3-element Vector{Int64}:
 2
 5
 8
























    

  • Construction d'un tableau par syntaxe de compréhension
    • 



















In [47]:






A = [10i+j for i in 0:2, j in 0:3]
display(A)
print("\nA = $A\n\n")
@show A
@show A[1,1]
@show A[1,2]
@show length(A)
@show size(A)
@show size(A, 1)
@show size(A, 2)
@show eltype(A)
;

























3×4 Matrix{Int64}:
  0   1   2   3
 10  11  12  13
 20  21  22  23












A = [0 1 2 3; 10 11 12 13; 20 21 22 23]

A = 












[0 1 2 3; 10 11 12 13; 20 21 22 23]
A[1, 1] = 0
A[1, 2] = 1
length(A) = 12
size(A) = (3, 4)
size(A, 1) = 3
size(A, 2) = 4
eltype(A) = Int64

























    

  • Matrice adjointe
    • 



















In [48]:






A = [1 2 3
     4 5 6]
A'





















Out[48]:




3×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6




















In [49]:






display(A'A)
display(A*A')

























3×3 Matrix{Int64}:
 17  22  27
 22  29  36
 27  36  45












2×2 Matrix{Int64}:
 14  32
 32  77




















In [50]:






A = [2+3im 4+2im 1+3im
     1+2im 2+3im 5+6im]
A'





















Out[50]:




3×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 2-3im  1-2im
 4-2im  2-3im
 1-3im  5-6im




















In [51]:






display(A'A)
display(A*A')

























3×3 Matrix{Complex{Int64}}:
 18+0im  22-9im  28-1im
 22+9im  33+0im  38+7im
 28+1im  38-7im  71+0im












2×2 Matrix{Complex{Int64}}:
 43+0im  45+0im
 45+0im  79+0im
























    

  • Concaténation de tableaux
    • 



















In [52]:






A = [1 2 3
     4 5 6]
u = [10, 20, 30]
v = [100, 200]
;





















In [53]:






[A v]





















Out[53]:




2×4 Matrix{Int64}:
 1  2  3  100
 4  5  6  200




















In [54]:






[A; u']





















Out[54]:




3×3 Matrix{Int64}:
  1   2   3
  4   5   6
 10  20  30




















In [55]:






[A A]





















Out[55]:




2×6 Matrix{Int64}:
 1  2  3  1  2  3
 4  5  6  4  5  6




















In [56]:






[A; A]





















Out[56]:




4×3 Matrix{Int64}:
 1  2  3
 4  5  6
 1  2  3
 4  5  6




















In [57]:






hcat(A, A, A)





















Out[57]:




2×9 Matrix{Int64}:
 1  2  3  1  2  3  1  2  3
 4  5  6  4  5  6  4  5  6




















In [58]:






vcat(A, A, A)





















Out[58]:




6×3 Matrix{Int64}:
 1  2  3
 4  5  6
 1  2  3
 4  5  6
 1  2  3
 4  5  6
























    

  • Les composantes sont stockées par colonnes et peuvent être accédées par un indice simple même pour des tableaux multidimensionnels
    • 

  • eachindex et CartesianIndices permettent d'itérer sur tous les indices d'un tableau (indice simple pour le premier, multi-indice pour le second)
    • 

  • Un CartesianIndex peut être converti en Tuple
    • 

  • ... est l'opérateur splat (voir la documentation)
    • 



















In [59]:






A = [1 2 3; 4 5 6]

for k ∈ eachindex(A)
    println("A[$k] = ", A[k])
end
println()

for k ∈ CartesianIndices(A)
    println("A[$k] = A$(Tuple(k)...) = ", A[k])
end
println()

for i in 1:size(A, 1)
    for j in 1:size(A, 2)
        println("A[$i, $j] = ", A[i, j])
    end
end
println()

for i in 1:size(A, 1), j in 1:size(A, 2)
    println("A[$i, $j] = ", A[i, j])
end

























A[1] = 1












A[2] = 4
A[3] = 2
A[4] = 5
A[5] = 3
A[6] = 6

A[CartesianIndex(1, 1)] = A11 = 1
A[CartesianIndex(2, 1)] = A21 = 4
A[CartesianIndex(1, 2)] = A12 = 2
A[CartesianIndex(2, 2)] = A22 = 5
A[CartesianIndex(1, 3)] = A13 = 3
A[CartesianIndex(2, 3)] = A23 = 6

A[1, 1] = 1












A[1, 2] = 2
A[1, 3] = 3
A[2, 1] = 4
A[2, 2] = 5
A[2, 3] = 6

A[1, 1] = 1
A[1, 2] = 2
A[1, 3] = 3
A[2, 1] = 4
A[2, 2] = 5
A[2, 3] = 6

























    

  • Une matrice (tableau d'ordre 2) n'est pas un vecteur de vecteurs
    • 



















In [60]:






B = [[10i+j for i in 0:2] for j in 0:3]
display(B)
@show B isa Matrix 
;

























4-element Vector{Vector{Int64}}:
 [0, 10, 20]
 [1, 11, 21]
 [2, 12, 22]
 [3, 13, 23]












B isa Matrix = false




































    

  • Les tableaux sont mutables (si le type est cohérent), les tuples ne le sont pas !
    • 



















In [61]:






v = [1, 2]
@show typeof(v)
@show eltype(v)
@show v
v[1] = 3
@show v
# v[1] = 1.5 # error since v is of type Vector{Int64}
# @show v
t = (1, 2)
@show t[1] ;
# t[1] = 3 # error since t is not mutable

























typeof(v) = Vector{Int64}
eltype(v) = Int64
v = [1, 2]
v = [3, 2]
t[1] = 1

























    

  • Définition d'un tableau vide et construction a posteriori
    • 





Remarquer la différence entre push! et append! (il est rappelé que pour les deux, ! exprime que par convention, et non par obligation, l'argument est modifié)




















In [62]:






v = [] # void vector containing any type (not optimal!)
@show typeof(v)
push!(v, 1)
append!(v, [2.5, 3.6]) # pushes each item of the vector [2.5, 3.6]
@show v
push!(v, [4.1, 5.2]) # pushes a vector as new item
@show v ;
@show v[end]
@show v[end-1] ;

























typeof(v) = Vector{Any}
v = 












Any[1, 2.5, 3.6]
v = Any[1, 2.5, 3.6, [4.1, 5.2]]
v[end] = [4.1, 5.2]
v[end - 1] = 3.6

























    

  • Pour de meilleures performances, il est préférable, si possible, de préciser le type des éléments du vecteur lors de la construction
    • 



















In [63]:






v = Float64[]
@show typeof(v)
push!(v, 1)
append!(v, [2//3, π]) # pushes each item of the vector [2//3, π] and convert on-the-fly to Float64
@show v
# push!(v, [4.3, 5.7]) # throws an error since [4.3, 5.7] is not of type Float64
;

























typeof(v) = Vector{Float64}
v = [1.0, 0.6666666666666666, 3.141592653589793]

























    

  • Référence d'un tableau
    • 



















In [64]:






A = [10i+j for i in 1:2, j in 1:3]
@show A
B = A # B is a reference to A (both point to the same memory location)
@show B == A # same content
@show B === A # same memory location
B .= 99 # .= assigns all the coefficients of the array to 99
@show B
@show A # A is of course also affected
;

























A = [11 12 13; 21 22 23]
B == A = true
B === A = true
B = [99 99 99; 99 99 99]
A = [99 99 99; 99 99 99]

























    

  • Copie d'un tableau
    • 



















In [65]:






A = [10i+j for i in 1:2, j in 1:3]
@show A
B = copy(A) # B and A are two distinct arrays
@show B == A # same content
@show B === A # not the same memory location
B .= 99 # .= assigns all the coefficients of the array to 99
@show B
@show A # A unaffected
;

























A = [11 12 13; 21 22 23]
B == A = true
B === A = false
B = [99 99 99; 99 99 99]
A = [11 12 13; 21 22 23]

























    

  • Tranche en tant que l-value ≡ référence
    • 



















In [66]:






A = [10i+j for i in 1:6, j in 1:8]
display(A)
A[1:2:end, 2:2:end] .= 0 # components of odd line and even column are set to 0
display(A)
;

























6×8 Matrix{Int64}:
 11  12  13  14  15  16  17  18
 21  22  23  24  25  26  27  28
 31  32  33  34  35  36  37  38
 41  42  43  44  45  46  47  48
 51  52  53  54  55  56  57  58
 61  62  63  64  65  66  67  68












6×8 Matrix{Int64}:
 11   0  13   0  15   0  17   0
 21  22  23  24  25  26  27  28
 31   0  33   0  35   0  37   0
 41  42  43  44  45  46  47  48
 51   0  53   0  55   0  57   0
 61  62  63  64  65  66  67  68
























    

  • Tranche en tant que r-value ≡ copie
    • 



















In [67]:






A = [10i+j for i in 1:6, j in 1:8]
display(A)
B = A[1:2:end, 2:2:end] # B is a copy of the submatrix
display(B)
B .= 0
display(B)
display(A) # A remains unaffected
;

























6×8 Matrix{Int64}:
 11  12  13  14  15  16  17  18
 21  22  23  24  25  26  27  28
 31  32  33  34  35  36  37  38
 41  42  43  44  45  46  47  48
 51  52  53  54  55  56  57  58
 61  62  63  64  65  66  67  68












3×4 Matrix{Int64}:
 12  14  16  18
 32  34  36  38
 52  54  56  58












3×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0












6×8 Matrix{Int64}:
 11  12  13  14  15  16  17  18
 21  22  23  24  25  26  27  28
 31  32  33  34  35  36  37  38
 41  42  43  44  45  46  47  48
 51  52  53  54  55  56  57  58
 61  62  63  64  65  66  67  68
























    

  • Référence en tant que r-value via view
    • 



















In [68]:






A = [10i+j for i in 1:6, j in 1:8]
display(A)
B = view(A, 1:2:size(A,1), 2:2:size(A,2)) # B is a view of the submatrix (reference not a copy)
display(B)
B .= 0
display(B)
display(A) # A is modified
;

























6×8 Matrix{Int64}:
 11  12  13  14  15  16  17  18
 21  22  23  24  25  26  27  28
 31  32  33  34  35  36  37  38
 41  42  43  44  45  46  47  48
 51  52  53  54  55  56  57  58
 61  62  63  64  65  66  67  68












3×4 view(::Matrix{Int64}, 1:2:5, 2:2:8) with eltype Int64:
 12  14  16  18
 32  34  36  38
 52  54  56  58












3×4 view(::Matrix{Int64}, 1:2:5, 2:2:8) with eltype Int64:
 0  0  0  0
 0  0  0  0
 0  0  0  0












6×8 Matrix{Int64}:
 11   0  13   0  15   0  17   0
 21  22  23  24  25  26  27  28
 31   0  33   0  35   0  37   0
 41  42  43  44  45  46  47  48
 51   0  53   0  55   0  57   0
 61  62  63  64  65  66  67  68
























    

  • Autres façons de définir des tranches par sélection d'indices
    • 



















In [69]:






A = [10i+j for i in 1:6, j in 1:8]
A[isodd.(1:end), iseven.(1:end)] .= 0 ; display(A)
A = [10i+j for i in 1:6, j in 1:8]
A[findall(isodd,A)] .= 0 ; display(A) ;

























6×8 Matrix{Int64}:
 11   0  13   0  15   0  17   0
 21  22  23  24  25  26  27  28
 31   0  33   0  35   0  37   0
 41  42  43  44  45  46  47  48
 51   0  53   0  55   0  57   0
 61  62  63  64  65  66  67  68












6×8 Matrix{Int64}:
 0  12  0  14  0  16  0  18
 0  22  0  24  0  26  0  28
 0  32  0  34  0  36  0  38
 0  42  0  44  0  46  0  48
 0  52  0  54  0  56  0  58
 0  62  0  64  0  66  0  68
























    

  • Expliquer les lignes suivantes
    • 



















In [70]:






A = [10i+j for i in 1:6, j in 1:8]
B = view(A, collect(1:size(A,1)) .% 3 .== 0, collect(1:size(A,2)) .% 5 .!= 0)
B .= 0
A





















Out[70]:




6×8 Matrix{Int64}:
 11  12  13  14  15  16  17  18
 21  22  23  24  25  26  27  28
  0   0   0   0  35   0   0   0
 41  42  43  44  45  46  47  48
 51  52  53  54  55  56  57  58
  0   0   0   0  65   0   0   0




















In [71]:






A = rand(5, 5) ; display(A)
A[findall(x -> x < 1/2, A)] .= 0
A

























5×5 Matrix{Float64}:
 0.629737  0.48251   0.697925  0.68341     0.902767
 0.970592  0.425703  0.938402  0.512732    0.391903
 0.572262  0.276557  0.880927  0.00621013  0.289465
 0.651025  0.687269  0.978548  0.571009    0.562262
 0.157937  0.27258   0.713588  0.376401    0.0146206








Out[71]:




5×5 Matrix{Float64}:
 0.629737  0.0       0.697925  0.68341   0.902767
 0.970592  0.0       0.938402  0.512732  0.0
 0.572262  0.0       0.880927  0.0       0.0
 0.651025  0.687269  0.978548  0.571009  0.562262
 0.0       0.0       0.713588  0.0       0.0




















In [72]:






n, m = 7, 9
B = rand(n, m)
u = rand(m)
display(B*u)
subi = collect(1:n) .% 3 .!= 0 ; @show subi
subj = collect(1:m) .% 3 .== 0 ; @show subj
B[subi,subj]*u[subj]

























7-element Vector{Float64}:
 2.1784480846131173
 2.2422262996043503
 2.0366621622328287
 2.0622759962259187
 2.2764233260826416
 1.8064111486377548
 0.927613589982037












subi = Bool[1, 1, 0, 1, 1, 0, 1]












subj = Bool[0, 0, 1, 0, 0, 1, 0, 0, 1]









Out[72]:




5-element Vector{Float64}:
 0.6338022054600351
 0.7072326593247753
 0.8209544084710862
 0.3885537023137302
 0.2351861237828015
























[Exercice 5]¶
Écrire une fonction qui calcule la comatrice d'une matrice explicitement à partir du déterminant de sous-matrices.


On rappelle que la composante $C_{ij}$ de la comatrice de $A$ est obtenue comme le déterminant de la sous-matrice obtenue à partir de $A$ en supprimant sa $i^\textrm{ème}$ ligne et sa $j^\textrm{ème}$ colonne, multiplié par $(-1)^{i+j}.$ Le déterminant d'une petite matrice dense est calculé par la fonction det disponible dans la bibliothèque LinearAlgebra.


















In [73]:






using LinearAlgebra
function comatrix(A)
    C = zero(A)
    ### BEGIN SOLUTION
    for ij ∈ CartesianIndices(A)
        i, j = Tuple(ij)
        C[ij] = (-1)^(i+j)*det(A[1:end .!= i, 1:end .!= j])
    end
    ### END SOLUTION
    return C
end





















Out[73]:




comatrix (generic function with 1 method)




















In [74]:






n = 10
A = rand(n, n)
C = comatrix(A)
display(C'A)
@assert C'A ≈ det(A)*LinearAlgebra.I

























10×10 Matrix{Float64}:
 -0.0118234     1.07504e-17  -7.68942e-18  …   4.51631e-18  -2.87112e-18
 -1.22014e-18  -0.0118234    -1.72461e-18     -1.94408e-18  -1.4214e-18
 -2.7005e-18   -8.47232e-19  -0.0118234        1.18278e-18  -1.7136e-18
  5.02304e-20  -1.16077e-18  -1.22284e-18      6.63492e-19   1.19226e-18
 -3.95394e-18  -6.12742e-18  -3.54564e-18     -3.30749e-18  -1.24493e-18
  1.2887e-18    7.9194e-18    1.71291e-18  …   3.59662e-18   7.14169e-19
  5.14955e-18   8.34019e-18   9.2066e-18       3.19522e-18   2.36739e-18
 -6.65131e-18  -1.11913e-17  -7.39803e-18     -1.83992e-18  -4.25048e-18
 -7.72255e-18  -9.44905e-18  -8.38531e-18     -0.0118234    -7.52197e-18
 -2.88097e-17  -2.6799e-17   -3.49374e-17     -1.67336e-17  -0.0118234
























Création de structures personnalisées¶






















Julia n'est pas un langage orienté objet.
Il supporte à la place la distribution multiple (multiple dispatch),
dont la motivation est exposée éloquemment dans la documentation de Julia :




The choice of which method to execute when a function is applied is called dispatch. Julia allows the dispatch process to choose which of a function's methods to call based on the number of arguments given, and on the types of all of the function's arguments. This is different from traditional object-oriented languages, where dispatch occurs based only on the first argument, which often has a special argument syntax, and is sometimes implied rather than explicitly written as an argument. Using all of a function's arguments to choose which method should be invoked, rather than just the first, is known as multiple dispatch. Multiple dispatch is particularly useful for mathematical code, where it makes little sense to artificially deem the operations to "belong" to one argument more than any of the others: does the addition operation in x + y belong to x any more than it does to y? The implementation of a mathematical operator generally depends on the types of all of its arguments. Even beyond mathematical operations, however, multiple dispatch ends up being a powerful and convenient paradigm for structuring and organizing programs.




















In [75]:






struct Dog
  name::String
end

struct Cat
  name::String
end

greet(a::Dog, b::Dog) = "$(a.name): Wouf wouf $(b.name), content de te revoir mon vieux"
greet(a::Dog, b::Cat) = "$(a.name): Wouf wouf wouf wouf !"
greet(a::Cat, b::Dog) = "$(a.name) s'éloigne discrètement de $(b.name)"
greet(a::Cat, b::Cat) = "$(a.name) commence à ronronner"

milou = Dog("Milou")
rantanplan = Dog("Rantanplan")
miaouss = Cat("Miaouss")
mrMoustache = Cat("MrMoustache")

@show greet(milou, rantanplan)
@show greet(milou, mrMoustache)
@show greet(mrMoustache, miaouss)
@show greet(mrMoustache, rantanplan);

























greet(milou, rantanplan) = "Milou: Wouf wouf Rantanplan, content de te revoir mon vieux"
greet(milou, mrMoustache) = "Milou: Wouf wouf wouf wouf !"
greet(mrMoustache, miaouss) = "MrMoustache commence à ronronner"
greet(mrMoustache, rantanplan) = "MrMoustache s'éloigne discrètement de Rantanplan"

























Pour conclure cette séance,
illustrons la puissance de la distribution multiple avec deux autres exemples.
Premièrement, créons un type spécifique pour la matrice identité,
et redéfinissons les opérations + et * sur les matrices pour qu'elles fonctionnent avec ce nouveau type.


















In [76]:






struct IdentityMatrix end

# This is required to redefine + and *
import Base: +,*

function +(A::Matrix, Id::IdentityMatrix)
    result = copy(A)
    for i in 1:size(A, 1)
        result[i, i] += 1
    end
    return result
end

+(Id::IdentityMatrix, A::Matrix) = +(A, Id)
*(A::Matrix, Id::IdentityMatrix) = copy(A)
*(Id::IdentityMatrix, A::Matrix) = copy(A)

A = ones(3, 3)
Id = IdentityMatrix()

# We can now write operations such as
@show A
@show A + Id
@show Id + A
@show A * Id;

























A = [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]
A + Id = 












[2.0 1.0 1.0; 1.0 2.0 1.0; 1.0 1.0 2.0]
Id + A = [2.0 1.0 1.0; 1.0 2.0 1.0; 1.0 1.0 2.0]
A * Id = [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0]

























Remarquer que nous n'avons pas eu besoin de modifier l'implémentation interne des matrices.
On notera également qu'avec notre implémentation de la matrice identité,
le produit A*I est pratiquement gratuit en temps de calcul.
Si la matrice identité avait été stockée comme une matrice pleine,
beaucoup de puissance de calcul aurait été gaspillée dans la multiplication matricielle.


















In [77]:






n = 5000
A = rand(n, n)
I1 = IdentityMatrix()
I2 = [(i == j)*1.0 for i in 1:n, j in 1:n]

@time A*I1;
@time A*I2;

























  0.133289 seconds (3 allocations: 190.738 MiB, 34.73% gc time)
  2.961209 seconds (114 allocations: 190.743 MiB, 5.11% gc time, 0.35% compilation time)

























[Exercice 6]¶
L'objectif de cet exercice est de créer une structure pour représenter les nombres duaux.
Les nombres duaux sont une construction mathématique qui étend
les réels en introduisant un nouvel élément, souvent
noté $\varepsilon$, tel que $\varepsilon^2 = 0$.
Un nombre dual est composé d'une partie réelle et d'une partie duale,
où la partie duale se comporte comme une quantité infinitésimale.
Ils seront très utiles plus loin dans le cours,
lorsque nous étudierons la différentiation automatique, une technique
pour calculer efficacement des dérivées numériquement.
Compléter la définition de structure ci-dessous en ajoutant les opérations -, * et / aux nombres duaux.


















In [78]:






import Base: +, -, *, /, show

struct Dual
    real::Float64
    dual::Float64
end

+(x::Dual, y::Dual) = Dual(x.real + y.real, x.dual + y.dual)
### BEGIN SOLUTION
*(x::Dual, y::Dual) = Dual(x.real*y.real, (x.dual*y.real + x.real*y.dual))
-(x::Dual, y::Dual) = Dual(x.real - y.real, x.dual - y.dual)
/(x::Dual, y::Dual) = Dual(x.real/y.real, (y.real*x.dual - x.real*y.dual)/y.real^2)
### END SOLUTION

show(io::IO, x::Dual) = print(io, x.real, x.dual ≥ 0 ? " + $(x.dual)" : " - $(-x.dual)", " ε")

# Example usage
a = Dual(3.0, 2.0)  # 3 + 2ε
b = Dual(1.0, 4.0)  # 1 + 4ε

@show a + b
@show a - b
@show a * b
@show a / b;

























a + b = 4.0 + 6.0 ε
a - b = 2.0 - 2.0 ε
a * b = 3.0 + 14.0 ε
a / b = 3.0 - 10.0 ε





















In [79]:






@assert a + b == Dual(4, 6)
@assert a - b == Dual(2, -2)
@assert a * b == Dual(3, 14)
@assert a / b == Dual(3, -10)