Examen final¶

  • Ce notebook est à soumettre sur Educnet avant 11h30.

  • L’examen comporte trois exercices indépendants. Dans chaque exercice, les cellules peuvent éventuellement dependre des cellules précèdentes.

  • Afin de faciliter l'évaluation de votre code, ne pas changer les signatures des fonctions à implémenter.

  • La cellule ci-dessous importe les bibliothèques utilisées dans ce notebook. Ce sont a priori des bibliothèques déjà présentes dans l'installation de base de Julia ou déjà installées en TD. Si une ou plusieurs d'entre elles manquent sur votre machine, vous êtes invités à les installer au préalable dans le gestionnaire de bibliothèques d'une console.

In [1]:
using ForwardDiff
using LaTeXStrings
using LinearAlgebra
using Plots
using Polynomials
using Random

Plots.default(titlefontsize=12,
              xlabelfontsize=10,
              ylabelfontsize=10,
              legendfontsize=10,
              xtickfontsize=10,
              ytickfontsize=10)

[Exercice 1] Intégration par interpolation¶

Le but de cet exercice est de calculer l'intégrale de la fonction $f \colon [-5, \,5] \to \mathbb{R} $ donnée par

$$ f (x) = \log(1 + |x| - x^3 + x^4). $$

Pour cela, nous vous proposons d'interpoler cette fonction par un polynome sur l'intervalle $[-5, \, 5]$, puis d'intégrer ce dernier.

  1. Ecrire une fonction interpolate(func, d, a) qui retourne le polynôme $p$ (de type Polynomial de la bibliothèque Polynomials.jl) résultant de l'interpolation polyonmiale de la fonction func sur l'intervalle $[-a, a]$ de degré $d$. Pour cela, vous pourrez utiliser la fonction fit de la bibliothèque Polynomials.jl, et choisir les nœuds d'interpolation de la manière qui vous semble la mieux adaptée.
In [2]:
f(x) = log(1 + abs(x) - x^3 + x^4)

function interpolate(func, d, a)
    ### BEGIN SOLUTION
    chebyshev_nodes(k, n) = cospi((2k - 1) / (2n))
    x = a * chebyshev_nodes.(1:d+1, d+1)
    y = func.(x)
    p = fit(x, y)
    ### END SOLUTION
    return p
end;
In [3]:
@assert typeof(interpolate(f, 20, 5)) <: Polynomial
@assert length(interpolate(f, 20, 5).coeffs) == 21
@assert begin xs = LinRange(-1, 1, 200); interpolate(cos, 20, 1).(xs) .- cos.(xs) |> maximum end ≤ 1e-14
@assert begin xs = LinRange(-1, 1, 200); func = x -> x^3; p = interpolate(func, 20, 1); sum(abs2, p.(xs) - func.(xs)) end ≤ 1e-14
  1. Illustrer sur un même graphique la fonction f définie ci-dessus, ainsi que son interpolation de degré 20.
In [4]:
# ### BEGIN SOLUTION
x = LinRange(-5, 5, 201)
p = interpolate(f, 20, 5)
plot(x, f, label=L"Function $f$")
plot!(x, p.(x), label="Interpolation")
# ### END SOLUTION
Out[4]:
No description has been provided for this image

  1. Écrire une fonction integrate_via_interp(func, d, a) qui retourne une approximation de l'intégrale $$ \int_{-a}^a \mathrm{func}(x) \, \mathrm{d} x\, , $$ en calculant d'abord une interpolation de degré d de la fonction passée en argument, puis en intégrant cette dernière de manière exacte. Pour cela, on pourra utiliser la fonction integrate de la bibliothèque Polynomials.jl, qui permet de calculer la primitve d'un polynôme:

        P = integrate(p)
In [5]:
function integrate_via_interp(func, d, a)
    ### BEGIN SOLUTION
    p = interpolate(func, d, a)
    P = integrate(p)
    return P(a) - P(-a)
    ### END SOLUTION
