# Une application de la méthode de Monte-Carlo à la résolution numérique des équations de Fredholm

M. Gatesoupe, J. Guilloud

<figure><img src="image_1.png" /></figure>

## To cite this version:

M. Gatesoupe, J. Guilloud. Une application de la méthode de Monte-Carlo à la résolution numérique des équations de Fredholm. Annales de l'ISUP, 1961, X (4), pp.243-265. hal-04092453

HAL Id: hal-04092453

https://hal.science/hal-04092453

Submitted on 9 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# UNE APPLICATION DE LA MÉTHODE DE MONTE-CARLO
A LA RÉSOLUTION
NUMÉRIOUE DES ÉQUATIONS DE FREDHOLM

M. GATESOUPE et J. GUILLOUD

## I - INTRODUCTION -

L'idée directrice de la méthode de Monte-Carlo est d'associer à un calcul numérique un jeu de hasard dont la règle, adaptée au problème à résoudre, conduit pour chaque partie à la détermination de la valeur d'une variable aléatoire, l'espérance mathématique de ces valeurs étant la quantité cherchée. La moyenne obtenue pour un grand nombre N de parties est donc une approximation de cette quantité. Cette méthode semble présenter un réel intérêt lorsqu'on ne dispose pas de méthode simple de calcul numérique classique. Elle est devenue pratiquement utilisable grâce à la mise en service de machines électroniques à grande rapidité de calcul.

La méthode de Monte-Carlo a été appliquée, de diverses façon, à la résolution des équations linéaires et à l'inversion des opérations linéaires. Nous allons montrer ici comment l'une des méthodes classiques d'inversion des matrices s'adapte à la résolution des équations intégrales de Fredholm, c'est-à-dire à l'inversion d'opérateurs intégraux dans certains espaces fonctionnels.

Soit donc l'équation générale de Fredholm

$$f(x) = g(x) + \int_{a}^{b} H(x,y) f(y) dy \quad (1)$$

dans laquelle $g(x)$ est une fonction donnée sur l'intervalle fermé $[a,b]$, $H(x,y)$ est une fonction donnée de deux variables continue sur le carré compact $[a,b] \times [a,b]$, $f(x)$ est une fonction inconnue continue sur l'intervalle fermé $[a,b]$.

Pour simplifier les notations nous supposerons $a = 0$, $b = 1$, cas auquel on peut toujours se ramener par changement de variable. Nous nous proposons de résoudre numériquement l'équation : pour
x donné $0 \le x \le 1$, nous voulons trouver une valeur approchée de $f(x)$. Soit $\mathcal{C} ([0,1])$ l'espace vectoriel des fonctions numériques continues sur $[0,1]$, c'est un espace de Banach pour la norme :

$$|f| = \sup_{x \in [0,1]} |f(x)|$$

Rappelons que l'algèbre $\mathbb{A}$ des opérateurs linéaires continus $L$ sur $\mathcal{C} ([0,1])$, est aussi un espace de Banach lorsqu'on le munit de la norme :

$$||L|| = \sup_{|f| \le 1} |Lf|;$$

la topologie ainsi définie étant d'ailleurs compatible avec la structure d'algèbre de $\mathbb{A}$.

Soit $I$ l'opérateur identique, élément unité de l'algèbre $\mathbb{A}$. On sait alors que si $||L|| < 1$, l'opérateur $(I-L)$ est inversible de façon unique dans $\mathbb{A}$ et :

$$(I-L)^{-1} = I+L+...+L^n+...$$

Si $L=H$ est un opérateur de $\mathbb{A}$ défini par le noyau $H(x,y)$ c'est-à-dire tel que :

$$(Hf)(x) = \int_{0}^{1} H(x,y) f(y)dy,$$

on démontre que :

$$||H|| = \sup_{x \in [0,1]} \int_{0}^{1} |H(x,y)| dy.$$

Donc si :

$$\sup_{x \in [0,1]} \int_{0}^{1} H(x,y) < 1 \quad (2)$$

ce que nous supposerons réalisé, l'équation de Fredholm, qui s'écrit en termes d'opérateurs :

$$(I-H)f = g$$

