CM 2 - Schémas “classiques” à un pas

%reset -f
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
# plt.rcdefaults()
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})
from scipy.optimize import fsolve

Considérons le problème de Cauchy

trouver une fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) définie sur un intervalle \(I=[t_0,T]\) telle que \(\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases}\) avec \(y_0\) une valeur donnée et supposons que l’on ait montré l’existence et l’unicité d’une solution \(y\) pour \(t\in I\).

Pour \(h>0\) soit \(t_n\equiv t_0+nh\) avec \(n=0,1,2,\dots,N\) une suite de \(N+1\) nœuds de \(I\) induisant une discrétisation de \(I\) en \(N\) sous-intervalles \(I_n=[t_n;t_{n+1}]\) chacun de longueur \(h=\frac{T-t_0}{N}>0\) (appelé le pas de discrétisation).

Pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\equiv y(t_n)\).

• L’ensemble de \(N+1\) valeurs \(\{t_0, t_1=t_0+h,\dots , t_{N}=T \}\) représente les points de la discrétisation.
  
• L’ensemble de \(N+1\) valeurs \(\{y_0, y_1,\dots , y_{N} \}\) représente la solution exacte discrète.
  
• L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,\dots , u_{N} \}\) représente la solution numérique obtenue en construisant une suite récurrente.

Les schémas que l’on va construire permettent de calculer (explicitement ou implicitement) \(u_{n+1}\) à partir de \(u_n, u_{n-1}, ..., u_{n-k}\) et il est donc possible de calculer successivement \(u_1\), \(u_2\),…, en partant de \(u_0\) par une formule de récurrence de la forme \(\begin{cases} u_0=y_0,\\ \vdots\\ u_\kappa=y_\kappa,\\ u_{n+1}=\Phi(u_{n+1},u_n, u_{n-1}, \dots, u_{n-k}),&\forall n=\kappa,\kappa+1,\dots,N-1. \end{cases}\)

Méthodes explicites et méthodes implicites

• Une méthode est dite explicite si la valeur \(u_{n+1}\) peut être calculée directement à l’aide des valeurs précédentes \(u_k\), \(k\le n\) (ou d’une partie d’entre elles).
  
• Une méthode est dite implicite si \(u_{n+1}\) n’est défini que par une relation implicite faisant intervenir la fonction \(\varphi\).

Méthodes à un pas et méthodes multi-pas

• Une méthode numérique pour l’approximation du problème de Cauchy est dite à un pas si pour tout \(n\in\mathbb{N}\), \(u_{n+1}\) ne dépend que de \(u_n\) et éventuellement de lui-même.
  
• Autrement, on dit que le schéma est une méthode multi-pas (ou à pas multiples).

1 Construction de schémas à un pas

Si nous intégrons l’EDO \(y'(t)=\varphi(t,y(t))\) entre \(t_n\) et \(t_{n+1}\) nous obtenons \( \int_{t_n}^{t_{n+1}}y'(t)\mathrm{d}t=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t \) c’est-à-dire \( y_{n+1}-y_n=\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t. \) On peut construire différents schémas selon la formule d’approximation utilisée pour approcher le membre de droite. Cette solution approchée sera obtenue en construisant une suite récurrente comme suit:

\(\begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\displaystyle\int_{t_n}^{t_{n+1}} \tilde f(t) \mathrm{d}t \quad \text{où $\tilde f(t)$ est un polynôme interpolant }\varphi(t,y(t)) \end{cases} \)

1.1 Schéma d’Euler explicite

Leonhard Euler

Si on remplace une fonction \(f\) par une constante égale à la valeur de \(f\) à borne gauche de l’intervalle \([a;b]\)
(polynôme qui interpole \(f\) au point \((a,f(a))\) et donc de degré \(0\)), on a \(\begin{align*} \tilde f(x)&=f(a)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(a). \end{align*}\) Cette formule est dite formule de quadrature du rectangle à gauche.

En utilisant cette formule pour approcher la fonction \(t\mapsto \varphi(t,y(t))\) on a \( \int_{t_n}^{t_{n+1}} \varphi(t,y(t)) \mathrm{d}t \approx h \varphi(t_n,y(t_n)) \) et on reconnait le schéma d’Euler progressif

\( \text{(EE)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases}\)

Il s’agit d’un schéma à 1 pas explicite car il permet d’expliciter \(u_{n+1}\) en fonction de \(u_n\).

