In [14]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('figure',figsize=(14,9))
In [15]:
from IPython.html.widgets import interactive
from IPython.display import display
Résolution de l'équation de la chaleur pour mettre en valeur les erreurs :¶
$\displaystyle \frac{\partial u}{\partial t} - \varepsilon \frac{\partial^2 u}{\partial x^2} = 0\quad\hbox{dans } ]-1,1[ \times ]0,1[\; .$¶
$\displaystyle u(x,0) = 0 \quad\hbox{dans }]-1,1[ $¶
$ \displaystyle u(-1,t) = u(1,t) = 1\quad\hbox{pour tous }t \in ]0,1[$¶
L'erreur est le "1" sur le bord du domaine (pour $x=+/-1$) : on fait une erreur d'ordre 1. Le $\varepsilon$ (petit) permet de remplacer l'intervalle $[-R,R]$ avec $R$ grand par l'intervalle $[-1,1]$. Normalement la solution $u$ (qui représente l'erreur) serait nulle si on n'avait pas "1" sur le bord mais "0". On teste donc l'influence de cette valeur sur le bord.¶
In [21]:
T=1.0
eps=0.01
N=1000
J=100
x=np.linspace(-1.0,1.0,J+1)
A=2*np.eye(J-1)
for i in range(J-2):
A[i,i+1]=-1
A[i+1,i]=-1
dx=2.0/J
dt=T/N
ll=dt/(dx*dx)
print(dx,dt,ll)
b=np.matrix(np.zeros((J-1,1)))
b[0,0]=ll
b[J-2,0]=ll
Schéma explicite :¶
In [22]:
def schExp(eps,ll,J,N,A,b):
v=np.matrix(np.zeros((J-1,1)))
u=np.matrix(np.zeros((J-1,1)))
for n in range(1,N):
v=u-eps*ll*(A*u)+eps*b
u=v
return u
In [23]:
u=schExp(eps,ll,J,N,A,b)
x=np.linspace(-1.0,1.0,J+1)
w=np.zeros(J+1)
w[0]=1.0
w[J]=1.0
for i in range(1,J):
w[i]=-np.log(u[i-1])
plt.plot(x,w)
Out[23]:
Cette figure représente l'erreur en échelle logarithmique : elle est donc de l'ordre de $e^{-25}\sim 10^{-14}$ en $x=0$ pour $t=1$!¶
On voit que sur l'intervalle [-0.5,0.5], l'erreur est inférieure à $e^{-7.5}\sim 10^{-3}$, toujours pour $t=1$. Donc acceptable vu la discrétisation peut fine que nous utilisons.¶
Schéma implicite :¶
In [24]:
def schImp(eps,ll,J,N,A,b):
v=np.matrix(np.zeros((J-1,1)))
u=np.matrix(np.zeros((J-1,1)))
M=np.eye(J-1)+eps*ll*A
for n in range(1,N):
v=np.linalg.solve(M,u+eps*b)
u=v
return u
In [25]:
u=schImp(eps,ll,J,N,A,b)
x=np.linspace(-1.0,1.0,J+1)
w=np.zeros(J+1)
w[0]=1.0
w[J]=1.0
for i in range(1,J):
w[i]=-np.log(abs(u[i-1]))
plt.plot(x,w)
Out[25]:
