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\} \)
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Écrire le schéma sous la forme d’une suite définie par récurrence.
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)\).
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}\)
Étudier théoriquement 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
\(\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\).
%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));|q(x)|<1 ssi q'(x)=
q'(x)=0 ssi x=
Comme q(x_sommet)=
donc q(x)>0 pour tout x et q(x)<1 ssi x< 2\(\displaystyle x < 2\)
\(\displaystyle x - 1\)
\(\displaystyle 1\)
\(\displaystyle \frac{1}{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.
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\} \)
\(h<\frac{2}{\beta}\) sur l’intervalle \([0;1]\) impose \(N>\frac{\beta}{2}\).
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)
Étudier théoriquement l’ordre du schéma.
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}
---------------------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} \)
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));EDO
Solution générale:
En prénant en compte la CI on fixe la constante d'intégation:
Solution particulière:
Graphe:\(\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}}\)
\(\displaystyle \left\{ C_{1} : \frac{1}{2}\right\}\)
\(\displaystyle u{\left(t \right)} = \frac{1}{t^{2} + \frac{1}{2}}\)
Vérifier empiriquement l’ordre de convergence en affichant la courbe de convergence et en estimant la pente sur ce problème de Cauchy.
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)
