2023 Contrôle Continu

⏱ 2h30, deposer le notebook dans Moodle

1 EXERCICE 1 [4 points]: choix d’un schéma

Soit le problème de Cauchy

\( \begin{cases} y'(t)=\left(1-y(t)-\dfrac{1}{4+t}\right)\sin(2\pi y(t)), &t\in[0;10]\\ y(0)=y_0 \ \in\ ]0.7;0.9[ \end{cases} \)

On ne peut pas calculer analytiquement la solution exacte.

Pour calculer des solutions approchées on doit utiliser un schéma numérique. Quelle caractéristique doit-il avoir ? Justifier la réponse.

Q [4 points] Pour calculer des solutions approchées on doit utiliser un schéma numérique. Quelle caractéristique doit-il avoir ? Justifier la réponse.

🐸 Il s’agit d’un des problèmes de Cauchy proposés lors du CT 2023 de l’ECUE M61 Équations différentielles.

Bien qu’on ne puisse pas calculer analytiquement la solution exacte, nous pouvons tracer les lignes de courant pour avoir une idée du comportement de la solution en fonction de la donnée initiale. On en déduit que le problème de Cauchy est numériquement mal posé. Il faut alors utiliser un schéma d’ordre elevé.

Pour s’en convaincre nous pouvons aussi demander à odeint du module scipy de tracer la solution approchée pour différentes valeurs de \(y_0\).

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# variables globales
t0     = 0
tfinal = 10

phi = lambda t,y: np.sin(np.pi*2*y)*(1-1/(4+t)-y)

N = 200
tt = np.linspace(t0, tfinal, N + 1)

plt.figure(figsize=(10,5))

# On trace d'abord les lignes de courant
g1  = np.linspace(0,10,21)
g2  = np.linspace(0.5,1,21)
T,Y = np.meshgrid(g1,g2)  
DT, DY = 1, phi(T,Y)          # compute growth rate on the grid
M = np.hypot(DT,DY) # = sqrt(DT**2+DY**2)         # norm growth rate 
plt.streamplot(T,Y, DT/M, DY/M, color=M, linewidth=0.5, cmap=plt.cm.viridis, density=2, arrowstyle='->', arrowsize=1.5)
plt.grid()
plt.xlabel('t')
plt.ylabel('y')
plt.title('Lignes de courant');

# On trace quelque solution approchée
for y0 in np.linspace(0.7,0.9,7):
    uu = odeint(phi,y0,tt, tfirst=True)
    plt.plot(tt,uu,lw=4,label=f'$y_0={</span>y0<span class="sc">:g}$')

plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left');

2 Exercice 2 [7 points] : étude d’un schéma multipas

On notera \(\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)\). Considérons la méthode multipas

\( u_{n+1}=u_{n-1}+h\left( \varphi_{n+1}+\varphi_{n-1}\right) \)