def EE(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        uu[i+1] = uu[i]+h*phi(tt[i],uu[i])
    return uu

1.2 Schéma d’Euler implicite

Si on remplace une fonction \(f\) par une constante égale à la valeur de \(f\) à borne droite de l’intervalle \([a;b]\)
(polynôme qui interpole \(f\) au point \((b,f(b))\) et donc de degré \(0\)), on a \(\begin{align*} \tilde f(x)&=f(b)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f(b). \end{align*}\) Cette formule est dite formule de quadrature du rectangle à droite.

En utilisant cette formule pour approcher la fonction \(t\mapsto \varphi(t,y(t))\) on a \( \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx (t_{n+1}-t_n) \varphi(t_{n+1},y(t_{n+1})) = h \varphi(t_{n+1},y(t_{n+1})) \) et on obtient le schéma d’Euler rétrograde

\( \text{(EI)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi(t_{n+1},u_{n+1})& n=0,1,2,\dots N-1 \end{cases} \)

Il s’agit d’un schéma à 1 pas implicite car il ne permet pas de calculer directement \(u_{n+1}\) en fonction de \(u_n\) lorsque la fonction \(\varphi\) n’est pas triviale.

Pour calculer \(u_{n+1}\) il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, \(u_{n+1}\) est solution de l’équation \(x=u_n+h\varphi(t_{n+1},x)\), c’est-à-dire un zéro de la fonction (en général non linéaire) \(x\mapsto -x+u_n+h\varphi(t_{n+1},x).\)

def EI(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        tmp = fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i])
        uu[i+1] = tmp[0]
    return uu

1.3 Schémas d’Euler modifié (EM) ou de type Runge-Kutta (RK1_M)

Si on remplace une fonction \(f\) par une constante égale à la valeur de \(f\) au milieu de l’intervalle \([a;b]\)
(polynôme qui interpole \(f\) au point \(\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)\) et donc de degré \(0\)), on a \(\begin{align*} \tilde f(x)&=f\left(\tfrac{a+b}{2}\right)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=(b-a)f\left(\tfrac{a+b}{2}\right). \end{align*}\) Cette formule est dite formule de quadrature du rectangle ou du point milieu.

En utilisant cette formule pour approcher la fonction \(t\mapsto \varphi(t,y(t))\) on a \( \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx h \varphi\left(t_n+\frac{h}{2},y\left(t_n+\frac{h}{2}\right)\right) \) et on obtient \(\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi\left(t_n+\frac{h}{2},u_{n+1/2}\right)& n=0,1,2,\dots N-1 \end{cases}\) où \(u_{n+1/2}\) est une approximation de \(y(t_n+h/2)\).

Cependant, \(u_{n+1/2}\) est une valeur inconnue, on doit alors en calculer une approximation \(\tilde u_{n+1/2}\). Pour cela nous pouvons utiliser deux stratégies:

1. Utiliser une prédiction d’Euler progressive sur un demi-pas: \(\tilde u_{n+1/2}=u_n+(h/2) \varphi(t_{n},u_{n})\).

   On peut ainsi remplacer \(\varphi(t_{n}+h/2,u_{n+1/2})\) par \(\varphi(t_{n}+h/2,\tilde u_{n+1/2})\).

   Nous avons construit ainsi un nouveau schéma appelé schéma d’Euler modifié qui s’écrit comme suit et qui est à 1 pas explicite car il permet d’expliciter \(u_{n+1}\) en fonction de \(u_n\).

