diffusionPerturb.py
#
# Analyse de stabilité de l'opérateur discret du laplacien
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
from scipy.sparse import spdiags
from scipy.linalg import eig
import matplotlib
from matplotlib import pyplot as plt
matplotlib.rcParams['toolbar'] = 'None'
myColorMap = matplotlib.cm.jet
#
# -0- Opérateur discret du laplacien
# Nombre de noeuds = m
# Taille des mailles dx = L/(m-1)
# Deux conditions essentielles
# Taille de l'opérateur discret m-2
# Dimension de l'espace discret m-2
#
# Exemple classique : m = 11 (m = 18 is also nice)
#
m = 11; dx = 1.0/(m-1); beta = (dx*dx)/2
e = array([0,*ones(m-2),0])
D = spdiags([e,-2*e,e],[-1,0,1],m,m)
D = D.T/(dx*dx)
#
# -1- Valeurs propres de l'opérateur discret du laplacien
# dans le diagramme de stabilité d'Euler explicite
#
lam,eigen = eig(beta * D.toarray()[1:-1,1:-1])
plt.figure("Stability of Euler with central finite differences")
x,y = meshgrid(linspace(-3,1,1000),linspace(-1.5,1.5,1000))
z = x + 1j*y
gain = abs(1.0 + z)
plt.contourf(x,y,gain,[0,1,100],colors=('lightgreen','none'))
plt.contour(x,y,gain,[1],colors='black',linewidths=0.5)
plt.gca().axhline(y=0,color='gray',linewidth=0.5)
plt.gca().axvline(x=0,color='gray',linewidth=0.5)
plt.plot(lam.real,lam.imag,'or')
plt.xticks(arange(-2.5,1,0.5))
plt.yticks(arange(-1,1.5,0.5))
plt.axis('equal')
#
# -2- Vecteurs propres de l'opérateur discret laplacien
# On les trie et on conserver l'ordre avec argsort()
#
msqrt = int(ceil(sqrt(m-2)))
order = argsort(-lam.real)
plt.figure("Eigenvectors")
X = linspace(0,1,m); V = zeros(m)
for i in range(m-2):
iloc = order[i]
V[1:-1] = eigen[:,iloc]
plt.subplot(msqrt,msqrt,i+1)
plt.plot(X,V,'.r',markersize=7,linewidth=1)
plt.title("%6.4f" % lam[iloc].real)
plt.xlim((-0.1,1.1)); plt.ylim((-2.1,2.1)); plt.axis('off')
#
# -4- Quelques peturbations bien choisies :-)
#
x = linspace(0,1,100);
X = linspace(0,1,m);
A = sin(outer(X[1:5],arange(1,5))*pi)
plt.figure("Some nice perturbations")
for i in range(m-2):
u = sin((i+1)*pi*x); U = sin((i+1)*pi*X)
plt.subplot(msqrt,msqrt,i+1)
plt.plot(x,u,'-b',[0,1],[0,0],'-k',X,U,'.r',markersize=7,linewidth=1)
plt.xlim((-0.1,1.1)); plt.ylim((-1.1,1.1)); # plt.axis('off')
plt.xticks([]); plt.yticks([])
#
# -5- Analyse de stabilité des perturbations
#
plt.figure("Stability analysis")
for i in range(m-2):
u = sin((i+1)*pi*x); U = sin((i+1)*pi*X)
plt.subplot(msqrt,msqrt,i+1)
plt.plot(x,u,'-b',X,U,'.r',markersize=7,linewidth=1)
Unew = -beta * D @ U
gain = Unew[1]/U[1]
plt.plot(X,Unew,'.g',markersize=7,linewidth=1)
plt.title("%6.4f" % gain)
plt.xlim((-0.1,1.1)); plt.ylim((-2.1,2.1)); plt.axis('off')
plt.show()