In [ ]:
from IPython.core.display import HTML, display
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[ ]:

2022 CC corrigé

  • ⏱ 2H30
  • ✍ Deposer votre fichier .ipynb sur Moodle (dépôt de devoir)
In [ ]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

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

Exercice 1 : étude d'un schéma multipas

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $$ u_{n+1} = \alpha u_n + (1-\alpha) u_{n-1} +h\left(2\varphi_{n+1} -\frac{3\vartheta}{2}\varphi_{n}+\frac{\vartheta}{2} \varphi_{n-1}\right). $$

Q1 [2 points]
Pour quelles valeurs des paramètres $\alpha$ et $\vartheta$ la méthode est zéro-stable?

Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-\alpha r+(\alpha-1)=(r-1)\big(r-(\alpha-1)\big) $$ dont les racines sont $$ r_0=1,\quad r_1=\alpha-1. $$ La méthode est donc zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,1 \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1 \end{cases} $$ donc ssi $0\le\alpha<2$.

Q2 [3 points]
Quel est l'ordre de convergence de la méthode?

On rappelle qu'une méthode est convergente ssi elle est zéro-stable et consistante. Nous allons donc étudier la consistance (puis l'ordre de consistance = ordre de convergence) pour les seules valeurs de $\alpha$ pour lesquelles on a la zéro-stabilité.

  • On a $p=1$: c'est une méthode à $q=p+1=2$ pas.
  • La méthode est implicite car $b_{-1}=2\neq0$.

Pour que la méthode soit consistante il faut que $$\begin{cases} \displaystyle\sum_{j=0}^p a_j=a_0+a_1=1, \\ \displaystyle\sum_{j=0}^p (-j)a_j+\sum_{j=-1}^p b_j=1. \end{cases}$$

De plus, la première barrière de Dahlquist affirme qu'un schéma implicite à $q=2$ pas consistante et zéro-stable peut être au mieux d'ordre $\omega=q+2=4$.

Pour que la méthode soit

  • au moins d'ordre 2 il faut qu'elle soit consistante et que $\displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=1$,
  • au moins d'ordre 3 il faut qu'elle soit au moins d'ordre 2 et que $\displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=1$,
  • au moins d'ordre 3 il faut qu'elle soit au moins d'ordre 2 et que $\displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=1$.
In [ ]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))
In [ ]:
# j = 0...p
# q = p+1 # nb de pas
# ω # ordre de la méthode

CAS_GENERAL = True
p = 1

q = p+1 
ordre_max = q+2 if q%2==0 else q+1

print(f"{'★'*70}\nC'est une méthode à q = {q} pas d'ordre ω ≤ {ordre_max}")


if CAS_GENERAL:
    aa = sympy.symbols(f'a_0:{q}')
    bb = sympy.symbols(f'b_0:{q}')
    bm1 = sympy.Symbol('b_{-1}')
else : # cas particulier
    α = sympy.Symbol('α')
    ϑ = sympy.Symbol('ϑ')
    aa = [α, 1-α]
    bb = [-3*ϑ/2, ϑ/2]
    bm1 = 2

i=1
sa=sum( [-j*aa[j] for j in range(len(aa))] )
sb=bm1+sum( [bb[j] for j in range(len(bb))] )
print(f"{'★'*70}\nLa méthode est d'ordre ω ≥ {i} (= consistante) si ")
prlat(sum(aa).factor(),"=1" )
prlat((sa).factor()+(i*sb).factor(),"=1")

for i in range(2,ordre_max+1):
    sa=sum( [(-j)**i*aa[j] for j in range(len(aa))] )
    sb=bm1+sum( [(-j)**(i-1)*bb[j] for j in range(1,len(bb))] )
    print(f"{'★'*70}\nLa méthode est d'ordre ω ≥ {i} si elle est d'ordre {i-1} et, de plus, ")
    prlat((sa).factor()+(i*sb).factor(),"=1" )