a une solution unique qui est :
$$f = (I+H+.. + H^n+..) g.$$

Posons :

$$f_n = (I+H+.. + H^n) g \quad (3)$$

On a :

$$|f - f_n| \le |g| \frac{[||H||]^{n+1}}{1 - ||H||}$$

ce qui donne une majoration de l'erreur faite sur f en se limitant à n termes. La méthode de Monte-Carlo consiste à chercher un jeu dont chaque partie donne une détermination d'une variable aléatoire Z, dont la valeur moyenne soit $f_n(x_o)$, valeur approchée de $f(x_o)$.

## II - PRINCIPE DE LA METHODE -

Nous allons construire un opérateur K et une fonction $g_1$, respectivement à partir de H et g, de façon à avoir :

$$(I+H.. + H^n)g = K^{n+1}g_1 \quad \text{pour} \quad 0 \le x < 1$$

où K est un opérateur linéaire sur $\mathcal{C}'([0,2])$ espace des fonctions continues sur les intervalles semi-ouverts $[0,1[\text{et}] 1,2]$. $\varphi(x)$ étant une fonction quelconque de $\mathcal{C}'([0,2])$; définissons K par :

$$\begin{cases} K \varphi = (H+I_1)\varphi & \text{pour } 0 \le x < 1 \\ K \varphi = I\varphi = \varphi & \text{pour } 1 < x \le 2 \end{cases}$$

où $I_1$ est l'opérateur translation : $(I_1 f)(x) = f(x+1)$.

Par récurrence sur n on trouve immédiatement :

$$\begin{cases} K^{n+1} \varphi = [H^{n+1} + (I+H+..+H^n)I_1]\varphi & \text{pour } 0 \le x < 1 \\ K^{n+1} \varphi = I\varphi = \varphi & \text{pour } 1 < x \le 2 \end{cases}$$

Prenons pour fonctions $\varphi(x)$ la fonction $g_1(x)$ définie par :

$$\begin{cases} g_1(x) = 0 & \text{pour } 0 \le x < 1 \\ g_1(x) = g(x-1) & \text{pour } 1 < x \le 2. \end{cases}$$

On a donc pour $0 \le x < 1$
$$K^{n+1} g_1 = [H^{n+1} + (I+H+\dots+H^n) I_1] g = (I+H+\dots+H^n) g$$

car :

$$I_1 g_1(x) = g(x) \quad \text{pour} \quad 0 \le x < 1$$

ainsi :

$$f_n = K^{n+1} g_1 \quad \text{pour} \quad 0 \le x < 1. \quad (4)$$

### III - CONSTRUCTION D'UN JEU DONNANT UNE VARIABLE ALEATOIRE DE MOYENNE ($K^{n+1} g_1$) ($x_o$).

Donnons-nous (n+2) variables aléatoires en chaîne $A_0, A_1, \dots, A_{n+1}$ pouvant prendre des valeurs $0 \le A_i \le 2$ et définisons la chaîne par une densité de probabilité de passage $\mu(x,y)$, c'est-à-dire :

$$\mathrm{Pr} [y \le A_{k+1} \le y+dy \mid A_k = x] = \mu(x,y) dy = d\bar{\omega}_x(y),$$

où $\bar{\omega}_x$ est, pour $x$ fixé, une mesure positive de masse totale 1, ayant pour support l'intervalle $[0,2]$ :

$$\int_{0}^{2} d\bar{\omega}_{x}(y) = 1 \quad \text{pour tout} \quad x \in [0,2].$$

Revenons à l'opérateur K défini plus haut. Désignons par $\delta(y-\alpha)$ la distribution de. Dirac, opérateur qui à $f(y)$ fait correspondre la constante $f(\alpha)$, ce que l'on peut écrire symboliquement

$$f(\alpha) = \int_{a_1}^{a_2} \delta(y-\alpha) f(y)dy$$

où les bornes d'intégration encadrent strictement $\alpha$ ; alors on peut considérer K comme défini symboliquement par un noyau K(x,y) avec :

$$K(x,y) = H(x,y) \quad \text{dans le domaine} \quad A \quad \begin{cases} 0 \le x < 1 \\ 0 \le y < 1 \end{cases}$$