1. Étudier la zéro-stabilité de cette méthode.
2. La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l’ordre de convergence.
3. Vérifier empiriquement l’ordre de convergence sur le problème de Cauchy
   \( \begin{cases} y'(t) = (y(t)+2)\dfrac{2t^2+2t-1}{1+t}, &\forall t \in I=[0,2],\\ y(0) = -1, \end{cases} \)

On remarque qu’il s’agit de la méthode de Crank-Nicolson sur l’intervalle \([t_{n-1},t_{n+1}]\).

Q2 [1 points] Étudier la zéro-stabilité.

Le schéma est de la forme

\( u_{n+1} = a_0u_n+a_1u_{n-1}+hb_0\varphi_n+hb_1\varphi_{n-1}+hb_{-1}\varphi_{n+1}. \)

On a

• \(p=1\), i.e. c’est un méthode à \(q=p+1=2\) pas;
• \(a_0=0\) et \(a_1=1\),
• \(b_0=0\) et \(b_1=1\),
• \(b_{-1}=1\neq0\) (\(\implies\) la méthode est implicite).

Le premier polynôme caractéristique est

\( \varrho(r)= %r^{p+1}-\sum_{j=0}^p a_jr^{p-j}= r^2-a_0r-a_1= r^2-1 =(r-1)(r+1). \)

Les racines sont

\( r_0=1,\quad r_{1}=-1. \)

La méthode est donc zéro-stable car

\( \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \text{ multiplicité de $r_j$ > }1 \implies |r_j|<1 \end{cases} \)

Q2 [3 points] La méthode est elle consistente? Est elle convergente? Si oui, étudier théoriquement l’ordre de convergence.

La méthode étant zéro-stable, si elle est consistante, alors elle est convergente. Pour étudier la convergence introduisons la quantité
\( \begin{align*} \xi(i) &= \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j \\ &= (-0)^{i}a_0+(-1)^ia_1+i \Big( (1)^{i-1}b_{-1}+(-0)^{i-1}b_0+(-1)^{i-1}b_1 \Big) \end{align*} \)
avec \((-0)^0=1\).

Pour que la méthode soit consistante il faut que \(\xi(0)=\xi(1)=1\). De plus, pour que la méthode soit d’ordre \(\omega\), il faut que \(\xi(i)=1\) pour tout \(i\le\omega\).

De plus, \(\omega\) vérifie la première barrière de Dahlquist: le schéma étant implicite à \(q=2\) (pair) pas consistante et zéro-stable, \(\omega=q+2\le4\).

En calculant \(\xi\) on trouve que la méthode est d’ordre 2 (cf. calculs avec sympy ci-dessous).

%reset -f
import sympy 
 
aa_eval = [0,1]
bb_eval = [0,1]
bm1_eval = 1 


###########################################################################

from IPython.display import display, Math, Markdown

def xi(i,aa,bb,bm1):
    sa = sum( [ (-j)**i * aa[j] for j in range(q) ] )
    sb = bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,q)] )
    if i==1:
        sb += bb[0]
    return (sa)+(i*sb)

q = len(aa_eval) # nb de pas
p = q-1  # j = 0...p # coeffs du schema
omega = q+2 if q%2==0 else q+1

aa = sympy.symbols(f'a_0:{</span>q<span class="sc">}')
bb = sympy.symbols(f'b_0:{</span>q<span class="sc">}')
bm1 = sympy.Symbol('b_{-1}')

display(Markdown(f"C'est une méthode à $q={</span>p <span class="op">+</span> <span class="dv">1</span><span class="sc">}$ pas d'ordre $\\omega\\le{</span>omega<span class="sc">}$ avec"))

texte = ''.join([ f"{</span>sympy<span class="sc">.</span>latex(aa[j])<span class="sc">}={</span>sympy<span class="sc">.</span>latex(aa_eval[j])<span class="sc">},"+r"\quad " for j in range(p+1)])
texte += ''.join([ f"{</span>sympy<span class="sc">.</span>latex(bb[j])<span class="sc">}={</span>sympy<span class="sc">.</span>latex(bb_eval[j])<span class="sc">},"+r"\quad " for j in range(p+1)])
texte += f"{</span>sympy<span class="sc">.</span>latex(bm1)<span class="sc">}={</span>sympy<span class="sc">.</span>latex(bm1_eval)<span class="sc">}"
display(Math( texte ))  

        
xi_eval = [xi(i,aa_eval,bb_eval,bm1_eval) for i in range(omega+1)]
for i in range(omega+1):
    display(Math( f"\\xi({</span>i<span class="sc">}) = {</span>sympy<span class="sc">.</span>latex(xi(i,aa,bb,bm1))<span class="sc">} = {</span>sympy<span class="sc">.</span>latex(xi_eval[i])<span class="sc">}" ))
    if xi_eval[i]!=1:
        display(Markdown( f"La méthode est donc d'ordre $\\omega = {</span>i<span class="op">-</span><span class="dv">1</span><span class="sc">}$" ))   
        break

C’est une méthode à \(q=2\) pas d’ordre \(\omega\le4\) avec

\(\displaystyle a_{0}=0,\quad a_{1}=1,\quad b_{0}=0,\quad b_{1}=1,\quad b_{-1}=1\)