\( \text{(EM)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+h \varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

2. Utiliser une moyenne: \(\tilde u_{n+1/2}=\dfrac{u_n+u_{n+1}}{2}\).

   On peut ainsi remplacer \(\varphi(t_{n}+h/2,u_{n+1/2})\) par \(\varphi(t_{n}+h/2,\tilde u_{n+1/2})\).

   Nous avons construit ainsi un nouveau schéma qu’on appelera schéma RK1_M qui s’écrit comme suit et qui est à 1 pas implicite car il ne permet pas d’expliciter \(u_{n+1}\) en fonction de \(u_n\). Pour calculer \(u_{n+1}\) il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, \(u_{n+1}\) est solution de l’équation \(x=u_n+h\varphi(t_{n}+\frac{h}{2},\frac{u_n+x}{2})\), c’est-à-dire un zéro de la fonction (en général non linéaire) \(x\mapsto -x+u_n+h\varphi(t_{n}+\frac{h}{2},\frac{u_n+x}{2}).\)

\(\text{(RK1\_M)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+h \varphi\left(\frac{t_n+t_{n+1}}{2},\dfrac{u_n+u_{n+1}}{2}\right)& n=0,1,2,\dots N-1 \end{cases} \)

def EM(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        u_tilde = uu[i]+0.5*h*phi(tt[i],uu[i])
        uu[i+1] = uu[i]+h*phi(tt[i]+h/2,u_tilde)
    return uu

def RK1_M(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        tmp = fsolve(lambda x: -x+uu[i]+h*phi( (tt[i]+tt[i+1])/2,(uu[i]+x)/2 ), uu[i])
        uu[i+1] = tmp[0]
    return uu

1.4 Schéma du trapèze (ou de Crank-Nicolson)

Si on remplace une fonction \(f\) par le segment qui relie \((a,f(a))\) à \((b,f(b))\)
(polynôme qui interpole \(f\) en les points \((a,f(a))\) et \((b,f(b))\) et donc de degré \(1\)), on a \(\begin{align*} \tilde f(x)&=\dfrac{f(b)-f(a)}{b-a}(x-a)+f(a)\\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=\frac{b-a}{2}\left(f(a)+f(b)\right). \end{align*}\) Cette formule est dite formule de quadrature du trapèze.

En utilisant cette formule pour approcher la fonction \(t\mapsto \varphi(t,y(t))\) on a \( \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{2}\Big(\varphi(t_{n},y(t_{n}))+\varphi(t_{n+1},y(t_{n+1}))\Big) \) et on obtient le schéma du trapèze ou de Crank-Nicolson

\( \text{(CN)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+\dfrac{h}{2} \Big(\varphi(t_{n},u_{n}) + \varphi(t_{n+1},u_{n+1})\Big)& n=0,1,2,\dots N-1 \end{cases} \)

En fait, ce schéma fait la moyenne des schémas d’Euler progressif et rétrograde.

Il s’agit à nouveau d’un schéma à 1 pas implicite car il ne permet pas d’expliciter directement \(u_{n+1}\) en fonction de \(u_n\) lorsque la fonction \(\varphi\) n’est pas triviale. Pour calculer \(u_{n+1}\) il faudra utiliser un schéma pour le calcul du zéro d’une fonction quelconque. En effet, \(u_{n+1}\) est solution de l’équation \(x=u_n+\frac{h}{2}\left(\varphi(t_n,u_n)+\varphi(t_{n+1},x)\right)\), c’est-à-dire un zéro de la fonction (en général non linéaire) \(x\mapsto -x+u_n+\frac{h}{2}\Big(\varphi(t_n,u_n)+\varphi(t_{n+1},x)\Big).\)

def CN(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        tmp = fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i],uu[i])+phi(tt[i+1],x) ), uu[i])
        uu[i+1] = tmp[0]
    return uu

1.5 Schéma de Heun

Karl Heun

Pour éviter le calcul implicite de \(u_{n+1}\) dans le schéma du trapèze, nous pouvons utiliser une prédiction d’Euler progressive et remplacer le \(u_{n+1}\) dans le terme \(\varphi(t_{n+1},u_{n+1})\) par \(\tilde u_{n+1}=u_n+h \varphi(t_{n},u_{n})\). Nous avons construit ainsi un nouveau schéma appelé schéma de Heun. Plus précisément, la méthode de Heun s’écrit

\( \text{(H)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1}=u_n+h \varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+\dfrac{h}{2} \Big(\varphi(t_{n},u_{n}) + \varphi(t_{n+1},\tilde u_{n+1})\Big) & n=0,1,2,\dots N-1 \end{cases} \)

Il s’agit à nouveau d’un schéma à 1 pas explicite.

def heun(phi,tt,y0):
    h  = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i+1], uu[i] + h*k1  )
        uu[i+1] = uu[i] + 0.5*h*(k1+k2)
    return uu

1.6 Remarques

1. À première vue, il semble que les schémas d’Euler progressif et de Heun soient préférable aux schémas d’Euler rétrograde et de Crank-Nicolson puisque ces derniers ne sont pas explicites. Cependant, on verra que les méthodes d’Euler implicite et de Crank-Nicolson sont inconditionnellement A-stables. C’est aussi le cas de nombreuses autres méthodes implicites. Cette propriété rend les méthodes implicites attractives, bien qu’elles soient plus coûteuses que les méthodes explicites.

2. Pour la mise en application d’un schéma il faut aussi prendre en compte l’influence des erreurs d’arrondi. En effet, afin de minimiser l’erreur globale théorique, on pourrait être tenté d’appliquer une méthode avec un pas très petit, par exemple de l’ordre de \(10^{-16}\), mais ce faisant, outre que le temps de calcul deviendrait irréaliste, très rapidement les erreurs d’arrondi feraient diverger la solution approchée. En pratique il faut prendre \(h\) assez petit pour que la méthode converge, mais pas trop petit non plus pour que les erreurs d’arrondi ne donnent pas lieu à des résultats incohérents et pour que les calculs puissent être effectués en un temps raisonnable.

2 Annexes

2.1 Schémas basés sur la formule de quadrature de Simpson

2.1.1 Schéma de Simpson implicite

Si on remplace une fonction \(f\) par la parabole segment qui passe par \((a,f(a))\), \((b,f(b))\) et \(\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)\)
(polynôme qui interpole \(f\) en les points \((a,f(a))\), \((b,f(b))\) et \(\left(\tfrac{a+b}{2},f\left(\tfrac{a+b}{2}\right)\right)\) et donc de degré \(2\)), on a \(\begin{align*} \tilde f(x)&=f(a)\frac{(x-b)\left(x-\tfrac{a+b}{2}\right)}{(a-b)\left(a-\tfrac{a+b}{2}\right)} +f\left(\tfrac{a+b}{2}\right)\frac{(x-a)(x-b)}{\left(\tfrac{a+b}{2}-a\right)\left(\tfrac{a+b}{2}-b\right)} +f(b)\frac{(x-a)\left(x-\tfrac{a+b}{2}\right)}{(b-a)\left(b-\tfrac{a+b}{2}\right)} \\ \int_a^b f(x)\mathrm{d}x &\approx \int_a^b \tilde f(x)\mathrm{d}x=\frac{h}{6} \left( f(a)+4f\left(\tfrac{a+b}{2}\right)+ f(b)\right). \end{align*}\) Cette formule est dite formule de Simpson.