$$K(x,y) = \delta[y-1-x] \quad \text{dans le domaine} \quad B \quad \begin{cases} 0 \le x < 1 \\ 1 < y \le 2 \end{cases}$$

$$K(x,y) = 0 \quad \text{dans le domaine} \quad C \quad \begin{cases} 1 < x \le 2 \\ 0 \le y < 1 \end{cases}$$

$$K(x,y) = \delta(y-x) \quad \text{dans le domaine} \quad D \quad \begin{cases} 1 < x \le 2 \\ 1 < y \le 2 \end{cases}$$
<figure><img src="image_3.png" /><figcaption>Fig. 1</figcaption></figure>

Cherchons une fonction de 2 variables z (x,y) définie dans le carré [0,2] x [0,2] sauf peut-être en certains points formant un ensemble de mesures de Lebesgue (à 2 dimensions) nulle, et vérifiant en dehors de ces points :

$$z(x,y) \cdot \mu(x,y) = K(x,y). \quad (5)$$

Considérons la mesure $V_x$ associée à la distribution $K_x(y) = K(x,y)$ telle que :

$$(Kf)(x) = \int_{0}^{2} f(y) dV_{x}(y).$$

On a donc :

$$Z(x,y) d\omega_x(y) = dV_x(y) \quad (6)$$

Introduisons la variable aléatoire :

$$Z = z(A_o, A_1, \dots, z_n(A_n, A_{n+1}), g_1(A_{n+1})) \quad (7)$$

associée à la chaîne $A_o A_1 \dots A_{n+1}$, et montrons que sa valeur moyenne conditionnelle est $(K^{n+1} g_1)(x_o)$ lorsque l'on a fixé $A_o = x_o$.

$$E(Z | A_o = x_o) = \int_0^2 \dots \int_0^2 z(x_o, x_1) \dots z(x_n, x_{n+1}) g_1(x_{n+1}) \, d\overline{\omega}_{x_o}(x_1) \dots \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad d\omega_{x_n}(x_{n+1}).$$

Donc d'après (6) :

$$E(Z | A_o = x_o) = \int_0^2 \dots \int_0^2 g_1(x_{n+1}) dV_{x_o}(x_1) \dots dV_{x_n}(x_{n+1}) \quad (8)$$
or on a :

$$ (K g_1) (x) = \int_0^2 g_1(y) \, dV_x(y); $$

il est immédiat par récurrence sur n que :

$$ (K^{n+1} g_1) (x) = \int_0^2 \dots \int_0^2 g_1(x_{n+1}) \, dV_x(x_1) \dots dV_x(x_{n+1}) $$

et d'après (8) :

$$ E(Z | A_o = x_o) = (K^{n+1} g_1) (x_o) = f_n(x_o) $$

si l'on a fixé $A_o = x_o$ avec $0 \le x_o < 1$.

## IV - CHOIX DU JEU ET DEROULEMENT D'UNE PARTIE -

La détermination de $z(x,y)$ et $\mu(x,y)$, liées par :

$$ z(x,y) \mu(x,y) = K(x,y) $$

laisse place à un choix. Nous allons exposer deux solutions possibles. La première est plus facile à programmer que la seconde, mais il n'est pas certain que ce soit la meilleure ou la plus simple possible.

### Première solution :

s étant un nombre choisi entre 0 et 1 (par exemple $\frac{1}{2}$) nous prenons :

$$ \begin{array}{ll} \mu(x,y) = (1-S) & \text{dans } A \\ \mu(x,y) = S \delta(y-1-x) & \text{dans } B \\ \mu(x,y) = 0 & \text{dans } C \\ \mu(x,y) = \delta(y-x) & \text{dans } D \end{array} \quad \text{et} \quad \left\{ \begin{array}{l} z(x,y) = \frac{1}{1-S} H(x,y) \text{ dans } A \\ z(x,y) = \frac{1}{S} \text{ pour } y-1-x = 0 \\ z(x,y) = 1 \text{ ailleurs} \end{array} \right. $$
<figure><img src="image_4.png" /><figcaption>Fig. 2</figcaption></figure>

<figure><img src="image_5.png" /></figure>

Une partie consiste à choisir $A_o = x_o$ avec $0 \le x_o < 1$, à tirer la chaîne $A_o, A_1, \dots, A_{n+1}$, et à calculer la détermination :

$$Z = z(A_o, A_1) \dots z(A_n, A_{n+1}) g_1(A_{n+1})$$

dont on fera la moyenne pour un grand nombre de parties. Le premier tirage se fait suivant la loi de probabilité de densité $\mu(A_o, y)$. Deux cas se présentent :

α) probabilité S de tirer la valeur $A_1 = A_o + 1$.

