Examen de rattrapage

• 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 cellulle ci-dessous importe les bibliothèques utilisées dans ce notebook. Si une ou plusieurs d'entre elles manquent sur votre machine, vous êtes invités à les installer au préalable.

using ForwardDiff
using LaTeXStrings
using LinearAlgebra
using Plots
using Polynomials

[Exercice 1] Quadrature et fonction Gamma

1. Écrire une fonction composite_midpoint(u, a, b, n) permettant d'approximer l'intégrale $$ I := \int_a^b u(x) \, \mathrm d x $$ par la méthode des points milieux composites avec n intervalles de longueur égale. Vérifier la fonction avec quelques cas simples en utilisant @assert.

function composite_midpoint(u, a, b, n)
    ### BEGIN SOLUTION
    Δ = (b - a) / n
    mids = a .+ (0:n-1) * Δ .+ Δ/2
    return Δ * sum(u, mids)
    ### END SOLUTION
end;

# Add additional tests here
@assert composite_midpoint(sin, 0, π, 100_000) ≈ 2

### BEGIN HIDDEN TESTS
@assert composite_midpoint(x -> 1, 0, 1, 100_000) ≈ 1
@assert composite_midpoint(x -> x, 0, 1, 100_000) ≈ 1/2
@assert composite_midpoint(x -> x^2, 0, 1, 100_000) ≈ 1/3
### END HIDDEN TESTS

2. On rappelle que la fonction Gamma est définie (pour $s>0$) par $$ \Gamma(s) \;=\; \int_{0}^{\infty} x^{s-1} e^{-x}\,\mathrm dx. $$ Écrire une fonction Gamma_approx(s; n=1000, L=10) qui retourne une approximation de $\Gamma(s)$ en intégrant sur $[0,L]$ avec la règle des points milieux composites et n intervalles.

function Gamma_approx(s; n=1000, L=10)
    ### BEGIN SOLUTION
    composite_midpoint(x -> x^(s-1) * exp(-x), 0, L, n)
    ### END SOLUTION
end;

En utilisant cette fonction, écrire dans la cellule suivante un code pour afficher des approximations des valeurs suivantes :

• $\Gamma(1) = 1$
• $\Gamma(2) = 1$
• $\Gamma(3) = 2$
• $\Gamma(1/2) = \sqrt{\pi}$.

### BEGIN SOLUTION
@show Gamma_approx(1)
@show Gamma_approx(2)
@show Gamma_approx(3)
@show Gamma_approx(.5)
@show √π;
### END SOLUTION

Gamma_approx(1) = 0.9999504336048892

Gamma_approx(2) = 0.9995047691053153
Gamma_approx(3) = 1.9944612237751818
Gamma_approx(0.5) = 1.7118894543851457
√π = 1.7724538509055159

3. Étudier numériquement l’erreur de cette approximation pour $s=1/2$ en fonction de n, pour une valeur de L égale à 10. Utiliser comme valeur de référence $\sqrt{\pi}$, et tracer l’erreur en fonction de n en échelle log-log.

### BEGIN SOLUTION
I_exact = √π
ns = 2 .^(5:20)
L = 10
I_mid = [Gamma_approx(0.5; n=n, L=L) for n in ns]
p = plot(title="Approximation of Γ(1/2)")
plot!(ns, abs.(I_mid .- I_exact),
      xaxis=:log, yaxis=:log, label=L"Error for $L = 10$", lw=2, xlabel="n")
### END SOLUTION

4. Estimer l’ordre de convergence de la méthode des points milieux composite pour l'intégrande $x^{s-1} e^{-x}$ avec $s=1/2$, sur base du comportement de l’erreur en fonction de $n$ pour $L = 10$.

exponent = 0  # valeur à changer après avoir déterminé l'ordre de convergence