En utilisant cette formule pour approcher la fonction \(t\mapsto \varphi(t,y(t))\) on a \( \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t\approx \frac{h}{6} \left( \varphi(t_n,y(t_n))+4\varphi\left(t_n+\frac{h}{2},y\left(t_n+\frac{h}{2}\right)\right)+ \varphi(t_{n+1},y(t_{n+1})) \right) \) et on obtient
\(\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},u_{n+1/2}\right)+\varphi(t_{n+1},u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases}\) où \(u_{n+1/2}\) est une approximation de \(y(t_n+h/2)\).

Cependant, \(u_{n+1/2}\) est une valeur inconnue, on doit alors en calculer une approximation \(\tilde u_{n+1/2}\). Pour cela nous pouvons utiliser une prédiction d’Euler progressive sur un demi-pas: \(\tilde u_{n+1/2}=u_n+(h/2) \varphi(t_{n},u_{n})\). On peut ainsi remplacer \(\varphi(t_{n}+h/2,u_{n+1/2})\) par \(\varphi(t_{n}+h/2,\tilde u_{n+1/2})\).

Nous avons construit ainsi un nouveau schéma appelé schéma de Simpson implicite qui s’écrit

\( \text{(SI)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)+\varphi(t_{n+1},u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} \)

Il s’agit d’un schéma à 1 pas implicite car il ne permet pas d’expliciter \(u_{n+1}\) en fonction de \(u_n\).

2.1.2 Schéma de Simpson explicite

Pour éviter le calcul implicite de \(u_{n+1}\) dans le schéma de Simpson implicite, nous pouvons utiliser une prédiction d’Euler progressive et remplacer le \(u_{n+1}\) dans le terme \(\varphi(t_{n+1},u_{n+1})\) par \(\check u_{n+1}=u_n+h \varphi(t_{n},u_{n})\). Nous avons construit ainsi un nouveau schéma qu’on appellera schéma de Simpson explicite et qui s’écrit

\( \text{(SE)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ \check u_{n+1}=u_n+h \varphi(t_{n},u_{n}),\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)+\varphi(t_{n+1},\check u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} \)

Il s’agit à nouveau d’un schéma à 1 pas explicite.

2.1.3 Schéma de Simpson explicite-variante

