%reset -f
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Contrôle Continu 2019¶
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\} $$
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))
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
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()))
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$.
%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));
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}$.
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
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$
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("---------------------")
Exercice : solution exacte¶
Calculer la solution exacte du problème de Cauchy $$ \begin{cases} y'(t)=-2t\left(y(t)\right)^2, & 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}$.
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));
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.
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)