β) probabilité (1-S) équirépartie sur (0,1) de tirer $A_1 \in (0,1)$.

- si (α) se produit la partie est terminée, car la loi $\mu(A_o + 1, y) = \delta(y - A_o - 1)$ du tirage de $A_2$ est la certitude de tirer $A_2 = A_o + 1$ etc. $A_1 = A_2 = \dots = A_{n+1} = A_o + 1$ ; la valeur de Z est alors :

$$Z = z(A_o, A_o + 1) g_1(A_o + 1) = \frac{1}{S} g(A_o).$$

- si (β) se produit, on fait le deuxième tirage avec la loi $\mu(A_1, y)$ analogue à la précédente. Deux cas sont possibles :

α) tirage de $A_2 = A_1 + 1$ (probabilité S).

β) tirage de $A_2 \in (0,1)$ (probabilité (1-S) équirépartie sur (0,1)).

- si (α) se produit, la partie est terminée avec :

$$Z = \frac{1}{S} z(A_o, A_1) g(A_1);$$

- Si (β) se produit, on fait un nouveau tirage, etc.
De façon générale si la chaîne $A_0, A_1, ... A_p$ avec $p \le n$ n'est pas sortie de l'intervalle $[0,1)$ et si le $(p+1)$ième tirage donne $A_{p+1} = A_p + 1$, la partie est terminée avec :

$$Z = \frac{1}{S} Z(A_0, A_1), Z(A_{p-1}, A_p) g(A_p).$$

La partie peut ne se terminer qu'au $(n+1)$ième tirage :

- si $A_{n+1} = 1 + A_n$ la détermination est :

$$Z = \frac{1}{S} Z(A_0, A_1) .. Z(A_{n-1}, A_n) g(A_n),$$

- si au contraire $A_{n+1}$ est encore dans l'intervalle $[0,1)$ la détermination est $Z = 0$ car $g_1(A_{n+1}) = 0$.

Il resterait à résoudre le problème du choix de S en fonction de $H(x,y)$, $g(x)$ et $n$, pour obtenir la meilleure convergence.

**Seconde solution :**

Un autre type de jeu consiste à prendre :