★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
C'est une méthode à q = 2 pas d'ordre ω ≤ 4
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 1 (= consistante) si
$\displaystyle a_{0} + a_{1}=1$
$\displaystyle - a_{1} + b_{0} + b_{1} + b_{-1}=1$ 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 2 si elle est d'ordre 1 et, de plus,
$\displaystyle a_{1} - 2 \left(b_{1} - b_{-1}\right)=1$ 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 3 si elle est d'ordre 2 et, de plus,
$\displaystyle - a_{1} + 3 \left(b_{1} + b_{-1}\right)=1$ 
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★
La méthode est d'ordre ω ≥ 4 si elle est d'ordre 3 et, de plus,
$\displaystyle a_{1} - 4 \left(b_{1} - b_{-1}\right)=1$

Dans notre cas

  • la méthode est consistante ssi $\vartheta=\alpha$ car $a_0+a_1=\alpha+(1-\alpha)=1$ et $-a_1+b_0+b_1+b_{-1}=-(1-\alpha)-\frac{3\vartheta}{2}+\frac{\vartheta}{2}+2=\alpha-\vartheta+1$;
  • pour que la méthode sois d'ordre au moins 2 il faudrait avoir $1= a_1-2(b_1-b_{-1})=1-\alpha-2(\frac{\alpha}{2}-2)=5-2\alpha$ ce qui n'est pas compatible avec la condition de zéro-stabilité.

Conclusion: si $0\le\alpha=\vartheta<1$ le schéma est convergent à l'ordre 1, dans les autres cas il n'est pas convergente.

Q3 [3 points]
Nous allons maintenant fixer $\alpha$. Pour cela, écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour le paramètre $\alpha$. Pour $\vartheta$, choisir une valeur qui garantie la convergence.

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

In [ ]:
nom = "Faccanoni"
prenom = "Gloria"

NUM = list(range(1,16))
idx = sum([ord(c) for c in nom+prenom])%len(NUM)
my_alpha = NUM[idx]/sympy.S(8)
prlat(r'\alpha=',my_alpha)
my_alpha = float(my_alpha) # conversion en float
$\displaystyle \alpha=\frac{5}{8}$
  • Pour estimer les erreurs on définit la solution exacte (calculée en utilisant le module sympy ou à la main).
  • On calcule la solution approchée pour différentes valeurs de $N$. Pour le calcul de l'ordre de convergence on initialise la suite avec la solution exacte.
In [ ]:
# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y: t**2-y

##############################################
# 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)
consts = sympy.solve( sympy.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
solpar = solgen.subs(consts).simplify()
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 = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,len(tt) - 1):
        eq = lambda x : -x+my_alpha*uu[i]+(1-my_alpha)*uu[i-1]+h*(2*phi(tt[i+1],x)-3*my_alpha/2*phi(tt[i],uu[i])+my_alpha/2*phi(tt[i-1],uu[i-1]))
        temp = fsolve( eq ,uu[i])
        uu.append( temp[0] )
    return uu

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

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

for k in range(4):
    N += 20
    tt = linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = multipas(phi, tt, sol_exacte)
    err.append( norm(uu-yy,inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={N}')

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

ax2 = subplot(1,2,2)
ax2.loglog( H, err, 'r-o')
xlabel('$h$')
ylabel('$e$')
title(f'Multipas ordre = {polyfit(log(H),log(err),1)[0]:1.2f}')
ax2.grid(True)
$\displaystyle y{\left(t \right)} = t^{2} - 2 t + 2$

Exercice 2 : étude d'un système

Soit le problème de Cauchy $$\begin{cases} y'(t)=\eta z(t),\\ z'(t)=y(t)+\eta z(t),\\ y(0)=1,\\ z(0)=1. \end{cases}$$ Pour fixer $\eta$, écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour ce paramètre. Cela fixera aussi le nombre de points $N$.
In [ ]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

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


nom="Faccanoni"
prenom="Gloria"

L=list(range(-5,0))
idx=sum([ord(c) for c in nom+prenom])%len(L)
eta=L[idx]
display(Math(r'\eta='+sympy.latex(eta)))
N = -15*eta
print(f'{N = }')
$\displaystyle \eta=-1$ 
N = 15

Q1 [1 point]
Calculer la solution exacte.

In [ ]:
t0 = 0
tfinal = 10

y0 = 1
z0 = 1

Calculons la solution exacte avec sympy par exemple:

In [ ]:
t = sympy.Symbol('t')
y = sympy.Function('y')
z = sympy.Function('z')

phi1 = lambda t,y,z : eta*z
phi2 = lambda t,y,z : y+eta*z


edo1 = sympy.Eq( sympy.diff(y(t),t) , phi1(t,y(t),z(t)) )
edo2 = sympy.Eq( sympy.diff(z(t),t) , phi2(t,y(t),z(t)) )
# display(edo1)
# display(edo2)
solgen = sympy.dsolve([edo1,edo2],[y(t),z(t)])
# display(solgen)
consts = sympy.solve( [ sympy.Eq( y0, solgen[0].rhs.subs(t,t0)) , sympy.Eq( z0, solgen[1].rhs.subs(t,t0)) ] , dict=True)[0]
# display(consts)
solpar_1=solgen[0].subs(consts)
solpar_2=solgen[1].subs(consts)
display(solpar_1,solpar_2)
$\displaystyle y{\left(t \right)} = \left(- \frac{\sqrt{3} \sin{\left(\frac{\sqrt{3} t}{2} \right)}}{3} + \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}}$
$\displaystyle z{\left(t \right)} = \left(\frac{\sqrt{3} \left(- \frac{\sin{\left(\frac{\sqrt{3} t}{2} \right)}}{2} + \frac{\sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}}{2}\right)}{3} + \frac{\sqrt{3} \sin{\left(\frac{\sqrt{3} t}{2} \right)}}{2} + \frac{\cos{\left(\frac{\sqrt{3} t}{2} \right)}}{2}\right) e^{- \frac{t}{2}}$