### BEGIN SOLUTION
ns = 2 .^(5:22)
Ls = [2, 5, 10]
I_mid = [Gamma_approx(0.5; n=n, L=10) for n in ns]
log_n = log.(ns)
log_e = log.(abs.(I_mid .- I_exact))
exponent = fit(log_n, log_e, 1).coeffs[2]
### END SOLUTION

println("The error scales approximately as n^($exponent), up to a constant factor")

The error scales approximately as n^(-0.5004376705481719), up to a constant factor

5. (Bonus) Commenter dans la cellule suivante l'ordre de convergence obtenu au point précédent.

[Exercice 2] Méthode d'Euler symplectique et systèmes dynamiques

On considère dans cet exercice des systèmes dynamiques hamiltoniens de la forme $$ \begin{cases} \dot{q} = \partial_p H(q, p), \\ \dot{p} = - \partial_q H(q, p), \end{cases} $$ où $H(q,p)$ est la fonction hamiltonienne du système et où, à un temps donné, $p$ et $q$ sont des quantités scalaires. La méthode d'Euler symplectique est une méthode numérique spécifiquement adaptée à la résolution de tels systèmes. Cette méthode s'écrit $$ \begin{cases} p_{n+1} = p_n - h \, \partial_q H(q_n, p_n), \\ q_{n+1} = q_n + h \, \partial_p H(q_n, p_{n+1}), \end{cases} $$ où $h$ est le pas de temps. On verra que cette méthode conserve mieux l’énergie que les méthodes d'Euler explicite ou implicite.

1. Écrire une fonction symplectic_euler(q0, p0, H, h, N) qui intègre le système hamiltonien sur N pas de temps, en utilisant la bibliothèque ForwardDiff pour calculer les dérivées partielles de H.

   • La dérivée partielle par rapport à x d'une fonction (x, y) -> f(x, y) en un point (x0, y0) peut s'obtenir avec la commande

         ForwardDiff.derivative(x -> f(x, y0), x0)
     

     et la dérivée partielle par rapport à y peut s'obtenir d'une manière similaire.

   • H est une fonction H(q,p) retournant la valeur de la fonction hamiltonienne.

   • La fonction doit retourner un tuple (q,p), ce qui signifie que la dernière ligne de la fonction doit être return q, p, où q et p sont des vecteurs de taille N+1 contenant les valeurs de $q_n$ et $p_n$ pour $n=0, \dotsc,N$.

function symplectic_euler(q0, p0, H, h, N)
    q = zeros(N+1)
    p = zeros(N+1)
    ### BEGIN SOLUTION
    q[1], p[1] = q0, p0
    ∂H_∂q(q, p) = ForwardDiff.derivative(q -> H(q,p), q)
    ∂H_∂p(q, p) = ForwardDiff.derivative(p -> H(q,p), p)
    for n in 1:N
        p[n+1] = p[n] - h * ∂H_∂q(q[n], p[n])
        q[n+1] = q[n] + h * ∂H_∂p(q[n], p[n+1])
    end
    ### END SOLUTION
    return q, p
end;

@assert symplectic_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[1] |> typeof === Vector{Float64}
@assert symplectic_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[2] |> typeof === Vector{Float64}
@assert symplectic_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[1] |> length == 10^5 + 1
@assert symplectic_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[1][end] ≈ -1
@assert symplectic_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[2][end] |> abs ≤ 1e-5

2. Écrire une fonction explicit_euler(q0, p0, H, h, N) qui intègre le système hamiltonien sur N pas de temps, en utilisant la méthode d'Euler explicite classique et la bibliothèque ForwardDiff pour calculer les dérivées partielles de H.

   Pour cela, remarquons que le système peut s'écrire sous la forme standard $\dot X = f(X)$ avec $X = (q, p)$ et $f(X) = (\partial_p H(q,p), -\partial_q H(q,p))$.

