-- Simulation numérique en mécanique des fluides --
TP n°9
Méthode de Newton
éléments de correction
simon.marie@lecnam.net
import numpy as np
import matplotlib.pyplot as plt
#Option d'affichage et taille de police sur les figures:
fs=20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
#plt.rc('text', usetex=True)
%matplotlib inline
1 - Méthode de Newton scalaire
La méthode de Newton permet de trouver numériquement la solution de l'équation $f(x)=0$ en partant d'un point de départ $x_0$ et en itérant la suite:
$$ x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)} $$jusqu'à atteindre un critère de convergence.
Ecrire une fonction permettant de calculer la solution de $f(x)=0$ en connaissant sa dérivée que l'on notera $df$ ainsi que le point de départ que l'on notera $x_0$:
def newton(f,df,x0):
norm=10.
n=0
while (norm>1e-6):
n+=1
norm=abs(f(x0)/df(x0))
x0=x0-norm
print(n,norm)
return x0
Icila norme est définie par l'écart $| x_{n+1}-x_n |$ qui est un scalaire.
def f(x):
return x**2-2
def df(x):
return 2*x
Ici les fonctions sont très simples, on peut donc calculer la dérivée facielement.
df(6.)
On test donc notre première ilmplémentation de la méthode de Newton pour le calcul de $\sqrt{2}$ en partant de $2$:
newton(f,df,2.)
np.sqrt(2.)
Après seulement $5$ itérations on a donc convergé à $10^{-12}$ !!
2 - Méthode de Newton vectorielle¶
La méthode de Newton permet de trouver numériquement la solution de l'équation $F(x,y)=0$ en partant d'un point de départ $X_0=(x_0,y_0)$ et en itérant la suite:
$$ \mathbf{J(X_n)}d\mathbf{X_{n+1}}=-\mathbf{F(X_n)} $$avec:
$$ d\mathbf{X_{n+1}}=\mathbf{X_{n+1}-X_n} $$jusqu'à atteindre un critère de convergence.
1 - Système non linéaire
Considérons le système suivant:
$$ \left\{\begin{array}{lcr} 9x^2+36y^2+4z^2&=&36\\ x^2-2y^2-20z&=&0\\ x^2-y^2+z^2&=&0 \end{array}\right. $$La Jacobienne du système s'écrit: $$ J=\left[\begin{array}{lcr} 18x & 72y & 8z\\ 2x & -4y & -20 \\ 2x & -2y & 2z \end{array}\right] $$
On peut donc définir les fonctions associés facilement. Encore une fois, tout est bien connue ici. La seule nouveauté c'est que l'on résout un système, on a donc une formulation vectoriel faisant intervenir la fonction python solve de numpy.linalg pour résoudre un système linéaire à chaque itérations. Ainsi, résoudre un système non linéaire avec Newton, c'est résoudre plusieurs systèmes linéaires à la suite !!
def F(W):
F1=9*W[0]**2+36*W[1]**2+4*W[2]**2-36
F2=W[0]**2-2*W[1]**2-20*W[2]
F3=W[0]**2-W[1]**2+W[2]**2
return np.array([F1,F2,F3])
def J(W):
J=np.array([[18*W[0],72*W[1],8*W[2]],[2*W[0],-4*W[1],-20],[2*W[0],-2*W[1],2*W[2]]])
return J
def NEWTON(F,dF,W0):
norm=10.
n=0
while (norm>1e-8):
n+=1
DX=np.linalg.solve(dF(W0),-F(W0))
W0=W0+DX
norm=np.sum(DX**2)**0.5
print(n,norm)
return W0
Ici $x_{n+1}-x_n$ est un vecteur. En prendra la norme $L_2$ comme critère de convergence. Essayons en partant d'un triplet initial arbitraire:
W0=[1,1,1]
W=NEWTON(F,J,W0)
print(W)
F(W)
Ici également, en seulement 6 itérations la solution est obtenue à $10^{-11}$ près
2 - EDP non linéaire: L'équation de Burger Stationaire
Considérons maintenant l'EDP non linéaire suivante:
$$ u\dfrac{\partial u}{\partial x}=\nu\dfrac{\partial^2 u}{\partial x^2} $$sous forme conservative, on a:
$$ \dfrac{1}{2}\dfrac{\partial u^2}{\partial x}=\nu\dfrac{\partial^2 u}{\partial x^2} $$Avec les conditions $u(0)=0$ et $u(N)=1$
L'utilisation des schémas centré d'ordre 2 donne:
$$ \dfrac{1}{2}\dfrac{u_{i+1}^2-u_{i-1}^2}{2\Delta x}=\nu\dfrac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2} $$On a ici $F(u)=\dfrac{1}{2}\dfrac{\partial u^2}{\partial x}-\nu\dfrac{\partial^2 u}{\partial x^2}$
avec:
$\dfrac{d F}{d u}=u\dfrac{\partial .}{\partial x}-\nu\dfrac{\partial^2 .}{\partial x^2}$
La nouveauté ici, c'est que la dérivé de la fonction est plus compliqué et non connue à priori. Cette Jacobienne dépend du schéma utilisé et s'écrit sous la forme d'un "opérateur" comme nous avions écrit l'opérateur Laplacien dans le TP8.
def F(W,dx,nu):
nx=W.shape[0]
ff=np.zeros(int(nx))
ff[1:-1]=0.25*(W[2:]**2-W[:-2]**2)/dx-nu*(W[2:]-2.*W[1:-1]+W[:-2])/dx**2
return ff
def J(W,dx,nu):
nx=W.shape[0]
J2=-2*np.eye(nx)+np.diag(np.ones(nx-1),-1)+np.diag(np.ones(nx-1),1)
J1=0.5*np.diag(W[1:],1)-0.5*np.diag(W[:-1],-1)
J0=J1/dx-nu*J2/dx**2
return J0
def NEWTON(F,dF,W0,dx=0.1,nu=0.1,nmax=10000):
norm=10.
n=0
while (abs(norm)>1e-12 and n<nmax):
n+=1
DX=np.linalg.solve(dF(W0,dx,nu),-F(W0,dx,nu))
W0=W0+DX
norm=np.sum(DX**2)**0.5
print(n,norm)
return W0
On choisit donc un domaine de simulation sur lequel on recherche la solution de l'EDP précédente. On se donne un nombre de points, et un $\Delta x$. On fixe ensuite le coefficient $\nu$.
La quantité $\dfrac{U_0\Delta x}{\nu}$ est alors importante. Elle doit rester la plus faible possible. Cette quantité correspond au nombre de Reynolds basé sur l'echelle de la taille de maille. on parle de Reynolds de maille, on le note $R_{em}$. Ici on initialisera avec une fonction linéaire:
L0=1.
dx=0.001
nx=int(L0/dx)
x=np.linspace(0,L0,nx)
u0,un=0,1
W0=np.linspace(u0,un,nx)
On choisit 3 valeurs différentes pour la viscosité:
nu1,nu2,nu3=1.,0.1,0.01
W1=NEWTON(F,J,W0,dx,nu1)
W2=NEWTON(F,J,W0,dx,nu2)
W3=NEWTON(F,J,W0,dx,nu3)
Encore une fois on observe une convergence rapide obtenu en quelques itérations en remarquant que si le Reynolds de maille devient grand, la convergence est plus longue.
fig=plt.figure(figsize=(16,6))
plt.plot(x,W1,label='$\\nu_1=$'+str(nu1))
plt.plot(x,W2,label='$\\nu_2=$'+str(nu2))
plt.plot(x,W3,label='$\\nu_3=$'+str(nu3))
plt.ylim(0,1);
plt.grid(True);plt.legend(fontsize=fs)
On visualise le "Résidu" de la solution c'est à dire la fonction $F$ évaluée au point solution $W$. Normalement cette quantité doit être nulle. En pratique les endroits ou elle n'est pas strictement nulle indique que la solution est moins "fiable".
fig=plt.figure(figsize=(16,6))
plt.plot(x,F(W1,dx,nu1),label='Rem='+str(dx/nu1))
plt.plot(x,F(W2,dx,nu2),label='Rem='+str(dx/nu2))
plt.plot(x,F(W3,dx,nu3),label='Rem='+str(dx/nu3))
plt.grid(True);plt.legend(fontsize=fs)
On pourra essayer de changer le Reynolds de Maille et les conditions aux limites et observer que la convergence est mise en défaut si le Reynolds de maille devient grand.
3- Solutions analytiques¶
L'EDP (1) admet une solution analytique du type:
$$ u(x)=k\sqrt{2\nu}\tan\left(k\dfrac{\sqrt{2\nu} x}{2\nu} \right) $$Avec k compatible avec les conditions aux limites $U(x=L)=1$, soit:
$$ k\sqrt{2\nu}\tan\left(k\dfrac{\sqrt{2\nu} L}{2\nu} \right)=1 $$Pour trouver la constante compatible avec les conditions aux limites du problème, on peut alors utiliser une méthode de Newton !!!
def const(x,L,nu):
a=np.sqrt(2*nu)
return x*a*np.tan(x*L/a)-1
def dconst(x,L,nu):
a=np.sqrt(2*nu)
return a*np.tan(x*L/a)+x*L/(np.cos(x*L/a)**2)
def newton(f,df,x0,L,nu):
norm=10.
n=0
while (norm>1e-9 and n<1000):
n+=1
norm=abs(f(x0,L,nu)/df(x0,L,nu))
x0=x0-norm
#print(n,norm,k)
return x0
Pour calculer la solution analytique on prendra le cas $\nu=0.1$. Le principal, problème ici, c'est que l'équation de départ admet de multiples solutions. Seules certaines sont compatibles avec la solution numérique que donne la méthode de Newton.
nuth=nu3
kn=np.arange(0.001,2.,0.001)
k0=np.zeros(len(kn))
for i,k in enumerate(kn):
k0[i]=newton(const,dconst,k,L=L0,nu=nuth)
kpos=k0[k0>0]
kpos_init=kn[k0>0]
fig=plt.figure(figsize=(16,6))
plt.plot(kn[k0>0],k0[k0>0],'.')
plt.xlim(0,1);plt.ylim(-1,1);plt.grid(True);
k=newton(const,dconst,kpos_init[0],L=L0,nu=nuth)
a=np.sqrt(2*nuth)
Wth=k*a*np.tan(k*x/a)
print(k*a*np.tan(k*L0/a))
fig=plt.figure(figsize=(16,6))
plt.plot(x,W3,label='Newton $\\nu=$'+str(nu3))
plt.plot(x[::7],Wth[::7],'o',label='theorique $\\nu=$'+str(nuth)+' (k='+str(k)+')')
plt.ylim(0,1);
plt.grid(True);plt.legend(fontsize=fs)
Conclusion¶
Présentez une synthèse de vos résultats et des étapes de calcul de votre TP.
from IPython.core.display import HTML
style=open('../notebooks.css', "r").read()
HTML(style)