end;
In [6]:
funcs = [x -> 1, x -> x, x -> x^2, x -> x^3]
@assert integrate_via_interp.(funcs, 20, 1) ≈ [2, 0, 2/3, 0]
@assert integrate_via_interp(exp, 20, 1) ≈ exp(1) - exp(-1)
@assert integrate_via_interp(cos, 20, 5) ≈ sin(5) - sin(-5)
@assert integrate_via_interp(x -> 1 + 3x^3, 5, 1) ≈ 2.0
@assert abs(integrate_via_interp(f, 20, 5) - 34.070458) ≤ 0.5

[Exercice 2] Une méthode d'ordre élevé pour les équations non-linéaires¶

Le but de cet exercice est de mettre en œuvre une méthode de résolution pour l'équation non-linéaire scalaire suivante : $$ f(x) = 0, \qquad f \colon \mathbb R \to \mathbb R $$ Inspiré par la méthode de Newton-Raphson, qui est basée sur une approximation affine de $f$ à chaque itération, on a l'idée de définir un schéma itératif basé sur des approximations quadratiques de $f$. Plus précisément, on considère la méthode itérative qui consiste à définir, étant donné une approximation $x_k$, l'approximation $x_{k+1}$ de la manière suivante :

  • Considérons le polynôme de Taylor quadratique $p_k$ approximant $f$ autour de $x_k$: $$ p_k(x) = \frac{f''(x_k)}{2} (x - x_k)^2 + f'(x_k) (x - x_k) + f(x_k). $$

  • Si ce polynôme admet deux racines $r_1 \leqslant r_2$ (éventuellement égales), alors on pose $x_{k+1} = r_1$ si $|x_k - r_1| \leqslant |x_k - r_2|$, et $x_{k+1} = r_2$ sinon.

  • Si ce polynôme n'admet pas de racine, on considérera que la méthode ne fonctionne pas. La fonction à implémenter ci-dessous devra retourner nothing.

  1. Implémenter la méthode suggérée et retourner l'ensemble des itérations, y compris le point initial x passé en argument et l'approximation finale, dans un vecteur xs. On supposera que l'argument x est un nombre flottant (Float16, Float32, Float64 ou BigFloat), et on demande de retourner dans le vecteur xs des flottants du même type. Ici f est une fonction, df est sa dérivée, et ddf sa dérivée seconde.
In [7]:
function super_newton(x, f, df, ddf; maxiter=100, ε = 1e-12)
    xs = [x]
    ### BEGIN SOLUTION
    for n in 1:maxiter
        abs(f(x)) < ε && return xs
        a, b, c = ddf(x), df(x), f(x)
        Δ = b^2 - 2*a*c
        Δ < 0 && return nothing
        x1, x2 = x + (-b - √Δ) / a, x + (-b + √Δ) / a
        x = abs(x - x1) < abs(x - x2) ? x1 : x2
        push!(xs, x)
    end
    ### END SOLUTION
    return xs
end;
In [8]:
funs = x -> exp(x) - 2, exp, exp
@assert super_newton(1., funs...) |> typeof === Vector{Float64}
@assert super_newton(1.f0, funs...) |> typeof === Vector{Float32}
@assert super_newton(BigFloat(1.), funs...) |> typeof === Vector{BigFloat}
@assert super_newton(1., funs...)[end] ≈ log(2)
@assert super_newton(3., x -> x^3 - x - 6, x -> 3*x^2 - 1, x -> 6x)[end] ≈ 2
@assert super_newton(1., x -> sin(x) - .5, cos, x -> -sin(x))[end] ≈ π/6
@assert super_newton(1., x -> x^3, x -> 3x^2, x -> 6x) === nothing
  1. Faire un graphique illustrant l'évolution de l'erreur dans le cas où $f(x) = x^3 - 2$. Initialiser l'itération à $x_0 = 1$, et utiliser le format BigFloat avec une grande précision, ainsi qu'une très petite tolérance $ε = 10^{-200}$, afin d'illustrer la convergence aussi clairement que possible. Pour éviter d'encombrer le namespace global, on encapsulera le code du graphique dans une fonction.
In [9]:
# Augmentation de la précision pour le format BigFloat
setprecision(1000, base=10)

# Fonction dont on cherche une racine
f(x) = x^3 - 2.0