Notons qu’on aurait pu remplacer le \(u_{n+1}\) dans le terme \(\varphi(t_{n+1},u_{n+1})\) par une approximation utilisant \(\tilde u_{n+1/2}\) comme par exemple une prédiction d’Euler progressive à partir de \(t_n+h/2\), ce qui donne \(\hat u_{n+1}=\tilde u_{n+1/2}+\frac{h}{2}\varphi(t_{n}+h/2,\tilde u_{n+1/2})\):

\( \text{(SEV)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2}\varphi(t_{n},u_{n}),\\ \hat u_{n+1}=\tilde u_{n+1/2}+\frac{h}{2}\varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2}),\\ u_{n+1}=u_n+\frac{h}{6} \left(\varphi(t_n,u_n)+4\varphi\left(t_n+\frac{h}{2},\tilde u_{n+1/2}\right)+\varphi(t_{n+1},\hat u_{n+1})\right)& n=0,1,2,\dots N-1 \end{cases} \)

2.1.4 Schémas de Simpson multipas

Au lieu de travailler sur \([t_n,t_{n+1}]\) on peut répéter la même construction sur \([t_{n-1},t_{n+1}]\) (et donc intégrer entre \(t_{n-1}\) et \(t_{n+1}\)).

Le SI devient alors

\( \text{(SMI)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{1}=u_0+h\varphi(t_{0},u_{0}),\\ u_{n+1}=u_{n-1}+\frac{h}{3} \left(\varphi(t_{n-1},u_{n-1})+4\varphi(t_n,u_{n})+\varphi(t_{n+1},u_{n+1})\right)& n=1,2,3,\dots N-1 \end{cases} \)

On peut remplacer le \(u_{n+1}\) dans le terme \(\varphi(t_{n+1},u_{n+1})\) par une approximation utilisant une prédiction d’Euler progressive à partir de \(t_n\), ce qui donne \(\tilde u_{n+1}=u_{n}+h\varphi(t_{n},u_{n})\):

\( \text{(SME)}\quad \begin{cases} u_0=y(t_0)=y_0,\\ u_{1}=u_0+h\varphi(t_{0},u_{0}),\\ \tilde u_{n+1}=u_{n}+h\varphi(t_{n},u_{n}),\\ u_{n+1}=u_{n-1}+\frac{h}{3} \left(\varphi(t_{n-1},u_{n-1})+4\varphi(t_n,u_{n})+\varphi(t_{n+1},\tilde u_{n+1})\right)& n=1,2,3,\dots N-1 \end{cases} \)

2.2 Calcul exact des polynomes interpolants et leur intégration

Pour vérifier nos calculs d’interpolation puis intégration, nous pouvons utiliser le package de calcul formel SymPy.

On notera

• t_n : \(t_n\)
• t_nm1 : \(t_{n-1}=t_n-h\)
• t_np1 : \(t_{n+1}=t_n+h\)
• t_np12 : \(t_{n+1/2}=t_n+\frac{h}{2}\)
• phi_n : \(\varphi(t_n,y(t_n))\)
• phi_nm1 : \(\varphi(t_{n-1},y(t_{n-1}))\)
• phi_np1 : \(\varphi(t_{n+1},y(t_{n+1}))\)
• phi_np12 : \(\varphi(t_{n+1/2},y(t_{n+1/2}))\)

import sympy as sym
from IPython.display import display, Math, Markdown
sym.init_printing()

phi_n    = sym.Symbol(r'\varphi_{n}')
phi_np12 = sym.Symbol(r'\varphi_{n+\frac</span><span class="sc">{1}{2}}')
phi_np1  = sym.Symbol(r'\varphi_{n+1}')
phi_nm1  = sym.Symbol(r'\varphi_{n-1}')
sym.var('h,t,t_n')

t_np1  = t_n+h
t_np12 = t_n+h/2
t_nm1  = t_n-h

display(Markdown("**EE**:"))
p  = sym.interpolate([(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
EE = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(EE)))

display(Markdown("**EI**:"))
p  = sym.interpolate([(t_np1,phi_np1)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
EI = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(EI)))

display(Markdown("**EM-init**:"))
p  = sym.interpolate([(t_np12,phi_np12)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
PM = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(PM)))

display(Markdown("**CN**:"))
p  = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
CN = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(CN)))