$$ \left\{ \begin{array}{ll} \mu(x,y) = |H(x,y)| & \text{dans } A \\ \\ \mu(x,y) = [1-\alpha(x)]\delta(y-1-x) & \text{dans } B \\ \\ \mu(x,y) = 0 & \text{dans } C \\ \\ \mu(x,y) = \delta(y-x) & \text{dans } D \end{array} \right. \quad \text{et} \quad \left\{ \begin{array}{l} z(x,y) = \pm 1 \quad \text{signe } H(x,y) \text{ dans } A \\ \\ z(x,y) = \frac{1}{1-\alpha(x)} \quad \text{pour } y-1-x=0 \\ \\ z(x,y) = \pm 1 \quad \text{ailleurs,} \end{array} \right. $$

où l'on a posé $\alpha(x) = \int_{0}^{1} |H(x,y)| dy$. D'après l'hypothèse faite on est sûr que $o \le \alpha(x) < 1$

<figure><img src="image_6.png" /><figcaption>Fig. 3</figcaption></figure>
Le déroulement du jeu est analogue, mais la programmation est plus difficile.

Pour conclure nous dirons que ces méthodes de Monte-Carlo encore très récentes, tributaires du développement des puissantes machines à calculer actuelles, demandent une expérimentation systématique et l'étude statistique des résultats. Il sera alors possible de se rendre compte de leurs possibilités et limites, et de dégager des critères d'optimation du choix du jeu pour une bonne convergence.

## V - APPLICATION NUMÉRIQUE -

Proposons nous d'appliquer la méthode précédente à la résolution de l'équation intégrale :

$$f(x) = 2x + \int_{0}^{1} H(x,y) f(y) dy$$

où :

$$H(x,y) = \begin{cases} 2y(1-x) & \text{si } y < x \\ 2x(1-y) & \text{si } y \ge x \end{cases}$$

La solution formelle de cette équation s'obtient sans difficulté :

$$f(x) = 2 \frac{\sin \sqrt{2} \cdot x}{\sin \sqrt{2}}$$

Nous serons donc en mesure d'apprécier la rapidité de convergence et la précision de la méthode utilisée.

### Première solution :

Le tirage de la chaîne $A_1, A_2, \dots, A_{n+1}$ se fait sans difficulté : on compare un nombre $R_1$, uniformément distribué sur l'intervalle $(0,1)$, au nombre $s$ ; si $R_1 \le s$ la partie est terminée et on a :

$$A_1 = A_2 = \dots = A_{n+1} = A_0 + 1$$

si $R_1 > s$ on tire un deuxième nombre au hasard $R_2$ et en fait :

$$A_1 = R_2$$
On prend alors un nouveau nombre au hasard $R_1$ que l'on compare à s, etc.

On aboutit ainsi à l'organigramme I dont la programmation sur une machine électronique se fait sans difficulté.

Les nombres au hasard $R_1$ et $R_2$ peuvent être engendrés de différentes manières ; dans notre exemple numérique, nous avons adopté un processus multiplicatif défini par la relation :

$$ \lambda_1 = 5^{13} \times \lambda_{1-1} \pmod{2^{35}} $$

A partir de deux nombres $\lambda_0$ choisis de façon quelconque nous avons engendré deux suites de nombres $\lambda_1$ ; les nombres $R_1$ sont ceux de la première suite et les nombres $R_2$ ceux de la seconde (1).

**Deuxième solution :**

Nous avons :

$$ \alpha(x) = \int_{0}^{1} H(x,y) dy = 1 - x(1-x) $$

Le tirage de la chaîne $A_1, A_2, \dots, A_{n+1}$ est un peu plus compliqué : on compare d'abord un nombre $R_1$, uniformément distribué sur l'intervalle $(0,1)$, au nombre $A_0 (1 - A_0)$ ; si $R_1 \ge A_0 (1 - A_0)$ la partie est terminée et on a :

$$ A_1 = A_2 = \dots = A_{n+1} = A_0 + 1 $$

Dans le cas contraire, le nombre $A_1$ est un nombre distribué sur l'intervalle $(0,1)$, avec la densité $H(A_0,y)$ (Figure 4).

Pour trouver $A_1$ nous allons appliquer la méthode bien connue, dite de "REJECTION", dont le principe consiste à choisir au hasard un point dans le carré O P Q S, ce choix étant effectué suivant une distribution uniforme. Si ce point tombe en dessous du graphe de $H(A_0,y)$ on l'accepte et on prend son abscisse comme représentée :

(1) On notera que la distribution des $\lambda_i$ n'a rien de statistique. C'est une distribution arithmétique bien définie, dont les propriétés sont mal connues. Il est "très probable" qu'elle est équirépartie au sens arithmétique, mais cette propriété, supposée démontrée, n'entraîne pas à elle seule que le procédé employé soit valable. Ce procédé est cependant utilisé couramment, et l'expérience montre que les vérifications qu'on a pu faire soit satisfaisantes.
ORGANIGRAMME 1

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 4</figcaption></figure>

tative de la distribution H (A₀,y) ; sinon on rejette le point et on recommence.

Analytiquement cela peut se faire en prenant deux nombres au hasard R₂ et R₃ uniformément distribués sur l'intervalle (0,1) et en comparant R₃ à H (A₀, R₂).

Si R₃ < H (A₀, R₂) le point est en dessous du graphe, en pose

$$A_1 = R_2$$

Sinon on recommence le tirage des 2 nombres R₂ et R₃...

Ayant obtenu A₁, on tire un nouveau nombre R₁ que l'on compare à A₁(1 - A₁), etc.

L'organigramme du calcul est alors l'organigramme II.

Les nombres au hasard R₁, R₂, R₃ sont extraits de 3 suites de nombres λ₁, fabriqués suivant le processus multiplicatif défini dans la première partie.

Les calculs correspondants aux deux solutions ont été exécutés sur un Ordinateur IBM 704 pour x₀ = 0,125.

On trouvera dans les pages suivantes :

1) La représentation graphique des résultats obtenus au cours de 100 000 épreuves pour diverses valeurs de n :

a) Première solution (s = 0,5) page 258

