8.6 Analyse numérique des fonctions d’une variable

8.6.1 Interpolation linéaire, interpolation polynomiale

Considérons f une fonction de classe {C}^{n} sur l’intervalle [a,b], soit {x}_{1} < \mathop{\mathop{…}} < {x}_{n} des points de [a,b] et considérons l’unique polynôme P ∈ {ℝ}_{n−1}[X] vérifiant P({x}_{i}) = f({x}_{i}) pour 1 ≤ i ≤ n (polynôme d’interpolation de Lagrange).

Lemme 8.6.1 Pour tout x ∈ [a,b], \mathop{∃}ζ ∈]a,b[, f(x) − P(x) ={ (x−{x}_{1})\mathop{\mathop{…}}(x−{x}_{n}) \over n!} {f}^{(n)}(ζ).

Démonstration Si x est l’un des {x}_{i}, n’importe quel ζ convient. On peut donc supposer que f n’est pas l’un des {x}_{i}. Considérons alors g : t\mathrel{↦}f(t) − P(t) − λ{ (t−{x}_{1})\mathop{\mathop{…}}(t−{x}_{n}) \over n!} λ est choisi de telle sorte que g(x) = 0 (c’est évidemment possible). Alors g est de classe {C}^{n} et a n + 1 zéros distincts sur l’intervalle [a,b]. Des applications répétées du théorème de Rolle assurent que la fonction g' a n zéros distincts sur l’intervalle ]a,b[ (un entre deux zéros de g au sens strict), puis que g'' a n − 1 zéros distincts jusque {g}^{(n)} qui a au moins un zéro. Soit donc ζ tel que {g}^{(n)}(ζ) = 0. On a donc 0 = {g}^{(n)}(ζ) = {f}^{(n)}(ζ) − λ car {P}^{(n)} = 0 (vu que \mathop{deg} P ≤ n − 1) et { {d}^{n} \over d{t}^{n}} ((t − {x}_{1})\mathop{\mathop{…}}(t − {x}_{n})) = n!. on a donc λ = {f}^{(n)}(ζ) et en écrivant que g(x) = 0, on obtient f(x) − P(x) ={ (x−{x}_{1})\mathop{\mathop{…}}(x−{x}_{n}) \over n!} {f}^{(n)}(ζ).

Considérons le cas particulier où n = 2 et où l’on prend {x}_{1} = a et {x}_{2} = b ; c’est le cas de l’interpolation linéaire entre a et b où l’on remplace la courbe y = f(x) par la corde joignant les points (a,f(a)) et (b,f(b)). On a alors P(t) = f(a) +{ f(b)−f(a) \over b−a} (t − a). Si on appelle {M}_{2} ={\mathop{ sup}}_{t∈[a,b]}|f''(t)|, on obtient

Proposition 8.6.2 (erreur dans une interpolation linéaire). Soit f : [a,b] → ℝ de classe {C}^{2}, {M}_{2} ={\mathop{ sup}}_{t∈[a,b]}|f''(t)|. Alors

\mathop{∀}t ∈ [a,b], |f(t) −\left (f(a) +{ f(b) − f(a) \over b − a} (t − a)\right )|≤ {M}_{2}{ {(b − a)}^{2} \over 8}

Démonstration On a en effet

\begin{eqnarray*} |f(x) − P(x)|& =&{ (x − a)(b − x) \over 2} |f''(ζ)|≤ {M}_{2}{ (x − a)(b − x) \over 2} %& \\ & ≤& {M}_{2}{ {(b − a)}^{2} \over 8} %& \\ \end{eqnarray*}

car (x − a)(b − x) ≤{ {(b−a)}^{2} \over 4} pour x ∈ [a,b].

8.6.2 Dérivation numérique

Nous nous limiterons au calcul approché des dérivées d’ordre 1 et 2 d’une fonction numérique. Pour la dérivée d’ordre 1 nous utiliserons la formule f(x + h) = f(x) + hf'(x) + hε(h) avec {\mathop{lim}}_{h→0}ε(h) = 0. On en déduit que f'(x) est peu différent de {Δ}_{h}f(x) ={ f(x+h)−f(x) \over h} quand h est petit. Mathématiquement, plus h est petit, plus {Δ}_{h}f(x) est proche de f'(x). Mais qu’en est-il de la valeur calculée \overline{{Δ}_{h}f(x))} ? Le calcul de f(x + h) − f(x) se fait avec une erreur absolue de l’ordre de 2δf(x), où δ est la précision de l’instrument de calcul ( 1{0}^{−16} par exemple). On a donc |{Δ}_{h}f(x) −\overline{{Δ}_{h}f(x))}|≤{ 2δ|f(x)| \over h} , qui tend vers + ∞ quand h tend vers 0. Il faut donc trouver un compromis pour la valeur de h à choisir. Supposons f de classe {C}^{2}. Alors on a f(x + h) = f(x) + hf'(x) +{ {h}^{2} \over 2} f''(x) + {h}^{2}ε(h) et donc |f'(x) − {Δ}_{h}f(x)| est peu différent de { h \over 2} |f''(x)|. On a donc |f'(x) −\overline{{Δ}_{h}f(x))}|≤{ h \over 2} |f''(x)| +{ 2δ|f(x)| \over h} . Le deuxième membre est une fonction de h qui est minimale pour h = 2\sqrt{{ δ|f(x)| \over |f''(x)|} } et qui vaut alors 2\sqrt{δ|f(x)f''(x)|}. Avec les fonctions usuelles, on retiendra qu’il faut choisir un h de l’ordre de \sqrt{δ} (plutôt un peu trop grand, qu’un peu trop petit) et que l’on obtient alors une erreur de l’ordre de \sqrt{ δ}. Ainsi on choisira par exemple h = 1{0}^{−8} et on pourra espérer avoir 7 ou 8 chiffres significatifs dans le calcul de la dérivée (ce qui est le plus souvent largement suffisant, par exemple pour une étude de fonction).