function my_plot_error()
    p = plot(title=L"Evolution of the error for $f(x) = x^3 - 2$")
    ### BEGIN SOLUTION
    x = BigFloat(1)
    x∞ = cbrt(BigFloat(2.0))
    df(x) = 3x^2
    ddf(x) = 6x
    xs = super_newton(x, f, df, ddf; ε = 1e-200)
    yplot = abs.(xs .- x∞)
    plot!(yplot, yscale=:log10, label=nothing)
    scatter!(yplot, yscale=:log10, label=nothing)
    ### END SOLUTION
    return p
end
my_plot_error()
Out[9]:
No description has been provided for this image
  1. Supposons qu'une suite $(x_0, x_1, \dotsc)$, dont on ne connaît qu'un nombre fini d'éléments $(x_0, \dotsc, x_K)$, converge à l'ordre $q$ vers une limite $x_{*} \in \mathbb R$, c'est à dire que $$ \lim_{k \to \infty} \frac{|x_{k+1} - x_{\infty}|}{|x_{k} - x_{\infty}|^q} \in \mathbb R. $$ Écrire une fonction estimate_order(xs, x∞) permettant d'estimer l'ordre de convergence $q$.
In [10]:
function estimate_order(xs, x∞)
    ### BEGIN SOLUTION
    logerr = (log∘abs).(xs .- x∞)
    p = fit(logerr[1:end-1], logerr[2:end], 1)
    return p[1]
    ### END SOLUTION
end;
In [11]:
@assert estimate_order([1e-1, 1e-2, 1e-4, 1e-8, 1e-16], 0.) ≈ 2
@assert estimate_order([1e-1, 1e-3, 1e-9, 1e-27], 0.) ≈ 3
@assert estimate_order([1e-1, 1e-2, 1e-3, 1e-4], 0.) ≈ 1

# Estimation de l'ordre de convergence de la méthode `super_newton`
function get_order_super_newton()
    x∞ = cbrt(BigFloat(2.0))
    xs = super_newton(BigFloat(2.), x -> x^3 - 2, x -> 3x^2, x -> 6x; ε = 1e-200)
    q = round(Float64(estimate_order(xs, x∞)), sigdigits=5)
end
println("Order of the proposed method: ", get_order_super_newton())
### BEGIN HIDDEN TESTS
@assert round.(Int, get_order_super_newton()) == 3
### END HIDDEN TESTS

Order of the proposed method: 3.0022
























    

  1. Pour éviter de calculer des dérivées manuellement,
on souhaite utiliser la différentiation automatique.
Écrire une fonction super_newton_auto implémentant la méthode tout en exploitant la différentiation automatique.
Il peut être utile pour ce faire d'utiliser la fonction ForwardDiff.derivative,
et de réutiliser la méthode super_newton définie ci-dessus.
    2. 



















In [12]:






function super_newton_auto(x, f; maxiter=100, ε=1e-12)
    ### BEGIN SOLUTION
    df(x) = ForwardDiff.derivative(f, x)
    ddf(x) = ForwardDiff.derivative(df, x)
    return super_newton(x, f, df, ddf; maxiter=maxiter, ε=ε)
    ### END SOLUTION
end;





















In [13]:






@assert super_newton_auto(1., x -> exp(x) - 2) |> typeof === Vector{Float64}
@assert super_newton_auto(1., x -> exp(x) - 2)[end] ≈ log(2)
@assert super_newton_auto(1., x -> sin(x) - .5)[end] ≈ π/6
@assert super_newton_auto(3., x -> x^3 - x - 6)[end] ≈ 2

























[Exercice 3] Équation de Lotka-Volterra¶






















Dans le cadre d’une étude écologique sur une île isolée,
on aimerait étudier l’interaction entre deux populations animales :


    

  • 🐇 des lapins sauvages, qui se reproduisent rapidement en l’absence de prédateurs,
    • 

  • 🦊 des renards, qui se nourrissent exclusivement de lapins.
    • 



Leur dynamique est modélisée par les équations classiques de Lotka-Volterra :


