Soit le schéma de Runge-Kutta dont la matrice de Butcher est \(
\begin{array}{c|cc}
0 & 0 & 0\\
\alpha & \alpha & 0\\
\hline
& \frac{2\alpha-1}{2\alpha} & \frac{1}{2\alpha}
\end{array}
\)
Chaque sujet correspond à une valeur de \(\alpha\) choisie parmi les suivantes: \(
\left\{\frac{1}{6},\frac{1}{4},\frac{1}{3},\frac{5}{12},\frac{7}{12},\frac{2}{3},\frac{3}{4},\frac{5}{6},\frac{11}{12}\right\}
\)
Code
from fractions import Fraction
ALPHA=[Fraction(2,12),Fraction(3,12),Fraction(4,12),Fraction(5,12),
Fraction(7,12),Fraction(8,12),Fraction(9,12),Fraction(10,12),Fraction(11,12)]
for alpha in ALPHA:
print("c_2=a_{21}=",alpha,"\tb_1=",1-Fraction(1,2*alpha),"\tb_2=",Fraction(1,2*alpha))
c_2=a_{21}= 1/6 b_1= -2 b_2= 3
c_2=a_{21}= 1/4 b_1= -1 b_2= 2
c_2=a_{21}= 1/3 b_1= -1/2 b_2= 3/2
c_2=a_{21}= 5/12 b_1= -1/5 b_2= 6/5
c_2=a_{21}= 7/12 b_1= 1/7 b_2= 6/7
c_2=a_{21}= 2/3 b_1= 1/4 b_2= 3/4
c_2=a_{21}= 3/4 b_1= 1/3 b_2= 2/3
c_2=a_{21}= 5/6 b_1= 2/5 b_2= 3/5
c_2=a_{21}= 11/12 b_1= 5/11 b_2= 6/11
Exercice: écriture schéma
Écrire le schéma sous la forme d’une suite définie par récurrence.
Correction : écriture schéma
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\stackrel{\text{déf.}}{=} 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>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,... , t_{N}=T \}\) représente les points de la discrétisation.
- L’ensemble de \(N+1\) valeurs \(\{y_0, y_1,... , y_{N} \}\) représente la solution exacte.
- L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,... , u_{N} \}\) représente la solution numérique. Cette solution approchée sera obtenue en construisant une suite récurrente.
Le schéma qu’on va construire permet de calculer explicitement \(u_{n+1}\) à partir de \(u_n\) et il est donc possible de calculer successivement \(u_1\), \(u_2\),…, en partant de \(u_0\) par la formule de récurrence \(\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\alpha h,u_n+\alpha hK_1\right)\\
u_{n+1} = u_n + \frac{h}{2\alpha} \left((2\alpha-1)K_1+K_2\right) & n=0,1,\dots N-1
\end{cases}\)
Exercice : étude de la A-stabilité
Étudier théoriquement la A-stabilité.
Correction : étude de la A-stabilité
Définition de A-stabilité
Soit \(\beta>0\) un nombre réel positif et considérons le problème de Cauchy \(\begin{cases}
y'(t)=-\beta y(t), &\text{pour }t>0,\\
y(0)=y_0
\end{cases}\) où \(y_0\neq0\) est une valeur donnée. Sa solution est \(y(t)=y_0e^{-\beta t}\) donc \(\lim_{t\to+\infty}y(t)=0.\)
Si, sous d’éventuelles conditions sur \(h\), on a \(
\lim_{n\to+\infty} u_n =0,
\) alors on dit que le schéma est A-stable.
Le schéma donné appliqué à cette équation devient
Code
import sympy as sym
sym.init_printing()
sym.var('h,Alpha,beta,u_n,u_np1')
Phi=lambda y : -beta*y
K_1=Phi(u_n)
display(K_1)
K_2=Phi(u_n+Alpha*h*K_1)
display(K_2.expand().factor(u_n))
display(sym.Eq(u_np1,(u_n+h/(2*Alpha)*((2*Alpha-1)*K_1+K_2 )).factor()))
\(\displaystyle - \beta u_{n}\)
\(\displaystyle \beta u_{n} \left(A \beta h - 1\right)\)
\(\displaystyle u_{np1} = \frac{u_{n} \left(\beta^{2} h^{2} - 2 \beta h + 2\right)}{2}\)
c’est-à-dire \(\begin{cases}
u_0 = y_0 \\
u_{n+1} = \dfrac{2 - 2(\beta h) + (\beta h)^2}{2}u_n & n=0,1,\dots N-1
\end{cases}\) Par induction on obtient \(
u_{n}=\left(1 - (\beta h) + \frac{1}{2}(\beta h)^2 \right)^nu_0.
\) Notons que \(u_n\) ne dépend pas du choix de \(\alpha\).
Étant une suite géométrique, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \(
\left|1 - (\beta h) + \frac{1}{2}(\beta h)^2\right|<1.
\) Notons \(x\) le produit \(\beta h\) et \(q\) le polynôme \(q(x)=\frac{1}{2}x^2-x+1\).
Code
%matplotlib inline
import sympy as sym
from sympy.plotting import plot
sym.init_printing()
sym.var('x',positive=True)
q=1-1*x+x**2/2
print("|q(x)|<1 ssi", end=' ')
display(sym.solve(abs(q)<1))
print("q'(x)=")
dq=sym.diff(q,x)
display(dq)
print("q'(x)=0 ssi x=")
x_sommet=sym.solve(dq)[0]
display(x_sommet)
print("Comme q(x_sommet)=")
y_sommet=q.subs(x,x_sommet)
display(y_sommet)
print("donc q(x)>0 pour tout x et q(x)<1 ssi x<",x_sommet*2)
sym.plot(q,1,-1,(x,0,x_sommet*2));
\(\displaystyle \frac{1}{2}\)
donc q(x)>0 pour tout x et q(x)<1 ssi x< 2
Il s’agit de la parabole convexe de sommet \(\left(1,\frac{1}{2}\right)\).
Nous avons \(|q(x)|<1\) si et seulement si \(x<2\). La relation \(\lim\limits_{n\to+\infty}u_n=0\) est donc satisfaite si et seulement si \(
h <\frac{2}{\beta}.
\) Cette condition de stabilité limite le pas \(h\) d’avance en \(t\) lorsqu’on utilise le schéma.
Nous avons \(q(x)>0\) pour tout \(x\) donc cette convergence est monotone.
Exercice : implémentation et test
Implémenter le schéma et le tester avec le problème de Cauchy utilisé pour la A-stabilité sur l’intervalle \([0;1]\).
On choisira le nombre minimal de points pour avoir A-stabilité.
Chaque sujet correspond à une valeur de \(\beta\) choisie parmi les suivantes: \(
\left\{50,60,70,80,90,100\right\}
\)
Correction : implémentation et test
\(h<\frac{2}{\beta}\) sur l’intervalle \([0;1]\) impose \(N>\frac{\beta}{2}\).
Code
def RKalpha(phi,tt,y0):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i]+h*alpha , uu[i]+h*alpha*k1 )
uu.append( uu[i] + h*((2*alpha-1)*k1+k2 )/(2*alpha))
return uu
Code
from matplotlib.pylab import *
BETA=arange(50,101,10)
t0, y0, tfinal = 0, 1, 1
sol_exacte = lambda t : y0*exp(-beta*t)
phi = lambda t,y : -beta*y
i=0
for beta in BETA:
N=(int(beta/2)+1)
tt = linspace(t0,tfinal,N+1)
yy = [sol_exacte(t) for t in tt]
alpha=ALPHA[0]
uu = RKalpha(phi,tt,y0)
i+=1
figure(i,figsize=(20,5))
plot(tt,yy,lw=2,label=("Exacte"))
plot(tt,uu,'r*',label=("N="+str(N)+" beta="+str(beta)))
xlabel('$t$')
ylabel('$u$')
legend()
grid(True)
Exercice : étude de l’ordre
Étudier théoriquement l’ordre du schéma.
Correction : étude de l’ordre
- C’est un schéma explicite à deux étage: l’ordre est au plus 2.
- \(c_1=a_{11}+a_{12}\), \(c_2=a_{21}+a_{22}\) et \(1=b_1+b_2\) \(\implies\) il est consistante
- \(c_1 b_1 + c_2b_2=\frac{1}{2}\) \(\implies\) il est d’ordre \(2\)
Code
for alpha in ALPHA:
print("alpha=",alpha)
c=[0,alpha]
b=[1-1/(2*alpha),1/(2*alpha)]
A=[[0,0],[alpha,0]]
s=len(c)
ordre=[]
ordre.append(sum(b)==1)
for i in range(s):
ordre.append(sum(A[i])==c[i])
ordre.append(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print(set(ordre))
print("---------------------")
alpha= 1/6
{True}
---------------------
alpha= 1/4
{True}
---------------------
alpha= 1/3
{True}
---------------------
alpha= 5/12
{True}
---------------------
alpha= 7/12
{True}
---------------------
alpha= 2/3
{True}
---------------------
alpha= 3/4
{True}
---------------------
alpha= 5/6
{True}
---------------------
alpha= 11/12
{True}
---------------------
Exercice : solution exacte
Calculer la solution exacte du problème de Cauchy \(
\begin{cases}
y'(t)=-2ty^2(t), & t\in[0;3]\\
y(0)=2
\end{cases}
\)
Correction : solution exacte
On peut soit calculer la solution “à la main” soit utiliser sympy.
Il s’agit d’une EDO à variables séparables. La fonction \(y(t)=0\) pour tout \(t\) est solution de l’EDO mais elle ne vérifie pas la CI. Toute autre solution de l’EDO sera non nulle et se trouve formellement comme suit: \(
y'(t)=-2ty^2(t)
\quad\implies\quad
\frac{y'(t)}{y^2(t)}=-2t
\quad\implies\quad
\int y^{-2}\mathrm{d}y=-2\int t\mathrm{d}t
\quad\implies\quad
y(t)=\frac{1}{t^2+c}, \ c\in\mathbb{R}.
\) En imposant la CI on obtient \(2=1/C\) d’où l’unique solution du problème de Cauchy: \(y(t)=\frac{2}{2t^2+1}\).
Code
import sympy as sym
sym.init_printing()
sym.var('t,C1')
u=sym.Function('u')
f=-2*t*(u(t))**2
edo=sym.Eq(sym.diff(u(t),t),f)
print('EDO')
display(edo)
solgen=sym.dsolve(edo,u(t))
print('Solution générale:')
display(solgen)
t0=0
u0=2
consts = sym.solve( [solgen.rhs.subs(t,t0)-u0 ], dict=True)[0]
print("En prénant en compte la CI on fixe la constante d'intégation:")
display(consts)
solpar=solgen.subs(consts)
print('Solution particulière:')
display(solpar)
print('Graphe:')
%matplotlib inline
sym.plot(solpar.rhs,(t,0,3));
\(\displaystyle \frac{d}{d t} u{\left(t \right)} = - 2 t u^{2}{\left(t \right)}\)
\(\displaystyle u{\left(t \right)} = \frac{1}{C_{1} + t^{2}}\)
En prénant en compte la CI on fixe la constante d'intégation:
\(\displaystyle \left\{ C_{1} : \frac{1}{2}\right\}\)
\(\displaystyle u{\left(t \right)} = \frac{1}{t^{2} + \frac{1}{2}}\)
Exercice : convergence
Vérifier empiriquement l’ordre de convergence en affichant la courbe de convergence et en estimant la pente sur ce problème de Cauchy.
Correction : étude empirique de l’ordre de convergence
On calcule la solution approchée avec différentes valeurs de \(h_k=1/N_k\), à savoir \(N_k=2^3\), … \(2^{10}\). On sauvegarde les valeurs de \(h_k\) dans le vecteur H.
Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on sauvegarde toutes ces erreurs dans le vecteur err de sort que err[k] contient \(e_k=\max_{i=0,...,N_k}|y(t_i)-u_{i}|\) avec \(N_k=2^{k+1}\).
Pour afficher l’ordre de convergence on utilise une échelle logarithmique : on représente \(\ln(h)\) sur l’axe des abscisses et \(\ln(\text{err})\) sur l’axe des ordonnées ainsi, si \(\text{err}=Ch^p\), alors \(\ln(\text{err})=\ln(C)+p\ln(h)\). En échelle logarithmique, \(p\) représente donc la pente de la ligne droite \(\ln(\text{err})\).
Pour estimer l’ordre de convergence on estime la pente de la droite qui relie l’erreur au pas \(k\) à l’erreur au pas \(k+1\) en echelle logarithmique en utilisant la fonction polyfit basée sur la régression linéaire.
Code
t0, y0, tfinal = 0, 2, 3
sol_exacte = lambda t : 1/(t**2+1/y0)
phi = lambda t,y : -2*t*y**2
i=0
for alpha in ALPHA:
H = []
err = []
for k in range(8):
N = 2**(k+3)
tt = linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]
H.append(h)
yy = [sol_exacte(t) for t in tt]
uu = RKalpha(phi,tt,y0)
err.append(max([abs(uu[i]-yy[i]) for i in range(len(uu))]))
i+=1
figure(i,figsize=(20,5))
loglog(H,err, 'r-o')
title('Alpha={}, Ordre = {:1.4f}'.format(alpha,polyfit(log(H),log(err), 1)[0]))
xlabel('$h$')
ylabel('$e$')
#legend()
grid(True)