display(Markdown("**SI**:"))
p  = sym.interpolate([(t_np1,phi_np1),(t_np12,phi_np12),(t_n,phi_n)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
SI = sym.integrate(p,(t,t_n,t_np1)).simplify()
display(Math(r'\int_{t_n}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(SI)))

display(Markdown("**SMI**:"))
p  = sym.interpolate([(t_np1,phi_np1),(t_n,phi_n),(t_nm1,phi_nm1)], t)
display(Math(r't\mapsto \tilde f(t)='+sym.latex(p)))
SI = sym.integrate(p,(t,t_nm1,t_np1)).simplify()
display(Math(r'\int_{t_{n-1</span><span class="sc">}}^{t_{n+1</span><span class="sc">}}\tilde f(t) \mathrm{d}t='+sym.latex(SI)))

EE:

\(\displaystyle t\mapsto \tilde f(t)=\varphi_{n}\)

\(\displaystyle \int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\varphi_{n} h\)

EI:

\(\displaystyle t\mapsto \tilde f(t)=\varphi_{n+1}\)

\(\displaystyle \int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\varphi_{n+1} h\)

EM-init:

\(\displaystyle t\mapsto \tilde f(t)=\varphi_{n+\frac{1}{2}}\)

\(\displaystyle \int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\varphi_{n+\frac{1}{2}} h\)

CN:

\(\displaystyle t\mapsto \tilde f(t)=\frac{\varphi_{n+1} t}{h} - \frac{\varphi_{n+1} t_{n}}{h} + \varphi_{n} - \frac{\varphi_{n} t}{h} + \frac{\varphi_{n} t_{n}}{h}\)

\(\displaystyle \int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\frac{h \left(\varphi_{n+1} + \varphi_{n}\right)}{2}\)

SI:

\(\displaystyle t\mapsto \tilde f(t)=- \frac{\varphi_{n+1} t}{h} + \frac{\varphi_{n+1} t_{n}}{h} + \frac{2 \varphi_{n+1} t^{2}}{h^{2}} - \frac{4 \varphi_{n+1} t t_{n}}{h^{2}} + \frac{2 \varphi_{n+1} t_{n}^{2}}{h^{2}} + \frac{4 \varphi_{n+\frac{1}{2}} t}{h} - \frac{4 \varphi_{n+\frac{1}{2}} t_{n}}{h} - \frac{4 \varphi_{n+\frac{1}{2}} t^{2}}{h^{2}} + \frac{8 \varphi_{n+\frac{1}{2}} t t_{n}}{h^{2}} - \frac{4 \varphi_{n+\frac{1}{2}} t_{n}^{2}}{h^{2}} + \varphi_{n} - \frac{3 \varphi_{n} t}{h} + \frac{3 \varphi_{n} t_{n}}{h} + \frac{2 \varphi_{n} t^{2}}{h^{2}} - \frac{4 \varphi_{n} t t_{n}}{h^{2}} + \frac{2 \varphi_{n} t_{n}^{2}}{h^{2}}\)

\(\displaystyle \int_{t_n}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\frac{h \left(\varphi_{n+1} + 4 \varphi_{n+\frac{1}{2}} + \varphi_{n}\right)}{6}\)

SMI:

\(\displaystyle t\mapsto \tilde f(t)=\frac{\varphi_{n+1} t}{2 h} - \frac{\varphi_{n+1} t_{n}}{2 h} + \frac{\varphi_{n+1} t^{2}}{2 h^{2}} - \frac{\varphi_{n+1} t t_{n}}{h^{2}} + \frac{\varphi_{n+1} t_{n}^{2}}{2 h^{2}} - \frac{\varphi_{n-1} t}{2 h} + \frac{\varphi_{n-1} t_{n}}{2 h} + \frac{\varphi_{n-1} t^{2}}{2 h^{2}} - \frac{\varphi_{n-1} t t_{n}}{h^{2}} + \frac{\varphi_{n-1} t_{n}^{2}}{2 h^{2}} + \varphi_{n} - \frac{\varphi_{n} t^{2}}{h^{2}} + \frac{2 \varphi_{n} t t_{n}}{h^{2}} - \frac{\varphi_{n} t_{n}^{2}}{h^{2}}\)

\(\displaystyle \int_{t_{n-1}}^{t_{n+1}}\tilde f(t) \mathrm{d}t=\frac{h \left(\varphi_{n+1} + \varphi_{n-1} + 4 \varphi_{n}\right)}{3}\)