\(\displaystyle \xi(0) = a_{0} + a_{1} = 1.0\)

\(\displaystyle \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1\)

\(\displaystyle \xi(2) = a_{1} - 2 b_{1} + 2 b_{-1} = 1\)

\(\displaystyle \xi(3) = - a_{1} + 3 b_{1} + 3 b_{-1} = 5\)

La méthode est donc d’ordre \(\omega = 2\)

Q3 [3 points] Vérifier empiriquement l’ordre de convergence sur le problème de Cauchy \(\begin{cases} y'(t) = (y(t)+2)\dfrac{2t^2+2t-1}{1+t}, &\forall t \in I=[0,2],\\ y(0) = -1, \end{cases} \)

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


# variables globales
t0     = 0
tfinal = 2
y0     = -1

phi = lambda t,y: (y+2)*(2*t**2+2*t-1)/(1+t) 

##############################################
# solution exacte
##############################################
t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgen = sympy.dsolve(edo)
# display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
# display(consts)
solpar = solgen.subs(consts)
display(solpar)

sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

##############################################
# schéma (initialisation avec sol exacte pour ordre de convergence)
##############################################
def multipas(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0:2] = sol_exacte(tt[0:2])
    for i in range(1,len(tt) - 1):
        eq = lambda x : -x+uu[i-1]+h*(phi(tt[i+1],x)+phi(tt[i-1],uu[i-1]))
        temp = fsolve( eq ,uu[i])
        uu[i+1] = temp[0]
    return uu

##############################################
# ordre
##############################################
H = []
err = []
N = 10

plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

for k in range(6):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = multipas(phi, tt, sol_exacte)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()


ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Multipas ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

\(\displaystyle y{\left(t \right)} = \frac{- 2 t + e^{t^{2}} - 2}{t + 1}\)

3 Exercice 3 [10 points] : étude d’un schéma RK

Considérons le schéma dont la matrice de Butcher est donnée ci-dessous:

\( \begin{array}{c|cc} \sigma & \sigma & 0 \\ 1 & \vartheta & 1-\vartheta\\ \hline & \gamma & 1-\gamma \end{array} \)

1. Déterminer les coefficients \(\sigma\), \(\vartheta\), \(\gamma\) pour que le schéma soit d’ordre 3. Est-il d’ordre 4 ?
2. Étudier empiriquement l’ordre de convergence sur le problème de Cauchy
   \( \begin{cases} y'(t)=\dfrac{t^2+y^2(t)}{ty(t)}, & t\in[1;5]\\ y(1)=1 \end{cases} \)
3. Étudier théoriquement la A-stabilité du schéma optimal.

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()


from IPython.display import display, Math, Markdown
prlat = lambda *args : display(Math(''.join( sympy.latex(arg) for arg in args )))

Q1 [5 points] Pour quelles valeurs des paramètres le schéma est d’ordre 3? Et d’ordre 4?

Soit \(\omega\) l’ordre de la méthode.

C’est un schéma implicite à \(2\) étages (\(s=2\)) donc \(\omega\le2s=4\).

sigma = sympy.Symbol(r"\sigma")
theta = sympy.Symbol(r"\vartheta")
gamma = sympy.Symbol(r"\gamma")


def ordre_RK(s, A=None, b=None, c=None):

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    matrice_Butcher(s, A, b, c)
    
    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    ordre_1(s, A, b, c, j)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    ordre_2(s, A, b, c, j)

    display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
    ordre_3(s, A, b, c, j)    
    
    display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
    ordre_4(s, A, b, c, j)

    return None


def matrice_Butcher(s, A, b, c):
    But = sympy.Matrix(A)
    But = But.col_insert(0, sympy.Matrix(c))   
    last = [sympy.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sympy.Matrix(last).T)
    display(Math(sympy.latex(sympy.Matrix(But))))
    return None


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
    return None    


def ordre_2(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
    texte += r"\text{ doit être égale à }\frac{1}{2}"
    display(Math(texte))
    return None


def ordre_3(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j^2='
    texte += sympy.latex( sum([b[i]*c[i]**2 for i in range(s)]).simplify() )
    texte += r"\text{ doit être égale à }\frac{1}{3}"
    display(Math(texte))
    texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j='    
    somma = sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{6}"
    display(Math(texte))
    return None


def ordre_4(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j^3='
    texte += sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).simplify() )
    texte += r"\text{ doit être égale à }\frac{1}{4}"
    display(Math(texte))
    texte = r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
    somma = sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{8}"
    display(Math(texte))
    texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='
    somma = sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{12}"
    display(Math(texte))
    texte = r'\sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k='
    somma = sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{24}"
    display(Math(texte))
    return None


# APPLICATION A NOTRE CAS
c = [sigma,sympy.S(1)]
b = [gamma,1-gamma]
A = [[sigma,sympy.S(0)],[theta,1-theta]]
s = len(c)
display(Markdown(f"La méthode de Runge-Kutta est à {</span>s<span class="sc">} étages"))
ordre_RK(s,A,b,c)

La méthode de Runge-Kutta est à 2 étages

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}\sigma & \sigma & 0\\1 & \vartheta & 1 - \vartheta\\ & \gamma & 1 - \gamma\end{matrix}\right]\)