Une autre méthode de calcul de la dérivée qui donne des valeurs un peu plus précises (mais peut poser des problèmes de définition de la fonction aux bornes de l’intervalle) est de prendre comme valeur approchée de la dérivée l’expression { f(x+h)−f(x−h) \over 2h} (dérivée symétrique). L’erreur est alors en { {h}^{2} \over 6} |{f}^{(3)}(x)| +{ 2δ|f(x)| \over h} , elle est minimale pour un h de l’ordre de \root{3}\of{δ} et elle est alors de l’ordre de {δ}^{2∕3} (prendre h = 1{0}^{−5} pour obtenir environ 9 à 10 chiffres significatifs).

Pour le calcul de la dérivée seconde, reprenons la formule de Taylor-Young f(x + h) = f(x) + hf'(x) +{ {h}^{2} \over 2} f''(x) + {h}^{2}ε(h). On obtient alors {\mathop{lim}}_{h→0}{ f(x+h)+f(x−h)−2f(x) \over {h}^{2}} = f''(x). On utilisera comme valeur approchée de f''(x) l’expression {Δ}_{h}^{(2)}f(x) ={ f(x+h)+f(x−h)−2f(x) \over {h}^{2}} . Utilisons la même méthode pour évaluer l’erreur entre la valeur calculée de {Δ}_{h}^{(2)}f(x) et f''(x). Supposons f de classe {C}^{4}. On a la formule de Taylor Young à l’ordre 4, f(x + h) = f(x) + hf'(x) +{ {h}^{2} \over 2} f''(x) +{ {h}^{3} \over 6} {f}^{(3)}(x) +{ {h}^{4} \over 24} {f}^{(4)}(x) + {h}^{4}ε(h) d’où |{Δ}_{h}^{(2)}f(x) − f''(x)| est peu différent de { {h}^{2} \over 12} |{f}^{(4)}(x)|. On obtient donc une majoration du type |\overline{{Δ}_{h}^{(2)}f(x)} − f''(x)|≤{ {h}^{2} \over 12} |{f}^{(4)}(x)| +{ 3δ \over {h}^{2}} |f(x)|. L’erreur est minimale pour une valeur de h de l’ordre de \root{4}\of{δ}. On choisira donc un h de l’ordre de 1{0}^{−4} pour obtenir un résultat optimal pour les fonctions usuelles.

8.6.3 Recherche des zéros d’une fonction

On suppose dans toutes les méthodes qui suivent que l’on a effectué la séparation des zéros de la fonction f, c’est-à-dire que l’on a trouvé un intervalle [a,b] sur lequel f est strictement monotone avec f(a)f(b) < 0. On suppose également que f est suffisamment dérivable. On sait alors que f a un unique zéro r sur l’intervalle [a,b].

Méthode de dichotomie