function explicit_euler(q0, p0, H, h, N)
    q = zeros(N+1)
    p = zeros(N+1)
    ### BEGIN SOLUTION
    q[1], p[1] = q0, p0
    ∂H_∂q(q,p) = ForwardDiff.derivative(q -> H(q,p), q)
    ∂H_∂p(q,p) = ForwardDiff.derivative(p -> H(q,p), p)
    for n in 1:N
        q[n+1] = q[n] + h * ∂H_∂p(q[n], p[n])
        p[n+1] = p[n] - h * ∂H_∂q(q[n], p[n])
    end
    ### END SOLUTION
    return q, p
end;

@assert explicit_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[1] |> typeof === Vector{Float64}
@assert explicit_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[2] |> typeof === Vector{Float64}
@assert explicit_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[1] |> length == 10^5 + 1
@assert explicit_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^6, 10^6)[1][end] ≤ -1
@assert explicit_euler(1.0, 0.0, (q,p) -> .5 * (p^2 + q^2), π/10^5, 10^5)[2][end] |> abs ≤ 1e-5

3. Considérer la fonction hamiltonienne simple (oscillateur harmonique) $$ H(q,p) = \frac{1}{2} (p^2 + q^2), $$ avec q0=1, p0=0. Pour la dynamique correspondante, comparer l’évolution de l’énergie pour les méthodes d'Euler explicite et symplectique avec h=0.1 et N=100, en traçant $H(q_n, p_n)$ en fonction du temps pour les deux méthodes, sur un même graphique.

### BEGIN SOLUTION
H(q,p) = .5 * (p^2 + q^2)
q0, p0 = 1.0, 0.0
h, N = 0.1, 100

q_symp, p_symp = symplectic_euler(q0, p0, H, h, N)
q_exp, p_exp = explicit_euler(q0, p0, H, h, N)

# Tracé des énergies
t = 0:h:(N*h)
plot(t, H.(q_symp,p_symp), label="Symplectic Euler", lw=2)
plot!(t, H.(q_exp,p_exp), label="Euler explicite", lw=2)
xlabel!("Time"); ylabel!("Energy"); title!("Énergie - Oscillateur harmonique")
### END SOLUTION

4. (Bonus) Étudier l’effet du pas de temps h sur la conservation de l’énergie pour la méthode symplectique. Pour cela, tracer l’erreur en énergie maximale, c'est-à-dire la quantité $\max_n |H(q_n,p_n)-H(q_0,p_0)|$ où le maximum porte sur tous les $n$ tels que le temps correspondant appartienne à $[0, 1]$, en fonction de h avec échelle logarithmique sur les deux axes. Déterminer la valeur de order telle que l'erreur se comporte en $O(h^{order})$ lorsque h est suffisamment petit.

order = 0  # valeur à changer après avoir déterminé l'ordre de convergence

# ### BEGIN SOLUTION
Ns = 2 .^ (10:16)
hs = 1 ./ Ns

function energy_discrepancy(h, N)
    q_symp, p_symp = symplectic_euler(q0, p0, H, h, N)
    return maximum(abs, H.(q_symp, p_symp) .- H(q0,p0))
end

err_energy = energy_discrepancy.(hs, Ns)
order = fit(log.(hs), log.(err_energy), 1)[1]

plot(hs, err_energy, xaxis=:log, yaxis=:log, marker=:o, lw=2,
     xlabel="h", label="Energy discrepancy")
### END SOLUTION

println("The error scales approximately as h^($order) for small h")

The error scales approximately as h^(0.999897418037422) for small h

[Exercise 3] Chebyshev à la rescousse pour les problèmes de stabilité

Dans cet exercice, on s'intéresse à l'équation de la chaleur unidimensionnelle sous condition aux limites de Dirichlet homogène : $$ \tag{EDP} \partial_t u = \partial_x^2 u, \qquad u(0, t) = u(1, t) = 0, \qquad u(x, 0) = v(x), $$ où $v(x) = \sin(\pi x) - \sin(3\pi x) + 2 \sin(4\pi x)$ est une condition initiale donnée. On admet qu'il existe une unique solution classique $u(x, t)$ au problème (EDP), qui est donnée par $$ u(x, t) = \sin(\pi x) \mathrm e^{- \pi^2 t} - \sin(3\pi x) \mathrm e^{- 9 \pi^2 t} + 2 \sin(4\pi x) \mathrm e^{- 16 \pi^2 t}. $$ La condition initiale et la solution sont définies ci-dessous:

# Initial condition
u_init(x) = sinpi(x) - sinpi(3x) + 2sinpi(4*x);

# Exact solution
u_exact(x, t) = sinpi(x) * exp(-π^2*t) -
                sinpi(3x) * exp(-9*π^2*t) +
                2sinpi(4*x) * exp(-16*π^2*t);

Pour construire une méthode d'approximation numérique de la solution, nous introduisons une grille d'abcisses équidistantes : $$ 0 = x_0 < x_1 < \dots < x_{N-1} < x_N = 1, \qquad x_i = \frac{i}{N}. $$ Soient $(u_i)_{i \in \{0, \dots N\}}$ les valeurs de la solution à ces points. Aux limites, il est clair que $u_0(t) = u_N(t) = 0$ pour tout $t$. Aux points intérieurs, une approximation par différence finie centrale de la dérivée seconde donne $$ \partial_x^2 u(x_i, t) \approx \frac{1}{Δx^2} \Bigl( u_{i-1}(t) - 2u_i(t) + u_{i+1}(t) \Bigr), \qquad \Delta x = \frac{1}{N}, \qquad i \in \{1, \dotsc, N-1\}. $$ On en déduit, par l'équation (EDP) et par la condition initiale, que le vecteur $\mathbf u(t) = \bigl(u_1(t), \dotsc, u_{N-1}(t) \bigr)^T$ satisfait approximativement l'équation différentielle suivante: $$ \tag{ODE} \frac{\mathrm d}{\mathrm d t} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-2} \\ u_{N-1} \end{pmatrix} \approx \frac{\mathsf A_N}{Δx^2} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-2} \\ u_{N-1} \end{pmatrix}, \qquad \begin{pmatrix} u_1(0) \\ u_2(0) \\ u_3(0) \\ \vdots \\ u_{N-2}(0) \\ u_{N-1}(0) \end{pmatrix} = \begin{pmatrix} v(x_1) \\ v(x_2) \\ v(x_3) \\ \vdots \\ v(x_{N-2}) \\ v(x_{N-1}) \end{pmatrix}, $$ où $\mathsf A_N \in \mathbb R^{(N-1)\times (N-1)}$ est la matrice suivante, de dimensions (N-1) x (N-1): $$ \begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & 1 & -2 & 1 \\ & & \ddots & \ddots & \ddots & \\ & & & 1 & -2 & 1 \\ & & & & 1 & -2 \\ \end{pmatrix}. $$ L'équation (ODE) n'est satisfaite qu'approximativement par la solution exacte de l'EDP, mais elle constitue une équation différentielle bien posée dont la résolution fournit une approximation de la solution exacte aux points de discrétisation; c'est la méthode des différences finies.

1. Sans utiliser la bibliothèque LinearAlgebra, écrire une fonction dominant_eigenvalue(N; ε=1e-8) permettant de calculer la plus grande valeur propre (en valeur absolue) de la matrice $\mathsf A_N$. Une fonction vous est donnée pour construire $\mathsf A_N$, ainsi qu'une cellule de test. Si votre implémentation est correcte, il devrait apparaître clairement que la valeur propre dominante de $\mathsf A_N$ est bornée en valeur absolue par 4, uniformément en $N$.

A_mat(N) = SymTridiagonal(fill(-2., N-1), fill(1., N-2))