On remarque que $y(t)\xrightarrow[t\to+\infty]{}0$ et $z(t)\xrightarrow[t\to+\infty]{}0$

In [ ]:
func_1 = sympy.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sympy.lambdify(t,solpar_2.rhs,'numpy')

tt = linspace(t0,tfinal,N+1)
h  = tt[1]-tt[0]

figure(figsize=(12,5))
yy=func_1(tt)
zz=func_2(tt)
subplot(1,2,1)
plot(tt,yy,'-*',tt,zz,'-*')
legend([r'$t\mapsto y(t)$',r'$t\mapsto z(t)$'])
subplot(1,2,2)
plot(yy,zz,'-*')
xlabel(r'$y$')
ylabel(r'$z$')
axis('equal');

Dans la suite on notera $u_n\approx y(t_n)$ et $w_n\approx z(t_n)$.

Pour l'affichage des solutions exactes VS approchées, on peut écrire une fonction qui prend en paramètres la discrétisation tt, les deux listes solution exacte yy et zz, les deux listes solution approchée uu et ww et pour finir une chaîne de caractères s qui contient le nom du schéma.

NOTA BENE : on a juste fixé le nombre de points mais nous n'avons pas fixé la taille du domaine. On regardera la solution exacte pour choisir un $T$ final significatif.

In [ ]:
def affichage(tt,yy,zz,uu,ww,s):

    subplot(1,2,1)
    plot(tt,uu,'o--',tt,ww,'d--')
    plot(tt,yy,tt,zz)
    xlabel('t')
    legend([r'$u(t)$ approchée',r'$w(t)$ approchée','$y(t)$ exacte','$z(t)$ exacte'])
    title(f'{s} - y(t) et z(t)') 
    grid()

    subplot(1,2,2)
    plot(uu,ww,'o--')
    plot(yy,zz)
    xlabel('y')
    ylabel('z')
    legend(['Approchée','Exacte'])
    title(f'{s} - z(y)')
    grid()
    axis('equal');

Q2 [3 points]
Calculer la solution approchée obtenue par la méthode d'Euler modifié.
Afficher $t\mapsto y(t)$, $t\mapsto z(t)$ et $y\mapsto z(y)$ en comparant solution exacte et solution approchée.

$$ \text{(EM) } \begin{cases} u_0=1,\\ w_0=0,\\ \tilde u=u_n+\frac{h}{2}\varphi_1(t_n,u_n,w_n),\\ \tilde w=w_n+\frac{h}{2}\varphi_2(t_n,u_n,w_n),\\ u_{n+1}=u_n+h\varphi_1\left(t_n+\frac{h}{2},\tilde u,\tilde w\right),\\ w_{n+1}=w_n+h\varphi_2\left(t_n+\frac{h}{2},\tilde u,\tilde w\right). \end{cases} $$
In [ ]:
# x0, y0, h variables globales
def EM(phi1,phi2,tt):
    uu = [y0]
    ww = [z0]
    for i in range(len(tt)-1):
        uu_tmp = uu[i]+h/2*phi1(tt[i],uu[i],ww[i])
        ww_tmp = ww[i]+h/2*phi2(tt[i],uu[i],ww[i])
        uu.append(uu[i]+h*phi1(tt[i]+h/2,uu_tmp,ww_tmp))
        ww.append(ww[i]+h*phi2(tt[i]+h/2,uu_tmp,ww_tmp))
    return [uu,ww]