On a 3 conditions pour avoir consistance = pour être d’ordre 1

\(\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1\)

\(\displaystyle \sum_{j=1}^s a_{0j}=\sigma\text{ doit être égale à }\sigma\)

\(\displaystyle \sum_{j=1}^s a_{1j}=1\text{ doit être égale à }1\)

On doit ajouter 1 condition pour être d’ordre 2

\(\displaystyle \sum_{j=1}^s b_j c_j=\gamma \sigma - \gamma + 1\text{ doit être égale à }\frac{1}{2}\)

On doit ajouter 2 conditions pour être d’ordre 3

\(\displaystyle \sum_{j=1}^s b_j c_j^2=\gamma \sigma^{2} - \gamma + 1\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\gamma \sigma^{2} - \sigma \vartheta \left(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\text{ doit être égale à }\frac{1}{6}\)

On doit ajouter 4 conditions pour être d’ordre 4

\(\displaystyle \sum_{j=1}^s b_j c_j^3=\gamma \sigma^{3} - \gamma + 1\text{ doit être égale à }\frac{1}{4}\)

\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\gamma \sigma^{3} - \sigma \vartheta \left(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\text{ doit être égale à }\frac{1}{8}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\gamma \sigma^{3} - \sigma^{2} \vartheta \left(\gamma - 1\right) + \left(\gamma - 1\right) \left(\vartheta - 1\right)\text{ doit être égale à }\frac{1}{12}\)

\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=\gamma \sigma^{3} - \sigma^{2} \vartheta \left(\gamma - 1\right) + \sigma \vartheta \left(\gamma - 1\right) \left(\vartheta - 1\right) + \left(1 - \gamma\right) \left(\vartheta - 1\right)^{2}\text{ doit être égale à }\frac{1}{24}\)

On commence à resoudre le système dont les inconnues sont les paramètres \(\sigma, \gamma, \vartheta\) pour avoir au moins l’ordre 3:

cnst = sympy.solve( [ sum([b[i]*c[i] for i in range(s)])-1/sympy.S(2),
                    sum([b[i]*c[i]**2 for i in range(s)])-1/sympy.S(3),
                    sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])-1/sympy.S(6)
        ] )[0]

display(Math(sympy.latex(cnst)))


c_substituted = [expr.subs(cnst) for expr in c]
b_substituted = [expr.subs(cnst) for expr in b]
A_substituted = [[expr.subs(cnst) for expr in row] for row in A]

\(\displaystyle \left\{ \gamma : \frac{3}{4}, \ \sigma : \frac{1}{3}, \ \vartheta : 1\right\}\)

Remarque : le choix optimal conduit au même schéma du DM de 2020, l’étude a donc déjà été fait en TD et en TP !

Le schéma n’est pas d’ordre 4 :