function dominant_eigenvalue(N; ε=1e-8)
    ### BEGIN SOLUTION
    A = A_mat(N)
    N = size(A, 1) + 1
    x = ones(N-1) / √(N-1)
    λ = 0.
    while norm(A*x - λ*x) > ε
        x = A*x
        x /= √(x'x)
        λ = x'A*x
    end
    ### END SOLUTION
    return λ
end;

ns = 10:50
λs = dominant_eigenvalue.(ns)
plot(ns, λs, label=nothing)
scatter!(ns, λs, label="Dominant eigenvalue")

2. On s'intéresse tout au long du reste de cet exercice à une méthode de résolution de l'equation différentielle (ODE) basée sur une itération explicite du type $$ \mathbf u_{n+1} = p\left(\mathsf A_N \frac{Δt}{Δx^2} \right) \mathbf u_{n} \tag{ODE-Scheme} $$ où $p$ est un polynôme. Remarquons que pour $p(x) = 1 + x$, cette méthode coincide avec le schéma d'Euler explicite.

   Écrire une fonction solve_pde(N, Δt, T, p) permettant de résoudre approximativement l'équation (ODE) sur l'intervalle de temps $[0, T]$ par la méthode proposée, avec un pas de temps $\Delta t$ et un pas de discrétisation en espace $\Delta x = 1/N$. Votre fonction devra renvoyer une structure contenant un vecteur de temps ts et un vecteur de vecteurs us contenant la solution à ces temps.

   Hint (click to display)

   • La fonction push!(us, u) permet d'insérer u à la fin de la collection us.

         a = [1., 2., 3.]  # typeof(a) = Vector{Float64}
         vec = [a]  # typeof(vec) = Vector{Vector{Float64}}
         push!(vec, [2., 3., 4.])  # vec = [[1.0, 2.0, 3.0], [2.0, 3.0, 4.0]]
     

   • Pour $N$ très grand, il serait plus efficace de calculer directement les produits matrice-vecteur $p(\mathsf A_N \Delta t / \Delta x^2) \mathbf u_n$ sans jamais calculer de produits matrice-matrice, donc sans jamais former $p(\mathsf A_N \Delta t / \Delta x^2)$ explicitement. Cette approche n'est pas demandée dans cet exercice.

struct SolPoisson
    ts::Vector{Float64}
    us::Vector{Vector{Float64}}
end

# Polynomial `p` for explicit Euler
p_Euler(M) = (I + M)

function solve_pde(N, Δt, T, p)
    ### BEGIN SOLUTION
    ts = Δt * (0:ceil(Int, T/Δt))
    Δx = 1/N
    x = Δx*(1:N-1)
    u = u_init.(x)
    A = A_mat(N)
    us = [u]
    B = p(A*Δt/Δx^2)
    for i in 1:length(ts)-1
        u = B*u
        push!(us, u)
    end
    ### END SOLUTION
    return SolPoisson(ts, us)
end;

N = 20; Δt = 1/(2*N^2); T = .1;
sol = solve_pde(N, Δt, T, p_Euler)

@assert typeof(sol.ts) === Vector{Float64}
@assert typeof(sol.us) === Vector{Vector{Float64}}
@assert length(sol.ts) == length(sol.us)
@assert sol.ts[1:11] ≈ Δt * (0:10)
@assert sol.us[1] ≈ u_init.((1:N-1)/N)

function animate_solution(sol::SolPoisson)
    ts, us = sol.ts, sol.us
    N = length(us[1]) + 1
    Δx = 1/N
    x = Δx*(0:N)
    anim = @animate for (i, t) in enumerate(ts)
        plot(title="t = $(round(t, digits=3))")
        plot!(x -> u_exact(x, t), label="Exact solution")
        plot!(x, [0; us[i]; 0], label="Numerical solution")
        scatter!([x], [0; us[i]; 0], label=nothing)
        xlims!(0, 1)
        ylims!(-3, 4)
    end
    return gif(anim, fps=5, show_msg=false)
end

animate_solution(sol)

4. Pour éviter la divergence du schéma itératif (ODE-Scheme), il est suffisant que l'inégalité $|p(λ)| \leqslant 1$ soit satisfaite pour toute valeur propre $λ$ de la matrice $\mathsf A_N \frac{Δt}{Δx^2}$. En admettant que les valeurs propres de $\mathsf A_N$ sont toutes comprises dans l'intervalle $]-4, 0[$, déterminer une condition suffisante pour garantir la non-divergence de la méthode d'Euler explicite: $$ \text{\color{orange} Écrire la condition de non-divergence dans la cellule suivante.} $$

   Remarque (cliquer pour afficher)

   On parle ici de non-divergence du schéma, plutôt que de stabilité absolue, car on souhaite considérer l'inégalité $|p(λ)| \leqslant 1$ sans inégalité stricte, afin de simplifier la suite de l'exercice. Notons que les deux concepts sont fortement liés, et par abus de langage, on utilisera les termes non-divergence et stabilité de manière synonyme dans le reste de l'exercice.

5. À l'aide de la fonction animate_solution qui est donnée ci-dessus, illustrer la stabilité de la méthode d'Euler explicite quand la condition suffisante de non-divergence est satisfaite, et sa divergence pour certains paramètres tels qu'elle ne l'est pas.

N = 20
T = .1
# Fix parameter Δt such that the Euler scheme is *stable*
### BEGIN SOLUTION
threshold = 1/(2*N^2)
Δt = threshold / 2
sol = solve_pde(N, Δt, T, p_Euler)
animate_solution(sol)
### END SOLUTION

N = 20
T = .1
# Fix parameter Δt such that the Euler scheme is *unstable*
### BEGIN SOLUTION
threshold = 1/(2*N^2)
Δt = threshold * 2
sol = solve_pde(N, Δt, T, p_Euler)
animate_solution(sol)
### END SOLUTION

6. Il est possible de montrer que si $p$ est un polynôme de degré arbitraire tel que $p(0) = 1$ et $p'(0) = 1$, alors le schéma itératif (ODE-Scheme) associé est convergent, dans le sens où la solution exacte de l'équation différentielle est retrouvée dans la limite $Δt \to 0$. Il est donc naturel de chercher, parmi tous les polynômes de degré $m$ tels que $p(0) = 1$ et $p'(0) = 1$, celui qui garantit la plage de non-divergence la plus large possible, c'est à dire celui tel que la condition $$ \forall λ \in [-L, 0], \qquad |p(λ)| \leqslant 1, \tag{Condition} $$ soit satisfaite pour $L > 0$ aussi grand que possible. Il s'avère que ce problème admet une solution explicite; le polynôme optimal n'est autre que le polynòme de Chebyshev de degré $m$, à un changement de variables près : $$ p_{m}(x) = T_m \left(1 + \frac{x}{m^2} \right), \qquad 0 < \alpha \ll 1. $$ Sachant que les polynômes de Chebychev sont donnés par la relation de récurrence $$ \begin{align*} T_0(x) & = 1, \\ T_1(x) & = x, \\ T_{n+1}(x) & = 2 x\,T_n(x) - T_{n-1}(x), \end{align*} $$ écrire une fonction cheb(m, M) permettant de calculer $T_m(\mathsf M)$. La fonction devra être suffisamment générale pour traiter aussi bien le cas d'un argument scalaire que celui d'un argument matriciel. La fonction one est utile pour cela: one(x) retourne 1.0 si x est un scalaire flottant, et la matrice identité de la même taille que x si x est une matrice.

function cheb(m, M)
    if m == 0
        return one(M)
    elseif m == 1
        return M
    else
        return 2*M*cheb(m-1, M) - cheb(m-2, M)
    end
end

cheb (generic function with 1 method)

M = reshape(1.0:9.0, 3, 3) |> Matrix
@assert cheb(0, 3.) ≈ 1.
@assert cheb(0, randn(3, 3)) ≈ I
@assert cheb(1, 3.) ≈ 3.
@assert cheb(1, M) ≈ M
@assert cheb(2, 3.) ≈ 17.
@assert cheb(2, M) ≈ 2M^2 - I
@assert cheb(3, 3.) ≈ 99.
@assert cheb(3, M) ≈ 4M^3 - 3M

7. Ilustrer les polynômes $p_{m}$ pour $m \in \{3, 4, 5\}$ pour $x \in [-51, 1]$.

function p_m(m, x)
    ### BEGIN SOLUTION
    return cheb(m, I + x/m^2)
    ### BEGIN SOLUTION
end

### BEGIN SOLUTION
xplot = LinRange(-51, 1, 500)
plot(xplot, p_m.(3, xplot), xlims=(-51, 1), label=L"p_3")
plot!(xplot, x -> p_m(4, x), xlims=(-51, 1), label=L"p_4")
plot!(xplot, x -> p_m(5, x), xlims=(-51, 1), label=L"p_5")
### END SOLUTION
hline!([-1, 1], label=nothing)
ylims!(-2, 2)

8. Pour $m = 5$, écrire une condition suffisante garantissant la non-divergence du schéma itératif (ODE-Scheme) lorsque $p = p_{m}$. Pour cela, on pourra se rappeler que les polynômes de Chebyshev sont tous bornés par 1 en valeur absolue pour tout $x$ dans l'intervalle $[-1, 1]$. Veiller à vérifier que la condition écrite coincide avec celle du schéma d'Euler explicite pour $m = 1$. $$ \text{\color{orange} Écrire la condition de non-divergence dans la cellule suivante:} $$

9. À l'aide de la fonction animate_solution qui est donnée ci-dessus, illustrer la stabilité de la méthode (ODE-Scheme) avec $p = p_5$ quand la condition suffisante de non-divergence est satisfaite, et son instabilité pour un certain jeu de paramètres tels qu'elle ne l'est pas.

N = 40
T = .1
m = 5
# Fix parameter Δt such that the Euler scheme is *stable*
### BEGIN SOLUTION
p_Cheb = M -> p_m(5, M)
threshold = m^2/(2*N^2)
Δt = threshold / 2
sol = solve_pde(N, Δt, T, p_Cheb)
animate_solution(sol)
### END SOLUTION

N = 40
T = .1
m = 5
# Fix parameter Δt such that the Euler scheme is *unstable*
### BEGIN SOLUTION
p_Cheb = M -> p_m(5, M)
threshold = m^2/(2*N^2)
Δt = threshold * 2
sol = solve_pde(N, Δt, T, p_Cheb)
animate_solution(sol)
### END SOLUTION

10. (Bonus) Le bout de code ci-dessous illustre la région de stabilité absolute de la méthode d'Euler explicite. Modifier ce code afin d'illustrer la région de la stabilité absolue de la méthode ODE-Scheme lorsque $m = 5$ et $p = p_m$.

xlims = (-55, 4)
ylims = (-8, 8)
Plots.plot(size=(900, 300))
x = LinRange(xlims..., 400)
y = LinRange(ylims..., 400)
is_method_stable(z, p) = abs(1 + z) < 1
### BEGIN SOLUTION
m = 5
p_Cheb = M -> cheb(m, I + M/(m^2))
is_method_stable(z, p) = abs(p_Cheb(z)) < 1
### END SOLUTION
stable = is_method_stable.(x' .+ im*y, 4)
Plots.contourf!(x, y, stable, c=:Pastel1_3,
                aspect_ratio=:equal, levels=[0, 1, 2],
                xlabel=L"\Re(\Delta \lambda)",
                ylabel=L"\Im(\Delta \lambda)",
                colorbar=:none, xlims=xlims, ylims=ylims)
Plots.contour!(x, y, stable, levels=5, c=:red)
Plots.vline!([0], color=:gray, label=nothing)
Plots.hline!([0], color=:gray, label=nothing)

11. (Bonus) Lorsqu'elle est implémentée efficacement, la méthode (ODE-Scheme) présente un coût computationnel proportionnel au degré du polynôme $p$. Expliquer dans la cellule suivante l'intérêt d'utiliser les polynômes de Chebyshev.