[uu, ww] = EM(phi1,phi2,tt)

figure(figsize=(12,5))
affichage(tt,yy,zz,uu,ww,"Euler modifié")

Q3 [3 points]
Calculer la solution approchée obtenue par la méthode de Crank-Nicolson.
Afficher $t\mapsto y(t)$, $t\mapsto z(t)$ et $y\mapsto z(y)$ en comparant solution exacte et solution approchée.

$$ \text{(CN) } \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+\frac{h}{2}\left(\varphi_1(t_{n+1},u_{n+1},w_{n+1})+\varphi_1(t_{n},u_{n},w_{n})\right),\\ w_{n+1}=w_n+\frac{h}{2}\left(\varphi_2(t_{n+1},u_{n+1},w_{n+1})+\varphi_2(t_{n},u_{n},w_{n})\right). \end{cases} $$
In [ ]:
from scipy.optimize import fsolve

# x0,y0, h variables globales
def CN(phi1,phi2,tt):
	uu = [y0]
	ww = [z0]
	for i in range(len(tt)-1):
		sys = lambda z : [ -z[0]+uu[i]+h/2*phi1(tt[i],uu[i],ww[i])+h/2*phi1(tt[i+1],z[0],z[1]) , -z[1]+ww[i]+h/2*phi2(tt[i],uu[i],ww[i])+h/2*phi2(tt[i+1],z[0],z[1]) ]
		utemp,wtemp = fsolve( sys , (uu[i],ww[i]) ) 
		uu.append(utemp)
		ww.append(wtemp)
	return [uu,ww]

[uu, ww] = CN(phi1,phi2,tt)

figure(figsize=(12,5))
affichage(tt,yy,zz,uu,ww,"Crank Nicolson")

Dans ce cas particulier il est possible d'expliciter cette récurrence en exploitant la linéarité des $\varphi_1$ et $\varphi_2$: $$ \text{(CN) } \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+\frac{h}{2}\eta \left(w_{n+1}+w_{n}\right),\\ w_{n+1}=w_n+\frac{h}{2}\left(u_{n+1}+ u_{n}\right)+\frac{h}{2}\eta\left(w_{n+1}+\eta w_{n}\right) \end{cases} $$ ce qui donne le schéma explicite calculé ci-dessous par sympy:

In [ ]:
t_n = sympy.Symbol('t_n')
dt = sympy.Symbol('h')
t_np1 = t_n+dt
u_n = sympy.Symbol(r'u_{n}')
u_np1 = sympy.Symbol(r'u_{n+1}')
w_n = sympy.Symbol(r'w_{n}')
w_np1 = sympy.Symbol(r'w_{n+1}')

Eq1 = sympy.Eq (u_np1, u_n+dt/2*(phi1(t_n,u_n,w_n)+phi1(t_np1,u_np1,w_np1)) )
Eq2 = sympy.Eq (w_np1, w_n+dt/2*(phi2(t_n,u_n,w_n)+phi2(t_np1,u_np1,w_np1)) )
# prlat(Eq1)
# prlat(Eq2)
sol=sympy.solve([Eq1,Eq2],(u_np1,w_np1))

prlat(u_np1,"=",sol[u_np1].factor(u_n))
prlat(w_np1,"=",sol[w_np1].factor(w_n))
$\displaystyle u_{n+1}=- \frac{4 h w_{n} + u_{n} \left(h^{2} - 2 h - 4\right)}{h^{2} + 2 h + 4}$
$\displaystyle w_{n+1}=- \frac{- 4 h u_{n} + w_{n} \left(h^{2} + 2 h - 4\right)}{h^{2} + 2 h + 4}$
In [ ]:
# x0,y0, h variables globales
def CN(tt):
	uu = [y0]
	ww = [z0]
	for i in range(len(tt)-1):
		uu.append( -((h**2-2*h-4)*uu[i]+4*h*ww[i])/(h**2+2*h+4) )
		ww.append( -((h**2+2*h-4)*ww[i]-4*h*uu[i])/(h**2+2*h+4) )
	return [uu,ww]

[uu, ww] = CN(tt)

figure(figsize=(12,5))
affichage(tt,yy,zz,uu,ww,"Crank Nicolson")

Q5 [1 point]
Calculer la solution approchée obtenue par la fonction odeint du module scipy.
Afficher $t\mapsto y(t)$, $t\mapsto z(t)$ et $y\mapsto z(y)$ en comparant solution exacte et solution approchée.

In [ ]:
from scipy.integrate import odeint

pphi = lambda yy,t : [ phi1(t,yy[0],yy[1]), phi2(t,yy[0],yy[1])  ]

yy0  = [y0,z0]
uu,ww = odeint(pphi,yy0,tt).T

figure(figsize=(12,5))
affichage(tt,yy,zz,uu,ww,"Scipy")

Exercice 3 : étude de la A-stabilité d'un schéma

Soit le schéma $$\begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{4}\left(\varphi(t_{n},u_{n})+\varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2})\right),\\ 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}$$

Q1 [3 points] Étudier théoriquement la A-stabilité du schéma suivant

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)=1. \end{cases}$$ Sa solution est $y(t)=e^{-\beta t}$ donc $\lim\limits_{t\to+\infty}|y(t)|=0$.

Le schéma appliqué à ce problème de Cauchy s'écrit $$ \begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n-\frac{\beta h}{4}u_{n}-\frac{\beta h}{4}\tilde u_{n+1/2},\\ u_{n+1}=u_n-\beta h \tilde u_{n+1/2}& n=0,1,2,\dots N-1 \end{cases} $$ d'où $$ u_{n+1} = \left( 1-\beta h\frac{4-\beta h}{4+\beta h} \right)u_n $$

Vérifions ces calculs:

In [ ]:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math

# pour ne pas effacer l'affectation de "h", ici je vais l'appeler "dt" mais afficher "h"
dt = sym.Symbol('h')
beta = sym.Symbol(r'\beta')
u_np1 = sym.Symbol(r'u_{n+1}')
sym.var('u_n,x')

u_tilde = sym.solve (-x + u_n - beta*dt/4*u_n-beta*dt/4*x, x)[0]
display(Math(r'\tilde u='+sym.latex(u_tilde)))

u_np1 = (u_n-beta*dt*u_tilde).factor()
display(Math('u_{n+1}='+sym.latex(u_np1)))
$\displaystyle \tilde u=\frac{u_{n} \left(- \beta h + 4\right)}{\beta h + 4}$
$\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} h^{2} - 3 \beta h + 4\right)}{\beta h + 4}$

On note $x=\beta h$. Le schéma est A-stable ssi $|q(x)|<1$. On étudie alors la fonction $q\colon \mathbb{R}^+\to\mathbb{R}$ définie par $q(x)=\dfrac{x^2-3x+4}{4+x}$.

In [ ]:
from IPython.display import display, Math
import sympy as sym
sym.init_printing()

sym.var('x',nonnegative=True)

print('Dans R^+ on étudie la fonction')
q = (x**2-3*x+4)/(x+4)
display(Math('q(x)='+sym.latex(q)))


print('Limites aux bornes du domaine de definition:')
display(Math('q(0)='+sym.latex(q.subs(x,0))))
lim=sym.Limit(q,x,sym.oo)
display(Math(sym.latex(lim)+'='+sym.latex(lim.doit())))

print('On cherche les points stationnaires dans R^+:')
dq=(q.diff(x)).factor()
display(Math("q'(x)="+sym.latex(dq)))
print("On remarque que le numerateur de q' est une parabole et on a")
sol=sym.solve(dq)
display(Math("q'(x)=0 \iff x\in"+sym.latex(sol)))
print("De plus q'(x)>0 ssi x est supérieur à cette valeur, donc")
minimum=sol[0]
display(Math("x="+sym.latex(minimum)+"\\text{ est un minimum et}"))
qmin=q.subs(x,sol[0])
display(Math("q("+sym.latex(minimum)+")="+sym.latex(qmin)+"="+str(float(qmin))+">0"))
print("Donc q(x)>0 pour tout x.")
print("De plus, on a q(x)<1 ssi")
display(sym.solve(q<1))
print("Conclusion: -1<q(x)<1 pour tout x<4.")
sym.plot(q,1,-1,xlim=[-1,15],ylim=[-1.5,1.5]);