Soit c un point de ]a,b[. Alors, soit f(a) et f(c) sont de signe contraire, auquel cas r ∈]a,c[, soit f(a) et f(c) sont de même signe et dans ce cas r ∈]c,a[. En prenant c ={ a+b \over 2} et en itérant le procédé, on construit une suite de segments emboîtés [{a}_{n},{b}_{n}] tels que \mathop{∀}n ∈ ℕ, r ∈ [{a}_{n},{b}_{n}] et {b}_{n} − {a}_{n} ={ 1 \over 2} ({b}_{n−1} − {a}_{n−1}), soit {b}_{n} − {a}_{n} ={ 1 \over {2}^{n}} (b − a). Le théorème des segments emboîtés garantit alors que {\mathop{\mathop{⋂ }} }_{n∈ℕ}[{a}_{n},{b}_{n}] = r. Si l’on prend {a}_{n} comme valeur approchée de r par exemple, on a |r − {a}_{n}|≤{ 1 \over {2}^{n}} (b − a).

Méthode de Lagrange

La méthode de Lagrange consiste à assimiler sur le segment [a,b] la courbe y = f(x) à la droite passant par les points (a,f(a)) et (b,f(b)), c’est-à-dire à approcher f par la fonction P(x) = f(a) +{ f(b)−f(a) \over b−a} (x − a) et à prendre comme valeur approchée de r le réel \bar{r} vérifiant P(\bar{r}) = 0.

Etudions l’erreur commise dans cette méthode. Soit P(x) = f(a) +{ f(b)−f(a) \over b−a} (x − a) (interpolation linéaire). La majoration de l’erreur dans une interpolation linéaire nous garantit que, en posant {M}_{2} ={\mathop{ sup}}_{t∈[a,b]}|f''(t)|, on a |f(r) − P(r)|≤ {M}_{2}{ {(b−a)}^{2} \over 8} . Mais f(r) = 0 et { P(r) \over r−\bar{r}} ={ P(r)−P(\bar{r}) \over r−\bar{r}} ={ f(b)−f(a) \over b−a} = f'(c) pour un c ∈]a,b[. On en déduit que |\bar{r} − r| = \left |{ P(r) \over f'(c)} \right |≤{ {M}_{2}{(b−a)}^{2} \over 8{m}_{1}} si l’on pose {m}_{1} ={\mathop{ inf} }_{t∈[a,b]}|f'(t)| (que nous supposerons strictement positif). D’où la majoration de l’erreur dans la méthode de Lagrange

|\bar{r} − r|≤{ {M}_{2}{(b − a)}^{2} \over 8{m}_{1}}

On pourra par exemple combiner la méthode de dichotomie et la méthode de Lagrange pour trouver rapidement une bonne approximation de r. Remarquons de plus que si on suppose en outre que f est de concavité constante sur [a,b], alors on connait le signe de r −\bar{ r}. En effet

(f convexe croissante, ou f concave décroissante)) ⇒\bar{ r} < r

(f convexe décroissante, ou f concave croissante)) ⇒\bar{ r} > r

Méthode de Newton

La méthode de Newton consiste à assimiler la courbe y = f(x) à la tangente en un point c ∈ [a,b]. Cette tangente a pour équation y − f(c) = f'(c)(x − c). Elle coupe donc l’axe des x au point d’abscisse r' = c −{ f(c) \over f'(c)} .

Cherchons à majorer l’erreur commise |r − r'|. On a d’après la formule de Taylor Lagrange 0 = f(r) = f(c) + (r − c)f'(c) +{ {(r−c)}^{2} \over 2} f''(x) pour un certain x ∈]r,c[. De plus on a f(c) + (r' − c)f'(c) = 0. En soustrayant membre à membre les deux égalités on trouve (r' − r)f'(c) −{ {(r−c)}^{2} \over 2} f''(x) = 0, soit r' − r ={ {(r−c)}^{2} \over 2} { f''(x) \over f'(c)} et donc (avec les notations données ci dessus dans la méthode de Lagrange)

|r' − r|≤{ {M}_{2}{(b − a)}^{2} \over 2{m}_{1}}

Remarquons que la majoration de l’erreur obtenue est 4 fois plus grande que dans la méthode de Lagrange. Ce n’est évidemment pas que la méthode de Newton soit moins bonne que la méthode de Lagrange, c’est simplement la majoration de l’erreur qui est un peu moins fine. Remarquons que si l’on suppose en outre que f est de concavité constante sur [a,b], alors on connait le signe de r' − r. En effet

(f convexe croissante, ou f concave décroissante)) ⇒ r' > r

(f convexe décroissante, ou f concave croissante)) ⇒ r' < r

Les inégalités sont donc en sens contraire de celles obtenues par la méthode de Lagrange. Dans le cas où f est monotone et de concavité constante sur [a,b], la combinaison de la méthode de Lagrange et de la méthode de Newton fournit un encadrement de r, ce qui est meilleur qu’une majoration d’erreur.

