Semaine 2 : Interpolation et intégration

using LaTeXStrings
using LinearAlgebra
using Plots
using Polynomials

Interpolation et approximation

[Exercice 1] Interpolation de Tchebychev et phénomène de Runge

L'objectif de cet exercice est d'illustrer l'influence du choix des nœuds d'interpolation sur l'erreur d'interpolation.

1. La fonction fit du package Polynomials.jl peut être utilisée de la façon suivante pour calculer, étant donnés des tableaux x et y de même taille, le polynôme interpolateur associé :

       p = fit(x, y)
   

   En utilisant cette fonction, écrire une fonction

       get_interpolations(f, d)
   

   qui interpole la fonction f par un polynôme de degré d. La fonction doit renvoyer un tuple de structures Polynomial, correspondant aux nœuds équidistants (extrémités incluses) et aux nœuds de Tchebychev sur l'intervalle $[-1, 1]$.

   Indice (cliquer pour afficher)

   Pour calculer rapidement les nœuds de Tchebychev, on peut exploiter la macro @. (comme toujours, il est conseillé de se référer à la documentation d'une commande en tapant ? puis la commande dans la console). Cette commande évite d'utiliser des . après chaque fonction ou avant chaque opérateur.

   x = @. -cos(π*((0:n-1)+1/2)/n)
   

function get_interpolations(f, d)
    ### BEGIN SOLUTION
    n = d + 1
    x_equi = LinRange(-1, 1, n)
    x_cheb = @. -cos(π*((0:n-1)+1/2)/n)
    p_equi = Polynomials.fit(x_equi, f.(x_equi))
    p_cheb = Polynomials.fit(x_cheb, f.(x_cheb))
    ### END SOLUTION
    return p_equi, p_cheb
end

get_interpolations (generic function with 1 method)

p_test = Polynomial([1., 2., 3.])
@assert get_interpolations(cos, 5) |> length == 2
@assert get_interpolations(sin∘exp, 5)[1].coeffs |> length == 6
@assert get_interpolations(sin∘exp, 5)[2].coeffs |> length == 6
@assert get_interpolations(p_test, 2)[1] ≈ p_test
@assert get_interpolations(p_test, 2)[2] ≈ p_test
@assert get_interpolations(cos, 4)[1](0) ≈ 1
@assert get_interpolations(cos, 4)[2](0) ≈ 1

2. Fixons $d = 20$ et prenons $f$ égale à la fonction suivante $$ f(x) = \tanh\left(\frac{x+1/2}{\varepsilon}\right) + \tanh\left(\frac{x}{\varepsilon}\right) + \tanh\left(\frac{x-1/2}{\varepsilon}\right), \qquad \varepsilon = .02. $$ En utilisant votre fonction get_interpolations, calculer les polynômes interpolateurs dans ce cas, et afficher l'erreur $L^{\infty}$ correspondant aux nœuds équidistants et aux nœuds de Tchebychev.

   Indice (cliquer pour afficher)

   • Pour limiter les erreurs d'arrondi numérique, il est préférable que la fonction renvoie un type BigFloat, autrement dit

         f(x) = BigFloat(tanh((x+1/2)/ε) + tanh(x/ε) + tanh((x-1/2)/ε))
     

   • Pour calculer la norme infinie d'une fonction afin d'évaluer la précision de l'interpolation, on pourra exploiter la fonction norm(...,Inf) de la bibliothèque LinearAlgebra avec un échantillonnage suffisamment fin des valeurs de la fonction, ou exploiter la fonction maximum :

         maximum(abs, [1, -3, 2])  # = 3
     

     Noter que la conversion d'un nombre y de type BigFloat en Float64 se fait par convert(Float64, y) ou plus simplement ici Float64(y).

d, ε = 20, .02
f(x) = BigFloat(tanh((x+1/2)/ε) + tanh(x/ε) + tanh((x-1/2)/ε))

# Calculate L^∞ errors below
### BEGIN SOLUTION
X = LinRange(-1, 1, 500)
p_equi, p_cheb = get_interpolations(f, d)
round_error(x) = round(x, sigdigits=3)
error_inf_equi = maximum(round_error∘abs, f.(X) - p_equi.(X))
error_inf_cheb = maximum(round_error∘abs, f.(X) - p_cheb.(X))
### END SOLUTION

println("L^∞ error with equidistant nodes: ", error_inf_equi)
println("L^∞ error with Chebyshev nodes: ", error_inf_cheb)

L^∞ error with equidistant nodes: 172.0

L^∞ error with Chebyshev nodes: 0.7319999999999999999999999999999999999999999999999999999999999999999999999999987

@assert (Float64∘round)(error_inf_equi, sigdigits=1) == 200
@assert (Float64∘round)(error_inf_cheb, sigdigits=1) == 0.7

3. Tracer les polynômes interpolateurs superposés à la fonction f.

   Indice (cliquer pour afficher)

   Il peut être utile pour comparer les deux interpolations de limiter les valeurs minimale et maximale sur l'axe y à l'aide de l'option ylims = (ymin,ymax) dans une fonction de tracé plot, ou son équivalent se terminant par !. On rappelle que, par convention en Julia (et non par obligation), une fonction dont le nom se termine par ! modifie ses arguments. Dans le cas d'un graphe, la première commande initiant le graphe ne doit pas comporter de ! (plot) tandis que les suivantes incrémentant le même graphe se terminent par ! (plot!, scatter!, ...). Toute omission du ! est considéré comme un redémarrage à zéro du tracé.

X = LinRange(-1, 1, 500)
plot(X, f.(X), linewidth=4, label="f")

### BEGIN SOLUTION
plot!(X, p_equi.(X), linewidth=3, color=:green, label="Equidistant interpolation")
plot!(X, p_cheb.(X), linewidth=3, color=:red, label="Chebyshev interpolation")
plot!(xlims = (-1, 1), ylims = (-3.5, 3.5))
### END SOLUTION

Intégration numérique

[Exercice 2] Implémentation d'un intégrateur composite

La règle d'intégration de Milne s'écrit $$ \int_{-1}^{1} u(x) \, dx \approx \frac{2}{3} \left( 2 u\left(-\frac{1}{2}\right) - u(0) + 2 u\left(\frac{1}{2}\right) \right) $$

1. Écrire une fonction composite_milne(u, a, b, N), qui renvoie une approximation de l'intégrale $$ \int_{a}^{b} u(x) \, dx $$ obtenue en partitionnant l'intervalle d'intégration $[a, b]$ en $N$ cellules de taille égale, et en appliquant la règle de Milne dans chaque cellule.

function composite_milne(u, a, b, N)
    ### BEGIN SOLUTION
    Δ = (b - a) / N
    x₁ = a .+ Δ/4 .+ Δ*(0:N-1)
    x₂ = a .+ Δ/2 .+ Δ*(0:N-1)
    x₃ = a .+ 3Δ/4 .+ Δ*(0:N-1)
    2Δ/3 * u.(x₁) - Δ/3 * u.(x₂) + 2Δ/3 * u.(x₃) |> sum
    ### END SOLUTION
end

composite_milne (generic function with 1 method)

@assert (abs∘composite_milne)(x -> x, -1, 1, 10) < 1e-13
@assert composite_milne(x -> x, 1, 2, 10) ≈ 3/2
@assert composite_milne(x -> x^2, -1, 1, 1) ≈ 2/3
@assert composite_milne(x -> x^4, -1, 1, 1) ≈ 2/12

2. Prendre $u(x) = \cos(x)$, $a = -1$ et $b = 1$. Tracer à l'aide de scatter l'évolution de l'erreur, en valeur absolue, pour les valeurs de $N$ données, en utilisant une échelle logarithmique pour les deux axes.

u(x) = cos(x)
a, b = -1 , 1

# Number of intervals
Ns = (round∘^).(10, LinRange(0, 3, 20))

# Exact value of the integral
I_exact = 2sin(1)

### BEGIN SOLUTION
Is = composite_milne.(u, a, b, Ns)
errors = abs.(Is .- I_exact)
scatter(Ns, errors, label="Integration error")
### END SOLUTION

# Set log scale for both axes
plot!(xscale=:log10, yscale=:log10)

3. Estimer l'ordre de convergence par rapport à $N$, c'est-à-dire trouver $\gamma$ tel que $$ \lvert \widehat{I}_{N} - I \rvert \propto \beta N^{-\gamma}, $$ où $I$ désigne la valeur exacte de l'intégrale et $\widehat{I}_{N}$ désigne son approximation. Pour trouver $\beta$ et $\gamma$, utiliser la fonction Polynomials.fit pour trouver une approximation linéaire de la forme $$ \log \lvert \widehat{I}_{N} - I \rvert \approx \log (\beta) - \gamma \log(N). $$ Si votre calcul est correct, la fonction N -> β*N^(-γ) devrait donner une bonne approximation de l'erreur d'intégration.

# Calculate β and γ
### BEGIN SOLUTION
p = fit(log.(Ns), log.(errors), 1)
β = round(exp(p[0]), sigdigits=3)
γ = -round(p[1], sigdigits=3)
### END SOLUTION
plot!(N -> β*N^(-γ), label=L"%$β \times N^{%$γ}")

@assert round(β, sigdigits=1) ≤ .1
@assert round(β, sigdigits=1) ≥ 1e-3
@assert round(γ, sigdigits=1) == 4

[Exercice 3] Intégration gaussienne

Notre objectif dans cet exercice est d'écrire un programme pour calculer des intégrales de la forme $$ I[f] := \int_0^{\infty} f(x) \mathrm{e}^{-x} \, \mathrm{d} x $$ À cet effet, nous utiliserons les polynômes de Laguerre, qui sont des polynômes orthogonaux pour le produit scalaire suivant : $$ \langle f, g \rangle := \int_0^{\infty} f(x) g(x) \mathrm{e}^{-x} \, \mathrm{d} x $$ Ces polynômes peuvent être construits par l'algorithme de Gram-Schmidt.

1. En utilisant le fait que les polynômes de Laguerre satisfont la relation de récurrence $$ L_{k + 1}(x) = \frac{(2k + 1 - x)L_k(x) - k L_{k - 1}(x)}{k + 1}, \qquad L_0(x) = 1, \qquad L_1(x) = 1-x, $$ écrire une fonction laguerre(n) qui renvoie le polynôme de Laguerre de degré $n$.

function laguerre(n)
    if n == 0
        return Polynomial([1])
    elseif n == 1
        return Polynomial([1, -1])
    else
        k = n-1
        x = Polynomial([0, 1])
        return ((2k + 1 - x) * laguerre(k) - k*laguerre(k-1))/(k+1)
    end
end

laguerre (generic function with 1 method)

2. Écrire une fonction get_nodes_and_weights(n) qui renvoie les nœuds et les poids de la quadrature de Gauss-Laguerre à $n$ nœuds.

   Indice (cliquer pour afficher)

   • Rappeler que les nœuds de la quadrature sont les racines du polynôme de Laguerre de degré $n$. Pour les trouver, utiliser la fonction roots du package Polynomials.

         p = Polynomial([1, 0, -1])
         r = roots(p)  # r = [-1.0, 1.0]
     

   • Une fois les nœuds de la quadrature trouvés, les poids peuvent être obtenus à partir de la relation $$ \int_0^{\infty} q(x) \, \mathrm{e}^{-x} \, \mathrm{d} x = \sum_{i=1}^n w_i q(x_i), $$ qui doit être vraie pour tout polynôme $q$ de degré au plus $2n - 1$. En prenant $q = \ell_i$ comme le polynôme de Lagrange associé au nœud $i$, on obtient immédiatement que $$ w_i = \int_0^{\infty} \ell_i(x) \, \mathrm{e}^{-x} \, \mathrm{d} x, \qquad \ell_i = \prod_{\substack{j=1 \\ j \neq i}}^n \frac{x - x_j}{x_i - x_j}. $$

   • Pour construire les polynômes de Lagrange $\ell_i$, il peut être utile d'utiliser la fonction fromroots du package Polynomials.

         r = [-1.0, 1.0]
         p = fromroots(r)  # p = Polynomial(-1.0 + 1.0*x^2)
     

     Rappeler aussi que, pour un vecteur x, l'expression x[1:end .!= 5] renvoie le vecteur obtenu en supprimant le cinquième élément de x.

   • Pour calculer l'intégrale des polynômes de Lagrange contre le poids exponentiel, rappeler que $$ \int_0^{\infty} x^n \mathrm{e}^{-x} \, \mathrm{d}x = n! $$

function get_nodes_and_weights(n)
    ### BEGIN SOLUTION
    nodes = roots(laguerre(n))
    weights = zero(nodes)
    for i in 1:n
        ℓ = fromroots(nodes[1:end .!= i])
        ℓ /= ℓ(nodes[i])
        weights[i] = factorial.(0:n-1)'ℓ.coeffs
    end
    return nodes, weights
    ### END SOLUTION
end

get_nodes_and_weights (generic function with 1 method)

@assert get_nodes_and_weights(5) |> length == 2
@assert get_nodes_and_weights(5)[1] |> length == 5
@assert get_nodes_and_weights(5)[2] |> length == 5
@assert get_nodes_and_weights(1)[1] ≈ [1.0]
@assert get_nodes_and_weights(1)[2] ≈ [1.0]
@assert get_nodes_and_weights(3)[1] .|> laguerre(3) |> abs∘sum < 1e-10
@assert get_nodes_and_weights(5)[1] .|> laguerre(5) |> abs∘sum < 1e-10
@assert get_nodes_and_weights(5)[2] |> sum ≈ 1

3. Écrire une fonction integrate_laguerre(f, n), qui renvoie une approximation de $I[f]$ obtenue par intégration de Gauss-Laguerre avec $n$ nœuds.

function integrate_laguerre(f, n)
    ### BEGIN SOLUTION
    nodes, weights = get_nodes_and_weights(n)
    return f.(nodes)'weights
    ### END SOLUTION
end

integrate_laguerre (generic function with 1 method)

@assert integrate_laguerre(x -> x, 5) ≈ 1
@assert integrate_laguerre(x -> x^2, 5) ≈ 2
@assert integrate_laguerre(x -> x^3, 5) ≈ 6
@assert integrate_laguerre(x -> exp(-x), 15) ≈ 1/2
@assert integrate_laguerre(x -> exp(-2x), 15) ≈ 1/3

4. Pour $n = 5$, on calcule numériquement que le degré d'exactitude est égal à 9.

n = 5
for i in 1:9
    correct = integrate_laguerre(x -> x^i, n) ≈ factorial(i)
    println("f = x^$i, Rule exact? ", correct)
    @assert correct
end

f = x^1, Rule exact? true

f = x^2, Rule exact? true
f = x^3, Rule exact? true
f = x^4, Rule exact? true
f = x^5, Rule exact? true
f = x^6, Rule exact? true
f = x^7, Rule exact? true
f = x^8, Rule exact? true
f = x^9, Rule exact? true

5. Poser $f(x) = \sin(x)$, et tracer l'erreur d'intégration en fonction de $n$, en utilisant des échelles appropriées pour les axes x et y.

### BEGIN SOLUTION
ns = 1:20
f(x) = sin(x)
I_exact = 1/2
Ih = integrate_laguerre.(f, ns)
plot(ns, abs.(Ih .- I_exact), yscale=:log10, xlabel=L"n", ylabel="Error")
scatter!(ns, abs.(Ih .- I_exact))
### END SOLUTION

6. L'objectif de cette question est d'étudier une autre manière de calculer les nœuds et les poids de la quadrature de Gauss-Laguerre.

   • Pour les nœuds, remarquer que la récurrence donnée ci-dessus pour les polynômes de Laguerre implique que $$ \begin{pmatrix} x L_0(x) \\ x L_1(x) \\ x L_2(x) \\ \vdots \\ x L_{n-2}(x) \\ x L_{n-1}(x) \end{pmatrix} = \underbrace{ \begin{pmatrix} \beta_0 & \alpha_1 \\ \alpha_1 & \beta_1 & \alpha_2 \\ & \alpha_2 & \beta_2 & \alpha_3 \\ & & \ddots & \ddots & \ddots \\ & & & \alpha_{n-2} & \beta_{n-2} & \alpha_{n-1} \\ & & & & \alpha_{n-1} & \beta_{n-1} \end{pmatrix} }_{\mathsf{A}} \underbrace{ \begin{pmatrix} L_0(x) \\ L_1(x) \\ L_2(x) \\ \vdots \\ L_{n-2}(x) \\ L_{n-1}(x) \end{pmatrix} }_{\mathbf{\boldsymbol{\Phi}}(x)} + \begin{pmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \\ \alpha_{n} L_{n}(x) \end{pmatrix} $$ avec $\alpha_k = k$ et $\beta_k = 2k + 1$. Ainsi, les racines de $L_n$ sont données par les valeurs propres de $\mathsf{A}$. Cela ouvre la voie à une autre approche pour calculer les nœuds d'intégration, en résolvant un problème aux valeurs propres.

   • Une fois calculés les nœuds $x_1, \dotsc, x_n$ de la quadrature de Gauss-Laguerre à $n$ nœuds, c'est-à-dire les racines de $L_n$, les poids peuvent être calculés par l'approche suivante. Notons d'abord que le polynôme de Lagrange associé à $x_i$ peut s'écrire $$ \ell_{i}(x) = \frac{1}{L_n'(x_i)} \frac{L_n(x)}{x - x_i}. $$ Fixons $i \in \{1,\dotsc,n\}$ et posons $q(x) = L_n'(x) \ell_i(x)$. Alors, puisque la quadrature de Gauss-Laguerre à $n$ nœuds doit être exacte pour tout polynôme de degré jusqu'à $2n-1$, on a $$ \int_0^{\infty} q(x) \, \mathrm{e}^{-x} \, \mathrm{d} x = \sum_{i=1}^n w_i q(x_i). $$ Par définition de $q$, le membre de droite vaut $w_i L_n'(x_i)$. D'autre part, par intégration par parties, le membre de gauche vaut $$ \begin{aligned} \int_0^{\infty} q(x) \, \mathrm{e}^{-x} \, \mathrm{d} x &= \int_0^{\infty} L_n'(x) \ell_i(x) \, \mathrm{e}^{-x} \, \mathrm{d} x \\ &= (L_n \ell_i) \Big\vert^{\infty}_0 + \int_0^{\infty} \bigl( L_n(x) \ell_i(x) - L_n(x) \ell_i'(x) \bigr) \mathrm{e}^{-x} \, \mathrm{d} x. \end{aligned} $$ Puisque $L_n(x)$ est orthogonal au sens $\mathrm{L}^2(\mathrm{e}^{-x})$ à tout polynôme de degré $n-1$ ou moins, l'intégrale du membre de droite est nulle. On a donc montré que $$ w_i L_n'(x_i) = - L_n(0) \ell_i(0) = \frac{L_n(0)^2}{L_n'(x_i) x_i}, $$ où l'on a utilisé l'expression de $\ell_i$ dans la dernière équation. En réarrangeant l'équation, on déduit que $$ w_i = \frac{L_n(0)^2}{L_n'(x_i)^2 x_i} $$ Cette relation permet de calculer rapidement les poids une fois les nœuds calculés. Notons également que cette relation reste vraie si l'on multiplie $L_n$ par un facteur constant. Autrement dit, on a $$ w_i = \frac{p_n(0)^2}{p_n'(x_i)^2 x_i}, \qquad p_n := (x - x_1) \dotsc (x - x_n). $$

   Utiliser ces approches pour calculer les nœuds et les poids de la quadrature de Gauss-Laguerre.

   Indice (cliquer pour afficher)

   • Utiliser LinearAlgebra.SymTridiagonal pour construire la matrice A.

         S = SymTridiagonal([1, 2, 3], [4, 5])  # S = [1 4 0; 4 2 5; 0 5 3]
     

   • Les valeurs propres et vecteurs propres d'une matrice peuvent être calculées par LinearAlgebra.eigen.

         E = eigen(S)
         λs = E.values  # eigenvalues from the smallest to the largest
         vs = E.vectors  # eigenvectors in columns
         λi = E.values[i]  # i-th eigenvalue
         vi = E.vectors[:, i]  # i-th eigenvector
     

     On notera d'une part que les valeurs propres sont triées par ordre croissant et sont renvoyées sous forme du vecteur E.values si E désigne le résultat de l'opérateur LinearAlgebra.eigen, et d'autre part que les vecteurs propres sont normalisés et sont les colonnes de la matrice renvoyée par E.vectors, si bien que le $i^{\textrm{ème}}$ vecteur propre est donné par E.vectors[:, i].

   • Utiliser Polynomials.fromroots et Polynomials.derivative pour construire et dériver le polynôme $p_n$.

function get_nodes_and_weights_bis(n)
    ### BEGIN SOLUTION
    αs = -Float64.(1:n-1)
    βs = Float64.(2(0:n-1) .+ 1)
    A = SymTridiagonal(βs, αs)
    nodes = eigen(A).values
    p = fromroots(nodes)
    dp = derivative(p)
    weights = (p(0)./dp.(nodes)).^2 ./ nodes
    return nodes, weights
    ### END SOLUTION
end

@show get_nodes_and_weights(5)[1]
@show get_nodes_and_weights_bis(5)[1]

(get_nodes_and_weights(5))[1] = [0.2635603197181409, 1.413403059106518, 3.5964257710407157, 7.085810005858862, 12.64080084427574]
(get_nodes_and_weights_bis(5))[1] = 

[0.2635603197181408, 1.413403059106515, 3.596425771040715, 7.08581000585883, 12.640800844275773]

5-element Vector{Float64}:
  0.2635603197181408
  1.413403059106515
  3.596425771040715
  7.08581000585883
 12.640800844275773

@assert get_nodes_and_weights_bis(5) |> length == 2
@assert get_nodes_and_weights_bis(5)[1] |> length == 5
@assert get_nodes_and_weights_bis(5)[2] |> length == 5
@assert get_nodes_and_weights_bis(1)[1] ≈ [1.0]
@assert get_nodes_and_weights_bis(1)[2] ≈ [1.0]
@assert get_nodes_and_weights_bis(3)[1] .|> laguerre(3) |> abs∘sum < 1e-10
@assert get_nodes_and_weights_bis(5)[1] .|> laguerre(5) |> abs∘sum < 1e-10
@assert get_nodes_and_weights_bis(5)[2] |> sum ≈ 1

Exercice supplémentaire (optionnel)

[Exercice 4] Résolution de l'équation d'Euler-Bernoulli par interpolation

L'objectif de cet exercice est d'explorer une application pratique de l'interpolation polynomiale. Plus précisément, nous allons implémenter une méthode numérique pour résoudre approximativement l'équation de la poutre d'Euler-Bernoulli avec des conditions aux limites de Dirichlet homogènes :

$$ u\in C^4([0,1]),\quad\left\{\begin{aligned} u''''(x) &= \varphi(x) \qquad \forall\, x\in(0,1),\\ u(0) &= u'(0) = u'(1) = u(1) = 0, \end{aligned}\right. $$ où $\varphi(x) = (2\pi)^4\cos(2\pi x)$ est un chargement transversal appliqué à la poutre. Pour résoudre l'équation numériquement, on approche le membre de droite $\varphi$ par son polynôme interpolateur $\widehat{\varphi}$, puis on résout l'équation exactement avec le membre de droite $\widehat{\varphi}$ à la place de $\varphi$.

1. Commençons par écrire une fonction fit_values_and_slopes(u₀, up₀, u₁, up₁) qui renvoie l'unique polynôme $p$ de degré 3 tel que $$ p(0) = u_0, \qquad p'(0) = up_0, \qquad p(1) = u_1, \qquad p'(1) = up_1. $$

function fit_values_and_slopes(u₀, up₀, u₁, up₁)
    # We look for polynomials p(x) = a₀ + a₁ x + a₂ x² + a₃ x³
    A = [1 0 0 0; 0 1 0 0; 1 1 1 1; 0 1 2 3]
    α = A\[u₀; up₀; u₁; up₁]
    return Polynomial(α)
end

# Test our code
p = fit_values_and_slopes(-1, -1, 1, 1)
plot(p, xlims=(0, 1))

2. Écrire une fonction approx(n) implémentant l'approche décrite ci-dessus pour résoudre l'EDP. La fonction doit renvoyer une approximation polynomiale de la solution basée sur une interpolation de degré $n$ du membre de droite aux points équidistants entre 0 et 1, bornes incluses.

   Indice (cliquer pour afficher)

   • On peut utiliser la fonction Polynomials.fit pour obtenir le polynôme interpolateur :

         p = fit(x, y)
     

     où x sont les nœuds d'interpolation, et y sont les valeurs de la fonction à interpoler.

   • Pour calculer la solution exacte avec un membre de droite polynomial, remarquer que toutes les solutions sont polynomiales, et sans conditions aux limites, la solution est unique modulo un polynôme cubique.

   • On peut utiliser la fonction integrate de la bibliothèque Polynomials.jl, qui calcule une primitive d'un polynôme :

         P = integrate(p)
     

   • Utiliser le format BigFloat pour limiter les erreurs d'arrondi.

         X = LinRange{BigFloat}(0, 1, n + 1)
     

# Right-hand side
φ(x) = (2π)^4 * cospi(2*x)

# Exact solution (for comparison purposes)
u(x) = cospi(2*x) - 1

function approx(n)
    X = LinRange{BigFloat}(0, 1, n + 1)
    ### BEGIN SOLUTION
    Y = φ.(X)
    p = fit(X, Y)
    uh = integrate(integrate(integrate(integrate(p))))
    ∂uh = derivative(uh)
    uh -= fit_values_and_slopes(uh(0), ∂uh(0), uh(1), ∂uh(1))
    return uh
    ### END SOLUTION
end

plot(approx(3), xlims=(0, 1), label="Exact solution")
plot!(u, xlims=(0, 1), label="Approximate solution")

3. Écrire une fonction estimate_error(n) qui approxime l'erreur, en norme $L^{\infty}$, entre la solution exacte et la solution approchée. On rappelle que la solution exacte est donnée par $$ u(x) = \cos(2\pi x) - 1. $$

function estimate_error(n)
    ### BEGIN SOLUTION
    un = approx(n)
    x_fine = LinRange(0, 1, 1000)
    un_fine, u_fine = un.(x_fine), u.(x_fine)
    return maximum(abs.(u_fine - un_fine))
    ### END SOLUTION
end

estimate_error (generic function with 1 method)

4. Tracer l'erreur pour $n$ dans la plage $\{ 5, \dotsc, 50 \}$. Utiliser une échelle logarithmique pour l'axe $y$.

# ### BEGIN SOLUTION
ns = 5:50
errors = estimate_error.(ns)
plot(ns, errors, marker = :circle, label=L"$L^{\infty}$ Error")
plot!(yaxis=:log, lw=2)
# ### END SOLUTION