ordre_RK(s,A_substituted,b_substituted,c_substituted)

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}\frac{1}{3} & \frac{1}{3} & 0\\1 & 1 & 0\\ & \frac{3}{4} & \frac{1}{4}\end{matrix}\right]\)

On a 3 conditions pour avoir consistance = pour être d’ordre 1

\(\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1\)

\(\displaystyle \sum_{j=1}^s a_{0j}=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{j=1}^s a_{1j}=1\text{ doit être égale à }1\)

On doit ajouter 1 condition pour être d’ordre 2

\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}\text{ doit être égale à }\frac{1}{2}\)

On doit ajouter 2 conditions pour être d’ordre 3

\(\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}\text{ doit être égale à }\frac{1}{6}\)

On doit ajouter 4 conditions pour être d’ordre 4

\(\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{5}{18}\text{ doit être égale à }\frac{1}{4}\)

\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{1}{9}\text{ doit être égale à }\frac{1}{8}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{1}{18}\text{ doit être égale à }\frac{1}{12}\)

\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=\frac{1}{18}\text{ doit être égale à }\frac{1}{24}\)

Q2 [3 points] Étudier empiriquement l’ordre de convergence sur le problème de Cauchy \( \begin{cases} y'(t)=\dfrac{t^2+y^2(t)}{ty(t)}, & t\in[1;5]\\ y(1)=1 \end{cases} \)

# variables globales
t0     = 1
tfinal = 5
y0     = 1

phi = lambda t,y: (t**2+y**2)/(t*y)

##############################################
# solution exacte
##############################################
import sympy 
sympy.init_printing()

t   = sympy.Symbol('t')
y   = sympy.Function('y')
edo = sympy.Eq( sympy.diff(y(t),t) , phi(t,y(t)) )
solgenLIST = sympy.dsolve(edo,y(t))
display(solgenLIST)
for solgen in solgenLIST: 
    #display(solgen)
    consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)
    #display(consts)
    if consts!=[]:
        solpar = solgen.subs(consts[0]).simplify()
display(solgen)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar = solgen.subs(consts).simplify()
display(solpar)


##############################################
# solution approchée
##############################################
sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

# Je note sigma_e etc pour sigma evaluée comme un float
sigma_e = 1/3
gamma_e = 1/(2*(1-sigma_e)) # pour ordre 2
theta_e = 1  