$$
\left\{
\begin{aligned}
 X'(s) &= \alpha \, X(s) - \beta \, X(s)\,Y(s), \\
 Y'(s) &= - \gamma \, Y(s) + \delta \, X(s)\,Y(s),
\end{aligned}
\right.
$$


où :


    

  • $s$ représente le temps.
    • 

  • $X(s)$ représente la population de lapins à l’instant $ s $,
    • 

  • $Y(s)$ représente la population de renards à l’instant $ s $,
    • 

  • $\alpha$ est le taux de reproduction naturelle des lapins en l’absence de renards,
    • 

  • $\beta$ est la contribution de chaque renard au taux de mortalité des lapins,
    • 

  • $\delta$ est la contribution de chaque lapin au taux de reproduction des renards,
    • 

  • $\gamma$ est le taux de mortalité naturelle des renards en l’absence de nourriture.
    • 



Ce modèle permet d’analyser les cycles prédateur-proie typiques dans les écosystèmes fermés.
Il est simple de montrer que l'unique état d'équilibre non trivial
est donné par $(X_{\rm eq}, Y_{\rm eq}) = (\gamma/\delta, \alpha/\beta)$.
Afin d'éliminer certains paramètres,
on a l'idée d'introduire les variables adimensionnées $X/X_{\rm eq}$ et $Y/Y_{\rm eq}$.
Pour simplifier encore davantage, on prendra comme nouvelle unité temporelle la constante de temps
associée à l'extinction exponentielle de renards en l'absence de lapins,
ce qui revient à poser $x(t) = X(t/γ)/X_{\rm eq}$ et $y(t) = Y(t/γ)/Y_{\rm eq}$.
L'évolution du couple $(x, y)$ est alors donnée par
$$
\left\{
\begin{aligned}
 x'(t) &= r \, x(t) - r \, x(t)\,y(t), \\
 y'(t) &= - y(t) + x(t)\,y(t),
\end{aligned}
\right.
\qquad r := \frac{\alpha}{\gamma}
\tag{Lotka-Volterra}
$$
On supposera dans la première partie de l'exercice que tous les paramètres α,β,γ,δ sont scalaires avec de plus $α = γ$,
de sorte que $r = 1$.
Dans la seconde partie, on adoptera un modèle plus réaliste dans lequel $α$ dépend de $x$.






















    

  1. En posant $\mathbf Z := (x, y)^\top$,
on peut réécrire l'équation de (Lotka-Volterra) sous la forme générique suivante :
$$
\mathbf Z'(t) = \mathbf f\bigl(t, \mathbf Z(s)\bigr).
\tag{Generic-ODE}
$$

Écrire la fonction $\mathbf f$ sous forme d'une fonction Julia f(t, Z).
    2. 



















In [14]:






function f(t, Z)
    ### BEGIN SOLUTION
    x, y = Z
    return [x - x*y, -y + x*y]
    ### END SOLUTION
end;





















In [15]:






@assert f(1, [0, 1]) ≈ [0, -1]
@assert f(1, [1, 0]) ≈ [1, 0]
@assert f(1, [5, 2]) ≈ [5 - 10, -2 + 10]

























    

  1. Écrire une fonction heun3(g, tₙ, Zₙ, Δ) implémentant un pas de temps de taille $\Delta$ de la méthode de Heun d'ordre 3 pour une équation différentielle générique de la forme $\mathbf Z' = \mathbf g(t, \mathbf Z)$.
Cette méthode est basée sur l'itération suivante:
$$
   \mathbf Z_{n+1} = \mathbf Z_n + \frac{\Delta}{4}\left(\mathbf k_1 + 3 \mathbf k_3 \right),
$$
où
\begin{align*}
\mathbf k_1 &= \ \mathbf g(t_n, \mathbf Z_n), \\
\mathbf k_2 &= \ \mathbf g\!\left(t_n + \frac{\Delta}{3}, \mathbf Z_n + \frac{\Delta}{3} \mathbf k_1\right), \\
\mathbf k_3 &= \ \mathbf g\!\left(t_n + \frac{2\Delta}{3}, \mathbf Z_n + \frac{2\Delta}{3} \mathbf k_2\right), \\
\end{align*}
La fonction devra renvoyer $\mathbf Z_{n+1}$.
    2. 



















In [16]:






function heun3(g, tₙ, Zₙ, Δ)
    ### BEGIN SOLUTION
    k₁ = g(tₙ,       Zₙ           )
    k₂ = g(tₙ + Δ/3, Zₙ + Δ/3 * k₁)
    k₃ = g(tₙ + 2Δ/3, Zₙ + 2Δ/3 * k₂)
    return Zₙ + Δ/4 * (k₁ + 3k₃)
    ### END SOLUTION
end;





















In [17]:






@assert heun3((t, Z) -> [1.],0., [0.], 1.) ≈ [1.0]
@assert heun3((t, Z) -> [t], 1., [0.], 1.)  ≈ [3/2]
@assert heun3((t, Z) -> [1.; t; t^2], 0., [2.; 2.; 2.], 1.) ≈ [3., 5/2, 7/3]

























    

  1. Écrire une fonction solve_ode(g, Z₀, Δ, niter) pour résoudre (Generic-ODE)
en effectuant niter itérations de la méthode de Heun d'ordre 3 avec un pas de temps fixe Δ.
Votre fonction devra renvoyer un vecteur de temps ts
(comprenant le temps initial, donc de taille niter + 1),
ainsi qu'un vecteur de vecteurs contenant la solution à ces temps.
    2. 



















In [18]:






function solve_ode(g, Z₀, Δ, niter)
    ### BEGIN SOLUTION
    ts = (0:niter) * Δ
    Z, Zs = Z₀, [Z₀]
    for t in ts[1:end-1]
        Z = heun3(g, t, Z, Δ)
        push!(Zs, Z)
    end
    return ts, Zs
    ### END SOLUTION
end;





















In [19]:






@assert solve_ode(f, [1.0, 0.0], 0.01, 100) |> length == 2
@assert solve_ode(f, [1.0, 0.0], 0.01, 100)[1] |> length == 101
@assert solve_ode(f, [1.0, 0.0], 0.01, 100)[1][101] ≈ 1
@assert solve_ode(f, [1.0, 0.0], 0.01, 100)[1] ≈ 0:.01:1
@assert solve_ode(f, [1.0, 0.0], 0.01, 100)[2] |> x -> sum(z -> z[2], x) == 0.
@assert solve_ode(f, [0.0, 1.0], 0.1, 1000)[2][end][2] ≤ 1e-3

























    

  1. Écrire une fonction my_plot_evolution(g, Z₀, Δ, niter) permettant d'illustrer sur un même graphique l'évolution des populations de lapins et de renards en fonction du temps,
pour une condition initiale Z₀ donnée.
Cette fonction prendra explicitement en argument une function g,
car on l'utilisera plus tard avec un autre modèle d'évolution des populations.
La cellule en dessous vous permettra de tester votre code.
    2. 



















In [20]:






function my_plot_evolution(g, Z₀, Δ, niter)
    x₀, y₀ = Z₀
    p = plot(title="(x₀, y₀) = ($x₀, $y₀)")
    ### BEGIN SOLUTION
    ts, Zs = solve_ode(g, Z₀, Δ, niter)
    xs, ys = first.(Zs), last.(Zs)
    plot!(ts, xs, label="Lapins")
    plot!(ts, ys, label="Renards")
    ### END SOLUTION
    return p
end;





















In [21]:






Δ, niter1, niter2 = .01, 1000, 3000
p1 = my_plot_evolution(f, [1.0, 0.0], Δ, niter1)
p2 = my_plot_evolution(f, [0.0, 1.0], Δ, niter1)
p3 = my_plot_evolution(f, [1.0, 1.0], Δ, niter1)
p4 = my_plot_evolution(f, [3.0, 1.0], Δ, niter2)
p5 = my_plot_evolution(f, [5.0, 2.0], Δ, niter2)
p6 = my_plot_evolution(f, [0.5, 0.5], Δ, niter2)
plot(p1, p2, p3, p4, p5, p6; layout=(2, 3), size=(900, 600))





















Out[21]:



No description has been provided for this image























    

  1. Dans la cellule suivante (qui ne doit en aucun cas être détruite), commenter les figures obtenues au point précédent pour les conditions initiales suivantes
    2. 



    

  • $(0.1, 0.0)$
    • 

  • $(0.0, 1.0)$
    • 

  • $(1.0, 1.0)$
    • 

  • $(3.0, 1.0)$
    • 











































    

  1. Écrire une fonction my_plot_phase(g, all_Z₀, Δ, niter) permettant d'illustrer,
sur un même graphique,
l'évolution des populations de lapins et de renards dans le plan $xy$ (lapins en abscisse et renards en ordonnée)
pour toutes les conditions initiales contenues dans le vecteur all_Z₀, supposé de type Vector{Vector{Float64}}.
La cellule suivante vous permettra de tester votre code.
    2. 



















In [22]:






function my_plot_phase(g, all_Z₀, Δ, niter)
    p = plot(title="Evolution in the phase space")
    ### BEGIN SOLUTION
    for Z₀ in all_Z₀
        x₀, y₀ = Z₀
        _, Zs = solve_ode(g, Z₀, Δ, niter)
        xs, ys = getindex.(Zs, 1), getindex.(Zs, 2)
        plot!(xs, ys, label="(x₀, y₀) = ($x₀, $y₀)")
    end
    ### END SOLUTION
    return p
end;





















In [23]:






all_Z₀ = [[3.0, 1.0],
          [5.0, 2.0],
          [0.5, 0.5]]

my_plot_phase(f, all_Z₀, .01, 3000)





















Out[23]:



No description has been provided for this image























    

  1. Suite à un problème technique, les données sur les conditions initiales de l’étude écologique ont été perdues.
Les populations animales ont néanmoins été observées au temps $t=1$,
où elles étaient données par $(x₁, y₁) = (2, 2)$.
En se basant sur ce point d'observation et sur le modèle d'évolution de (Lotka-Volterra),
estimer les populations au temps $t = 0$.
Il peut être utile pour cela d'utiliser la bibliothèque ForwardDiff.jl,
en particulier la fonction ForwardDiff.jacobian.
    2. 



















In [24]:






function recover_initial_data()
    Z₁ = [2., 2.]
    ### BEGIN SOLUTION
    obj(Z) = solve_ode(f, Z, .01, 100)[2][end] - Z₁
    Z₀ = [5., .1]
    for n in 1:10
        Z₀ -= ForwardDiff.jacobian(obj, Z₀)\obj(Z₀)
    end
    ### END SOLUTION
    return Z₀
end

Z₀ = recover_initial_data()
@show Z₀
my_plot_evolution(f, Z₀, .01, 100)

























Z₀ = [2.2095710044614583, 0.4962424586387766]




















Out[24]:



No description has been provided for this image



















In [25]:






@assert length(Z₀) == 2
### BEGIN HIDDEN TESTS
@assert round.(Z₀, digits=3) ≈ [2.21, 0.496]
### END HIDDEN TESTS

























    

  1. L'évolution de Lotka-Volterra avec $r = 1$ conserve la quantité $V(x, y) = x - \ln(x) + y - \ln(y)$ (pour mémoire le logarithme népérien est donné par log en Julia).
Écrire une fonction calculate_energy_mismatch retournant,
pour la fonction solve_ode avec le pas de temps donné en argument,
la déviation maximale de l'énergie sur l'intervalle de temps $[0, 1]$.
Plus précisément, la fonction devra retourner
$$
E(Δ) := \max_{n \in \mathbb N,  nΔ \in [0, 1]} \Bigl| V\bigl( x_n, y_n \bigr) - V\bigl( x(0), y(0) \bigr) \Bigr|
$$
    2. 



















In [26]:






V(z) = z[1] - log(z[1]) + z[2] - log(z[2])

function calculate_energy_mismatch(Z₀, Δ)
    ### BEGIN SOLUTION
    niter = floor(Int, 1/Δ)
    ts, Zs = solve_ode(f, Z₀, Δ, niter)
    E = maximum(abs, V.(Zs) .- V(Z₀))
    ### END SOLUTION
    return E
end;





















In [27]:






Z₀ = [2.0, 2.0]
@assert calculate_energy_mismatch(Z₀, .1) > 1e-8
@assert calculate_energy_mismatch(Z₀, .1) < 1e-3
@assert calculate_energy_mismatch(Z₀, .01) > 1e-12
@assert calculate_energy_mismatch(Z₀, .01) < 1e-6
@assert calculate_energy_mismatch(Z₀, .001) < 1e8

























    

  1. Il est possible de montrer que $E(Δ) \propto Δ^q$, de manière approximative.
Estimer l'exposant $q$.
Pour cela, approximer par une droite la fonction $\log(Δ) \mapsto \log\bigl(E(Δ)\bigr)$,
et déduire la valeur de $q$ du coefficient directeur de cette droite.
    2. 



















In [28]:






function estimate_q()
    ### BEGIN SOLUTION
    Z₀ = [2.0, 2.0]
    Δs = 1 ./ 2 .^(4:10)
    Es = [calculate_energy_mismatch(Z₀, Δ) for Δ in Δs]
    q = fit(log.(Δs), log.(Es), 1)[1]
    ### END SOLUTION
    return q
end;

println("Order of convergence: ", estimate_q())

























Order of convergence: 3.023790725383018
































In [29]:






@assert estimate_q() > 0
### BEGIN HIDDEN TESTS
@assert round.(Int, estimate_q()) == 3
### END HIDDEN TESTS

























    

  1. On considère à partir de maintenant un modèle plus réaliste,
incluant une saturation du taux de reproduction des lapins
lorsque leur population devient trop importante.
Concrètement, on considère que le taux de reproduction des lapins est désormais donné par $\widehat α(x) = α (1-\frac{x}{4})$,
où $α$ désigne le même paramètre scalaire que précédemment.
Ce choix mène aux équations d'évolution suivantes :
$$
\left\{
\begin{aligned}
 x'(t) &= \left( 1 - \frac{x(t)}{4} \right) x(t) - x(t)\,y(t), \\
 y'(t) &= - y(t) + x(t)\,y(t).
\end{aligned}
\right.
\tag{Modified-Lotka-Volterra}
$$
Avec ce modèle, le nombre de lapins en l'absence de renards ne diverge plus dans la limite $t \to \infty$.
Tout comme les équations de Lotka-Volterra classiques étudiées précédemment,
les équations (Modified-Lotka-Volterra) peuvent s'écrire sous la forme d'une équation différentielle générique $\mathbf Z' = \mathbf g(t, \mathbf Z)$.
Écrire la fonction $\mathbf g$ correspondante sous forme d'une fonction Julia.
    2. 



















In [30]:






function g(t, Z)
    ### BEGIN SOLUTION
    x, y = Z
    return [x*(1-x/4) - x * y,  y * (x - 1)]
    ### END SOLUTION
end;





















In [31]:






@assert g(1, [0, 1]) ≈ [0, -1]
@assert g(0, [4, 1]) ≈ [-4, 3]
@assert g(10, [8, 2]) ≈ [-24.0, 14.0]

























Il est possible de montrer qu'avec ce nouveau modèle,
pour toute condition initiale telle que $x₀ > 0$ et $y₀ > 0$, les populations de lapins et de renards tendent vers l'équilibre $(x_{\infty}, y_{\infty}) := \bigl( 1, \frac{3}{4} \bigr)$ dans la limite $t \to \infty$,
comme illustré ci-dessous.


















In [32]:






Δ, niter = .01, 8_000
p1 = my_plot_evolution(g, [1.0, .1], Δ, niter)
p2 = my_plot_evolution(g, [3.0, .1], Δ, niter)
p3 = my_plot_evolution(g, [0.5, .75], Δ, niter)
plot(p1, p2, p3; layout=(1, 3), size=(900, 300))





















Out[32]:



No description has been provided for this image























    

  1. On suppose maintenant que les populations animales de l'île sont à l'équilibre pour $t < 0$,
donc en particulier $x(0^-) = x_{\infty}$ et $y(0^-) = y_{\infty}$.
Afin de repeupler les îles voisines, des lapins sont extraits de l'île considérée au temps $t = 0$.
    

    Soucieux de préserver l'équilibre écologique de l'île,
les habitants souhaitent déterminer la proportion maximale de lapins pouvant être extraite
tout en s'assurant que la population de renards ne s'écarte jamais (pour aucun $t \geq 0$) de plus de 20% de sa valeur d'équilibre $y_{\infty} = \frac{3}{4}$.
    

    Pour cela, commencer par écrire une fonction max_fox_deviation afin d'estimer la déviation maximale de la population de renards,
en proportion (la fonction devra donc retourner un réel compris entre 0 et 1).
    

    2. 



















In [33]:






function max_fox_deviation(x₀)
    ### BEGIN SOLUTION
    x∞, y∞ = 1, .75
    Z₀ = [x₀, y∞]
    Δ, niter = .01, 10_000
    ts, Zs = solve_ode(g, Z₀, Δ, niter)
    ys = getindex.(Zs, 2)
    return maximum(abs, ys .- y∞) / y∞
    ### END SOLUTION
end;





















In [34]:






@assert max_fox_deviation(0.) ≈ 1.
@assert max_fox_deviation(1.) ≈ 0.
@assert max_fox_deviation(1.e-10) > 1.9
@assert round(max_fox_deviation(0.5), digits=3) ≈ 0.481

























    

  1. Estimer graphiquement la valeur minimale de $x₀$ qui est admissible. On pourra utiliser les fonctions hline et vline (ou leurs homologues hline! et vline!) pour tracer des droites horizontales ou verticales.
    2. 



















In [35]:






### BEGIN SOLUTION
x = LinRange(0, 1, 1000)
plot(x, max_fox_deviation.(x), label=nothing)
hline!([.2], label=nothing)
vline!([.8], label=nothing)
ylims!(0, 2.)
### END SOLUTION





















Out[35]:



No description has been provided for this image























    

  1. (Bonus) Ayant constaté à la question précédente l'existence d'un $x_0$ minimal
satisfaisant la contrainte sur la population de renard,
on se propose ici de déterminer cette valeur de manière très précise.
Pour cela, on suggère de procéder par une méthode de Newton-Raphson s'appuyant sur les nombres duaux.
Écrire une fonction newton_raphson_dual(x, f, maxiter=100; ε = 1e-12) renvoyant une racine de la fonction f en partant d'un point initial x.
    

    

 Indication (cliquer pour afficher)


    On rappelle que l'on peut obtenir simultanément la valeur et la dérivée d'une fonction f en x par la méthode suivante
    

    y = f(ForwardDiff.Dual(x, 1.))
fx = y.value  # renvoie f(x)
dfx = y.partials[1]  # renvoie f'(x)

    

    

    2. 



















In [36]:






function newton_raphson_dual(x, f; maxiter=100, ε = 1e-12)
    ### BEGIN SOLUTION
    for i in 1:maxiter
        y = f(ForwardDiff.Dual(x, 1.))
        x -= y.value/y.partials[1]
        norm(f(x)) < ε && return x
    end
    return nothing
    ### END SOLUTION
end;





















In [37]:






@assert newton_raphson_dual(1, x -> x^2 - 2) ≈ √2
@assert newton_raphson_dual(-1, x -> x^2 - 2) ≈ -√2
@assert newton_raphson_dual(1, x -> x^3 - 2) ≈ cbrt(2)
@assert newton_raphson_dual(2, x -> cos(x) - .5) ≈ acos(.5)

























    

  1. (Bonus) À l'aide des fonctions max_fox_deviation et newton_raphson_dual,
déterminer la valeur minimale de x₀ qui est acceptable,
et vérifier cette solution en traçant un graphe présentant l'évolution des populations pour le x₀ obtenu.
    2. 



















In [38]:






### BEGIN SOLUTION
x₀ = newton_raphson_dual(.8, x -> max_fox_deviation(x) - .2)
y₀ = .75

Δ, niter = .01, 10_000
my_plot_evolution(g, [x₀, y₀], Δ, niter)
hline!([.75*.8, .75*1.2], label="Maximum deviation allowed")
### END SOLUTION





















Out[38]:



No description has been provided for this image



















In [39]:






@assert x₀ > 0.5
### BEGIN HIDDEN TESTS
@assert round.(x₀, digits=3) ≈ 0.789
### END HIDDEN TESTS