b) Deuxième solution page 259
<figure><img src="image_9.png" /><figcaption>ORGANIGRAMME 2</figcaption></figure>
Les valeurs exactes, calculées avec 4 chiffres significatifs, sont :

$$f_1 (0,125) = 0,3320$$

$$f_2 (0,125) = 0,3510$$

$$f_3 (0,125) = 0,3550$$

$$f_4 (0,125) = 0,3559$$

$$f_5 (0,125) = 0,3560$$

$$f_\infty (0,125) = 0,3561 (0,356071)$$

Le temps de calcul est de l'ordre de 3 à 4 secondes pour 1 000 épreuves.

Nous avons ensuite fait varier la valeur du paramètre s dans la première solution. Les résultats de calcul, poursuivis jusqu'à 500 000 épreuves, sont représentés graphiquement aux pages 260,261, 262 et 263 respectivement pour s = 0,2 - 0,5 - 0,8 et 0,9.

On constate d'importantes fluctuations lorsque s est voisin des bornes de son intervalle de définition (0,2 et 0,9). Par contre la convergence est obtenue rapidement pour n = 0,5 et 0,8. C'est d'ailleurs ce que confirme le calcul de la variance :

$$\sigma^2 = \frac{1}{N-1} \left[ \sum_{i=1}^{N} (Z_i - \bar{Z})^2 \right]$$

En effet celle-ci est égale à :

0,24 pour s = 0,2

0,038 pour s = 0,5

0,036 pour s = 0,8

0,17 pour s = 0,9

Le graphique de la page 264 montre la variation de $\sigma^2$ en fonction de s. Il met en évidence un minimum de l'ordre de 0,016 pour s = 0,7.

En appliquant la formule connue :

$$\text{Prob} [|\bar{F} - \langle F \rangle| > h] \le \frac{\sigma^2}{nh^2}$$
on voit que l'on a environ 84 chances sur 100 d'obtenir une précision supérieure à 1 % en 1 000 épreuves, c'est-à-dire en 3 sec, 5 de calcul environ (une expérience a donné 0,3531, valeur proche de la solution exacte : 0,3561).

D'autre part, pour des raisons évidentes de rapidité, on a intérêt à choisir une valeur de s la plus grande possible. Nous avons représenté en page 265 la variation du temps de calcul de 10 000 épreuves en fonction de s. La forme de la courbe montre, qu'à partir de s = 0,7, le gain de temps obtenu en augmentant s est très peu appréciable.

Dans ces conditions la valeur s = 0,7 est donc particulièrement favorable.

A titre de comparaison, signalons que la variance $σ^2$ est égale à 0,083 par la deuxième méthode.

## REFERENCES

J. BASS - Nombres aléatoires, suites arithmétiques, méthode de Monte-Carlo.
Publ. de l'Inst. de statistique de l'Université de Paris IX, 1960, p. 289-325.

J.P. BERTRANDIAS - Calcul d'une intégrale au moyen de la suite $X_n = \text{An mod. 1.}$ Evaluation de l'erreur. Publ. de l'Inst. de statistique de l'Université de Paris, IX, 1960, p. 335-357.

J.H. CURTISS - "Monte-Carlo" Methods for iteration of lineare operators J. of Mathematics and Physics, 32, 1954, p. 209-232.

Symposium on Monte-Carlo Methods (J. Wiley and sons, 1954) - Articles de GE. Albert et J.H. Curtiss.
<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /></figure>

Variation du temps de calcul en fonction du paramètre $\delta$