Dans la pratique on ne s’arrête pas après avoir appliqué une fois la méthode de Newton, mais au contraire on applique de nouveau la méthode de Newton, mais cette fois ci au point r'. Cela revient à construire une suite ({x}_{n}) par récurrence par {x}_{o} = c et {x}_{n+1} = {x}_{n} −{ f({x}_{n}) \over f'({x}_{n})} . Les questions qui se posent naturellement sont

Il est clair que dans la mesure où les réponses aux questions (i) et (ii) sont oui, la réponse à la question (iii) est aussi oui, puisque si l’on appelle L la limite de la suite, on doit avoir L = L −{ f(L) \over f'(L)} , soit f(L) = 0. Il nous reste donc à répondre aux deux premières questions.

Nous nous placerons sous les hypothèses suivantes : f est de classe {C}^{2} sur [a,b], f' ne s’annule pas sur [a,b] et f'' est de signe constant sur [a,b] (donc f est strictement monotone et de concavité constante). Dans un premier temps nous supposerons pour faire les raisonnements que f est croissante convexe sur [a,b]. Prenons {x}_{o} = c > r (par exemple c = b). On va alors montrer que \mathop{∀}n, {x}_{n} ∈ [r,b] et que la suite ({x}_{n}) est décroissante, ce qui permettra de répondre positivement aux questions (i) et (ii). Pour cela considérons la fonction g définie par g(x) = x −{ f(x) \over f'(x)} ={ xf'(x)−f(x) \over f'(x)} . On a g'(x) ={ f(x)f''(x) \over f'{(x)}^{2}} ≥ 0 donc g est croissante sur [a,b]. Supposons que {x}_{n} ∈ [r,b]. Comme g(r) = r, on a {x}_{n+1} = g({x}_{n}) ≥ g(r) = r. De plus f étant strictement croissante, elle est positive sur [r,b] et donc {x}_{n+1} = {x}_{n} −{ f({x}_{n}) \over f'({x}_{n})} ≤ {x}_{n}, d’où r ≤ {x}_{n+1} ≤ {x}_{n} ≤ b. On en déduit que les réponses aux questions (i) et (ii) sont positives et fournissent donc un moyen d’approximation de r. Remarquons qu’il est fondamental pour cela d’avoir choisi un c > r, car g est décroissante sur [a,r] et donc si {x}_{o} = c < r, {x}_{1} = g({x}_{o}) ≥ g(r) = r, mais plus rien ne garantit que {x}_{1} appartient toujours à [a,b]. Dans les cas de monotonie ou de concavité différents on a les conclusions suivantes

(f convexe croissante ou f concave décroissante) : choisir {x}_{o} > r ; la suite ({x}_{n}) est décroissante et converge vers r

(f convexe décroissante ou f concave croissante) : choisir {x}_{o} < r ; la suite ({x}_{n}) est croissante et converge vers r

Il nous reste à savoir avec quelle vitesse la suite converge. On a g'(r) = 0. Puisque g est continue, si l’on se donne K < 1, il existe h > 0 tel que \mathop{∀}x ∈ [r,r + h], |g'(x)| < K. Alors soit N tel que n > N ⇒ {x}_{n} ∈ [r,r + h]. On a alors pour n > N, |{x}_{n+1} − r| = |g({x}_{n}) − g(r)|≤ K|{x}_{n} − r|, d’où pour n > N, |{x}_{n} − r|≤ {K}^{n−N}|{x}_{N} − r|. On en déduit que la suite ({x}_{n} − r) est négligeable devant la suite {K}^{n} pour tout K < 1, et on a donc une convergence extrêmement rapide dès que l’on se rapproche de r. On peut d’ailleurs trouver un équivalent de |{x}_{n} − r| dans le cas où g''(r)\mathrel{≠}0 et montrer que |{x}_{n} − r|∼ {K}^{{2}^{n} } avec un K < 1. On a donc une convergence très rapide (de type hyperexponentiel) : en gros le nombre de décimales double à chaque itération dans la limite de la précision de la machine.

Méthode des approximations successives

Elle consiste à transformer une équation du type f(x) = 0 en une équation du type g(x) = x à laquelle on essayera d’appliquer une méthode de point fixe : on construit une suite {x}_{n+1} = g({x}_{n}) ; si cette suite converge, elle converge vers un point fixe de g. Le lecteur pourra remarquer que la méthode de Newton en est un cas particulier avec g(x) = x −{ f(x) \over f'(x)} .