Dans R^+ on étudie la fonction
$\displaystyle q(x)=\frac{x^{2} - 3 x + 4}{x + 4}$ 
Limites aux bornes du domaine de definition:
$\displaystyle q(0)=1$
$\displaystyle \lim_{x \to \infty}\left(\frac{x^{2} - 3 x + 4}{x + 4}\right)=\infty$ 
On cherche les points stationnaires dans R^+:
$\displaystyle q'(x)=\frac{x^{2} + 8 x - 16}{\left(x + 4\right)^{2}}$ 
On remarque que le numerateur de q' est une parabole et on a
$\displaystyle q'(x)=0 \iff x\in\left[ -4 + 4 \sqrt{2}\right]$ 
De plus q'(x)>0 ssi x est supérieur à cette valeur, donc
$\displaystyle x=-4 + 4 \sqrt{2}\text{ est un minimum et}$
$\displaystyle q(-4 + 4 \sqrt{2})=\frac{\sqrt{2} \left(- 12 \sqrt{2} + \left(-4 + 4 \sqrt{2}\right)^{2} + 16\right)}{8}=0.3137084989847604>0$ 
Donc q(x)>0 pour tout x.
De plus, on a q(x)<1 ssi
$\displaystyle 0 < x \wedge x < 4$ 
Conclusion: -1<q(x)<1 pour tout x<4.

Q2 [2 points] Approcher la solution du problème de Cauchy suivant avec le schéma donné d'abord avec 120 points puis avec 130 points. Commenter le résultat. $$\begin{cases} y'(t)=-50y(t), &t\in[0;10]\\ y(0)=1 \end{cases}$$

Pour la A-stabilité il faut que $h=\frac{10-0}{N}$ soit $<\frac{4}{\beta}$ donc $N>\frac{(10-0)\beta}{4}=\frac{10\times50}{4}=125$. De plus, le problème est stiff.

In [ ]:
%reset -f

from matplotlib.pylab import *
from scipy.optimize import fsolve


def schema(tt,phi,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        u_tilde = fsolve ( lambda x : -x+uu[i]+h/4*phi(tt[i],uu[i])+h/4*phi(tt[i+1],x) , uu[i] )[0]
        uu.append( uu[i]+h*phi(tt[i]+h/2,u_tilde) )
    return uu


t0, y0, tfinal = 0, 1, 10
beta=50

phi = lambda t,y : -beta*y
sol_exacte = lambda t : exp(-beta*t)

figure(figsize=(10,5))

subplot(1,2,1)
N = 120
tt = linspace(t0,tfinal,N+1)
h  = tt[1]-tt[0]
    
yy = [sol_exacte(t) for t in tt]    
uu = schema(tt,phi,y0)
plot(tt,yy,'b-',label=("Exacte"))
plot(tt,uu,'ro',label=("RK"))
title(rf' $N$={N}, $h$={h:1.5f}')
xlabel('$t$')
ylabel('$u$')
legend() 
grid(True)

subplot(1,2,2)
N = 130
tt = linspace(t0,tfinal,N+1)
h  = tt[1]-tt[0]
    
yy = [sol_exacte(t) for t in tt]    
uu = schema(tt,phi,y0)
plot(tt,yy,'b-',label=("Exacte"))
plot(tt,uu,'ro',label=("RK"))
title(rf' $N$={N}, $h$={h:1.5f}')
xlabel('$t$')
ylabel('$u$')
legend() 
grid(True);


figure(figsize=(10,5))
H = []
err = []
N = 130
for k in range(7):
    N+=100
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu = schema(tt,phi,y0)
    H.append(h)
    err.append(max([abs(uu[i] - yy[i]) for i in range(len(yy))]))
    
omega=polyfit(log(H),log(err), 1)[0]
figure(figsize=(8,5))
plot(log(H), log(err), 'c-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title(rf'Ordre $\omega$={omega:1.2f}')
grid(True);

<Figure size 720x360 with 0 Axes>