def RK(tt,phi,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        sys = lambda X : [ -X[0]+phi( tt[i]+h*sigma_e , uu[i]+h*X[0]*sigma_e  ) , 
                           -X[1]+phi( tt[i]+h*theta_e , uu[i]+h*(X[0]*theta_e+X[1]*(1-theta_e))  )
                           ]
        K_start = phi(tt[i],uu[i])
        K1, K2 = fsolve( sys , [K_start,K_start ] )
        uu.append( uu[i]+h*(gamma_e*K1+(1-gamma_e)*K2) )
    return uu


##############################################
# ordre
##############################################
H = []
err = []

figure(figsize=(16,5))
ax1 = subplot(1,2,1)
ax2 = subplot(1,2,2)

N = 50
for k in range(6):
    N += 50
    tt = linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(tt,phi,y0)
    err.append( norm(uu-yy,inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte')
xlabel('$t$')
ylabel('$y$')
ax1.grid(True)
ax1.legend()

ax2.plot( log(H), log(err), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title(f'RK : ordre = {</span>polyfit(log(H),log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

\(\displaystyle \left[ y{\left(t \right)} = - t \sqrt{C_{1} + 2 \log{\left(t \right)}}, \ y{\left(t \right)} = t \sqrt{C_{1} + 2 \log{\left(t \right)}}\right]\)

\(\displaystyle y{\left(t \right)} = t \sqrt{C_{1} + 2 \log{\left(t \right)}}\)

\(\displaystyle \left\{ C_{1} : 1\right\}\)

\(\displaystyle y{\left(t \right)} = t \sqrt{2 \log{\left(t \right)} + 1}\)

Q3 [2 points] Étudier théoriquement la A-stabilité du schéma optimal.

La méthode appliquée à l’EDO \(y'(t)=-\beta y(t)\) avec \(\beta>0\) s’écrit: \(\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-(\beta h) \sigma K_1 \leadsto K_1=\frac{-\beta}{1+(\beta h) \sigma}u_n\\ K_2 = -\beta\left(u_n+h K_1 \vartheta + h K_2 (1-\vartheta)\right)\leadsto K_2=\frac{1+\frac{(\beta h)\vartheta}{1+(\beta h) \sigma}}{1+(\beta h)(1-\vartheta)}(-\beta u_n)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} \)

En particulier, dans le cas d’ordre 3, on a \(\begin{cases} u_0 = y_0 \\ K_1 = -\beta u_n-\beta \frac{h}{3}K_1 \leadsto K_1=\frac{-3\beta}{3+\beta h}u_n\\ K_2 = -\beta\left(u_n+hK_1\right)\leadsto K_2=\frac{2\beta^2h-3\beta}{3+\beta h}u_n\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases} \) L’étude ci-dessous montre que, dans ce cas, le schéma est A-stable ssi \(h<\frac{6}{\beta}\).

sympy.var('u_n, h, β, K1, K2,')  

g = lambda y : -β*y
eq1 = sympy.Eq( K1 , g(u_n+h*K1*sigma))
eq2 = sympy.Eq( K2 , g(u_n+h*(K1*theta+K2*(1-theta))) )
sol = sympy.solve([eq1,eq2],[K1,K2])
sol = { key:value.factor() <span class="cf">for</span> key,value <span class="kw">in</span> sol.items() }
display(sol)

RHS = u_n+h*(gamma*K1+(1-gamma)*K2).factor()
RHS = RHS.subs(sol).simplify().factor()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte)) 

print(r"Notons x=hβ>0. On trouve donc")
x = sympy.Symbol('x',positive=True)
RHS = RHS.subs(h*β,x).simplify()
texte = 'u_{n+1}=' + sympy.latex(RHS)
display(Math(texte)) 

print("Considérons le schéma d'ordre 3:")

display(cnst)

RHS = RHS.subs(cnst)
texte = r"\text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1</span><span class="sc">}}{u_n}="
texte += sympy.latex(RHS/u_n)
display(Math(texte)) 
sympy.solve([-1<RHS/u_n,RHS/u_n<1])

\(\displaystyle \left\{ K_{1} : - \frac{u_{n} β}{\sigma h β + 1}, \ K_{2} : \frac{u_{n} β \left(\sigma h β - \vartheta h β + 1\right)}{\left(\sigma h β + 1\right) \left(\vartheta h β - h β - 1\right)}\right\}\)

\(\displaystyle u_{n+1}=- \frac{u_{n} \left(\gamma \sigma h^{2} β^{2} - \gamma h^{2} β^{2} - \sigma \vartheta h^{2} β^{2} + \sigma h β + \vartheta h^{2} β^{2} - \vartheta h β + 1\right)}{\left(\sigma h β + 1\right) \left(\vartheta h β - h β - 1\right)}\)

Notons x=hβ>0. On trouve donc
Considérons le schéma d'ordre 3:

\(\displaystyle u_{n+1}=\frac{u_{n} \left(\gamma \sigma x^{2} - \gamma x^{2} - \sigma \vartheta x^{2} + \sigma x + \vartheta x^{2} - \vartheta x + 1\right)}{\left(\sigma x + 1\right) \left(- \vartheta x + x + 1\right)}\)

\(\displaystyle \left\{ \gamma : \frac{3}{4}, \ \sigma : \frac{1}{3}, \ \vartheta : 1\right\}\)

\(\displaystyle \text{La méthode est A-stable ssi }|q(x)|<1\text{ avec }q(x)=\frac{u_{n+1}}{u_n}=\frac{\frac{x^{2}}{6} - \frac{2 x}{3} + 1}{\frac{x}{3} + 1}\)

\(\displaystyle x < 6\)