laplacien discret (python)

[Brinicle]

Bonjour,

Je dois programmer une fonction renvoyant la matrice du laplacien discret en 2D (dans la but d'avoir ses valeurs et vecteurs propres). évidemment, mon maillage est constitué de triangles plus ou moins quelconque, je vous ai mis en image l'énoncé où se trouve la formule à implémenter.

J'ai testé sur un rectangle pour comparer avec la formule analytique. Bien évidement mon programme ne marche pas... les valeurs propres sont délirantes, mais bizarrement, les vecteurs propres semblent corrects

Je ne vois vraiment pas ce qui ne marche pas, j'ai mis mon code en pièce jointe.

Le laplacien discret est la fonction "M" dans le fichier "fonction_in_curve"
-----

[JPL]

Le codes doivent être intégrés dans le message grâce aux balises [Code]...[/Code] et non dans un zip. Merci

[Brinicle]

Le code est un peu long, c'est pour ça que je l'avais mis en zip... je le remet en texte :

Contient toutes les fonctions (il faut nommer la fichier "fonction_in_curve) :

import numpy as np
import scipy.spatial

def rect_tambour(L,H,U):
    """
    Paramètres
    ----------
    L : Largeur du tambour
    
    H : Hauteur du tambour
    
    U : Valeur propre maximale voulue

    Retourne
    -------
    Modes propres de valeur propre supérieure à -U.
    """
    vals = []
    val = 0
    k = 1
    l = 1
    while val >= -U:
        while val >= -U:
            #print(val,l,k)
            vals.append(val)
            val = -np.pi**2*((k/L)**2+(l/H)**2)
            l += 1
        l = 2
        k += 1
        val = -np.pi**2*((k/L)**2+(l/H)**2)
    return -np.sort(-np.array(vals)[1:])

def count_vp(tab,U):
    #compte les n valeurs propres upérieures au égales à -U dans le tableau tab
    return tab[tab>=-U]

def in_curve(f,fargs,shape,a):
    """
    Permet de trouver les point situés à l'intérieur de la coube délimitée par l'équation f(x,y,*fargs)=0
    
    Paramètres.
    ----------
    f : fonction de type f(x,y) définissant l'intérieur de la courbe par f(x,y,fargs)<=0
    Exemple, cercle de rayon R : f(x,y,R) = x**2+y**2-R
    
    fargs : Autre arguments de f. Exemple, cercle de rayon R : fargs=[R]
    
    shape : dimensions du tableau représentant l'espace où la coube est calculée.
    
    a : distance représentée par chaque case.
    
    eps : précision pour la détermination du bord

    Retourne
    -------
    Tableau de points.
    
    """
    points = [] # les points à l'intérieur de la courbe
    for j in range(shape[0]):
        for i in range(shape[1]):
            if f(i*a,j*a,*fargs)<0:
                points.append([i*a,j*a])

    return np.array(points)

def triang(points,a,f,fargs,bord):
    """
    Transforme un réseau carré en un réseau triangulaire.
    
    Paramètres
    ----------
    points : tableau de points.
    a : côté des carré.
    f : intérieur de la courbe.
    fargs : Autre arguments de f.
    bord : bord de la courbe.

    Retourne
    -------
    tableau de points.

    """
    tri_points = points.copy()
    tri_points[:,1] *= np.sqrt(3)
    tri_points2 = np.vstack((points,bord))
    tri_points2[:,1] *= np.sqrt(3)
    tri_points2[:,0] += a/2
    tri_points2[:,1] += np.sqrt(3)/2*a
    fin = np.vstack((tri_points,tri_points2))
    i = 0
    eps = 0.01
    while i < len(fin):
        if f(fin[i,0]+eps,fin[i,1]+eps,*fargs) > 0:
            fin = np.delete(fin,i,0)
            i -= 1
        i += 1
    return np.vstack((fin,bord)),len(fin),len(bord)

def tri_ang(points,ind,p0):
    # trie les points dans l'ordre trigo
    vec=np.arctan2((points-p0)[:,1],(points-p0)[:,0]) # les angles de chaque point
    values = []
    dtype = [('val',float),('n',int)]
    # on associe à chaque valeur un entier afin de retrouver à quel point ils appartiennent après le tri
    for i in range(len(vec)):
        values.append((vec[i],i))
    values = np.sort(np.array(values,dtype),order='val') # on trie les angles
    new_points = []
    new_ind = []
    # on utilise les entier afin de trier la liste de points et les indices :
    for tup in values:
        new_points.append(points[tup[1]])
        new_ind.append(ind[tup[1]])
    return np.array(new_points),new_ind





def M(points,tri,Nint,dl):
    
    indptr,ind = tri.vertex_neighbor_vertices
    N = len(points)
    M = np.zeros((N,N))
    
    for i in range(Nint):
        tot = 0
        nhb_ind = ind[indptr[i]:indptr[i+1]] # indices des points voisin du point d'indice k
        nhb = points[nhb_ind] # leurs coordonnées
        nhb,nhb_ind = tri_ang(nhb,nhb_ind,points[i])
        # tiers de la somme des aires : 
        A = 0 # pour les aires
        for j in range(len(nhb_ind)):
            vec = points[nhb_ind[j]]-points[i] # un veteur reliant le point à un voisin
            vec_av = points[nhb_ind[j-1]]-points[i] # un autre vecteur mais avec le vosin d'avant
            # on utilise le produit vectoriel pour calculer les aires des triangles : ||A^B||/2
            A += np.linalg.norm(np.cross(vec,vec_av))/2
        #print(A)
        A /= 3
        
        for j in range(len(nhb_ind)):
            if nhb_ind[j] < Nint:
                vec = points[nhb_ind[j]]-points[i] # un veteur reliant le point à son voisin d'indice 0
                vec_av = points[nhb_ind[j-1]]-points[i] # un autre vecteur mais avec le vosin d'avant
                if j+1 >= len(nhb_ind):
                    k = 0
                else:
                    k = j+1
                vec_ap = points[nhb_ind[k]]-points[i] # un autre vecteur mais avec le voisin d'après
                
                # on utilise le produit vectoriel pour calculer les sinus : ||A^B||/(||A|| ||B||)
                sin_alpha = np.linalg.norm(np.cross(vec_av,vec_av-vec))/(np.linalg.norm(vec_av)*np.linalg.norm(vec_av-vec))
                sin_beta = np.linalg.norm(np.cross(vec_ap-vec,vec_ap))/(np.linalg.norm(vec_ap)*np.linalg.norm(vec_ap-vec))
                
                # on utilise le produit scalaire pour calculer les cosinus : ||A.B||/(||A|| ||B||)
                cos_alpha = np.dot(vec_av,vec_av-vec)/(np.linalg.norm(vec_av)*np.linalg.norm(vec_av-vec))
                cos_beta = np.dot(vec_ap-vec,vec_ap)/(np.linalg.norm(vec_ap)*np.linalg.norm(vec_ap-vec))
                
                cotan_alpha = np.abs(cos_alpha/sin_alpha)
                cotan_beta = np.abs(cos_beta/sin_beta)
                #print(cotan_alpha,cotan_beta)
                #print(nhb,vec_av,vec_ap,vec,points[i])
                M[i,nhb_ind[j]] = -1/(2*A)*(cotan_alpha+cotan_beta)
        
                tot += cotan_alpha+cotan_beta
    
        M[i,i] = -tot/(2*A) # les valeurs pour la diagonale
    
    return M

def rect(x,y,L,H,x0=0,y0=0):
    if 0<x-x0<L and 0<y-y0<H:
        return -1
    else:
        return 1

def rect_bord(L,H,a,x0=0,y0=0):
    tab1 = np.arange(x0,L+x0,a)[:,np.newaxis]
    h = np.hstack((tab1,H*np.ones((len(tab1),1))+y0))
    b = np.hstack((tab1,np.zeros((len(tab1),1))+y0))
    tab2 = np.arange(y0+a,H+y0,a)[:,np.newaxis]
    g = np.hstack((np.zeros((len(tab2),1))+x0,tab2))
    d = np.hstack((L*np.ones((len(tab2),1))+x0,tab2))
    hp = np.array([[L+x0,H+y0]])
    bp = np.array([[L+x0,0]])
    return np.vstack((h,b,g,d,hp,bp))

Le code pour les valeurs et vecteurs propres (il faut l’exécuter dans un autre fichier) :

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import fonction_in_curve as tamb
import scipy.spatial


L = 1
H = 1

# cette partie prépare une maille triangulaire
dl = 0.05
sol = tamb.in_curve(tamb.rect,[L,H],(100,100),dl)
sol_tri,Nint,Nbord = tamb.triang(sol,dl,tamb.rect,[L,H],tamb.rect_bord(L,H,dl))

# plt.plot(sol_tri[:,0],sol_tri[:,1],linestyle="",marker="+",label="tri")
# plt.plot(sol[:,0],sol[:,1],linestyle="",marker="x")
# plt.legend()
# plt.show()

# triangulation
tri = scipy.spatial.Delaunay(sol_tri)

# enlever les "#" si dessous pour voir à quoi ressemble la maille
# plt.triplot(sol_tri[:,0],sol_tri[:,1],tri.simplices)
# plt.show()

M = tamb.M(sol_tri,tri,Nint,dl)

valp,vecp = np.linalg.eig(M)
valp = np.real(valp)
vecp = np.real(vecp)

# comparaison avec la solution exacte :
eps = 1e-10
T = 1000
U = np.arange(0,T,1)
valp = -np.sort(-valp[valp < -eps]) # on retire les valeurs propres (quasi-) nulles
NUsim = np.array([len(tamb.count_vp(valp,u)) for u in U]) # la simulation
NU = np.array([len(tamb.rect_tambour(L,H,u)) for u in U]) # exactes
plt.plot(U,NUsim,label='simulation')
plt.plot(U,NU,label='exactes')
plt.legend()
plt.show()


# plot 3D d'un vecteur propre
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(sol_tri[:,0],sol_tri[:,1],vecp[:,0],triangles=tri.simplices)
plt.show()

Désolé pour les morceaux de code en commentaire...

[Brinicle]

Bon, j'ai fini par trouver où se trouvaient les problèmes .

Il s'agissait de plusieurs problèmes :

- Je multipliais les colonnes par les colonnes par Ai alors qu'il fallait le faire suivant les lignes (pour que le produit matriciel avec un vecteur place ensuite 1/2Ai à la coordonnées i du vecteur U).

- La formule de l'énoncé est bonne (Uj-Ui), il y avait donc une erreur de signe dans mon code.

- Un problème de bord : il fallait compter l'aire du triangle même si un de ses points était sur le bord.


J'ai aussi optimisé le code (moins de boucles).

Code complet :

Code pour "fonction_in_curve" :

import numpy as np
import scipy.spatial
import scipy.sparse as sparse

def rect_tambour(L,H,U):
    """
    Paramètres
    ----------
    L : Largeur du tambour
    
    H : Hauteur du tambour
    
    U : Valeur propre maximale voulue

    Retourne
    -------
    Modes propres de valeur propre supérieure à -U.
    """
    vals = []
    val = 0
    k = 1
    l = 1
    while val >= -U:
        while val >= -U:
            val = -np.pi**2*((k/L)**2+(l/H)**2)
            if val >= -U:
                vals.append(val)
            l += 1
        l = 1
        k += 1
        val = -np.pi**2*((k/L)**2+(l/H)**2)
    return np.array(vals)

def count_vp(tab,U):
    #compte les n valeurs propres upérieures au égales à -U dans le tableau tab
    return tab[tab>=-U]

def in_curve(f,fargs,shape,a):
    """
    Permet de trouver les point situés à l'intérieur de la coube délimitée par l'équation f(x,y,*fargs)=0
    
    Paramètres.
    ----------
    f : fonction de type f(x,y) définissant l'intérieur de la courbe par f(x,y,fargs)<=0
    Exemple, cercle de rayon R : f(x,y,R) = x**2+y**2-R
    
    fargs : Autre arguments de f. Exemple, cercle de rayon R : fargs=[R]
    
    shape : dimensions du tableau représentant l'espace où la coube est calculée.
    
    a : distance représentée par chaque case.

    Retourne
    -------
    Tableau de points.
    
    """
    points = [] # les points à l'intérieur de la courbe
    for j in range(shape[0]):
        for i in range(shape[1]):
            if f(i*a,j*a,*fargs) < 0:
                points.append([i*a,j*a])

    return np.array(points)

def in_curve2(points,f,fargs):
    in_points = []
    for p in points:
        if f(p[0],p[1],*fargs) < 0:
            in_points.append(p)
    return np.array(in_points)
              

def triang(points,a,f,fargs,bord):
    """
    Transforme un réseau carré en un réseau triangulaire.
    
    Paramètres
    ----------
    points : tableau de points.
    a : côté des carré.
    f : intérieur de la courbe.
    fargs : Autres arguments de f.
    bord : bord de la courbe.

    Retourne
    -------
    tableau de points.

    """
    tri_points = points.copy()
    tri_points[:,1] *= np.sqrt(3)
    tri_points2 = np.vstack((points,bord))
    tri_points2[:,1] *= np.sqrt(3)
    tri_points2[:,0] += a/2
    tri_points2[:,1] += np.sqrt(3)/2*a
    fin = np.vstack((tri_points,tri_points2))
    i = 0
    eps = 0.01
    while i < len(fin):
        if f(fin[i,0]+eps,fin[i,1]+eps,*fargs) > 0:
            fin = np.delete(fin,i,0)
            i -= 1
        i += 1
    return np.vstack((fin,bord)),len(fin),len(bord)


def tri_ang(points,ind,p0):
    # trie les points dans l'ordre trigo
    vec=np.arctan2((points-p0)[:,1],(points-p0)[:,0]) # les angles de chaque point
    values = []
    dtype = [('val',float),('n',int)]
    # on associe à chaque valeur un entier afin de retrouver à quel point ils appartiennent après le tri
    for i in range(len(vec)):
        values.append((vec[i],i))
    values = np.sort(np.array(values,dtype),order='val') # on trie les angles
    new_points = []
    new_ind = []
    # on utilise les entier afin de trier la liste de points et les indices :
    for tup in values:
        new_points.append(points[tup[1]])
        new_ind.append(ind[tup[1]])
    
    return np.array(new_points),np.array(new_ind)



def Laplacien(points,tri,Nint):
    
    indptr,ind = tri.vertex_neighbor_vertices
    W = np.zeros((Nint,Nint))
    A = np.zeros((Nint,1))
    
    for i in range(Nint):
        tot = 0
        nhb_ind = ind[indptr[i]:indptr[i+1]] # indices des points voisin du point d'indice k
        nhb = points[nhb_ind] # leurs coordonnées
        nhb,nhb_ind = tri_ang(nhb,nhb_ind,points[i])
        
        for j in range(len(nhb_ind)):
            vec = nhb[j]-points[i] # un veteur reliant le point à son voisin d'indice 0
            vec_av = nhb[j-1]-points[i] # un autre vecteur mais avec le vosin d'avant
            if j+1 >= len(nhb_ind):
                k = 0
            else:
                k = j+1
            vec_ap = nhb[k]-points[i] # un autre vecteur mais avec le voisin d'après
            
            # on utilise le produit vectoriel pour calculer les aires des triangles : ||A^B||/2 :
            A[i] += 0.5/3*np.linalg.norm(np.cross(vec,vec_av))
            
            # on utilise le produit vectoriel et scalaire pour calculer les cotangentes : A.B/||A^B||
            cotan_alpha = np.dot(vec_av,vec_av-vec)/np.linalg.norm(np.cross(vec_av,vec_av-vec))
            cotan_beta = np.dot(vec_ap,vec_ap-vec)/np.linalg.norm(np.cross(vec_ap,vec_ap-vec))
            
            tot += cotan_alpha+cotan_beta
            
            if nhb_ind[j] < Nint:
                W[i,nhb_ind[j]] = 0.5*(cotan_alpha+cotan_beta)
        
        W[i,i] = -0.5*tot # les valeurs pour la diagonale
    
    return (1/A)*W


def rect(x,y,L,H,x0=0,y0=0):
    if 0<x-x0<L and 0<y-y0<H:
        return -1
    else:
        return 1

def rect_bord(L,H,a,x0=0,y0=0):
    tab1 = np.arange(x0,L+x0,a)[:,np.newaxis]
    h = np.hstack((tab1,H*np.ones((len(tab1),1))+y0))
    b = np.hstack((tab1,np.zeros((len(tab1),1))+y0))
    tab2 = np.arange(y0+a,H+y0,a)[:,np.newaxis]
    g = np.hstack((np.zeros((len(tab2),1))+x0,tab2))
    d = np.hstack((L*np.ones((len(tab2),1))+x0,tab2))
    hp = np.array([[L+x0,H+y0]])
    bp = np.array([[L+x0,0]])
    return np.vstack((h,b,g,d,hp,bp))

Code pour le rectangle :

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import fonction_in_curve as tamb
import scipy.spatial
import scipy.sparse.linalg as splin


L = 1
H = 1

dl = 0.04
sol = tamb.in_curve(tamb.rect,[L,H],(100,100),dl)
sol_tri,Nint,Nbord = tamb.triang(sol,dl,tamb.rect,[L,H],tamb.rect_bord(L,H,dl))

# plt.plot(sol_tri[:,0],sol_tri[:,1],linestyle="",marker="+",label="tri")
# plt.plot(sol[:,0],sol[:,1],linestyle="",marker="x")
# plt.legend()
# plt.show()

# triangulation
tri = scipy.spatial.Delaunay(sol_tri)
# plt.triplot(sol_tri[:,0],sol_tri[:,1],tri.simplices)
# plt.show()

M = tamb.Laplacien(sol_tri,tri,Nint)

valp,vecp = np.linalg.eig(M)
valp = np.real(valp)
vecp = np.real(vecp)

# comparaison avec la solution exacte :
eps = 1e-10
T = 1000
U = np.arange(0,T,1)
NUsim = np.array([len(tamb.count_vp(valp,u)) for u in U]) # la simulation
NU = np.array([len(tamb.rect_tambour(L,H,u)) for u in U]) # exactes
plt.plot(U,NUsim,label='simulation')
plt.plot(U,NU,label='exactes')
plt.legend()
plt.show()


# plot 3D d'un vecteur propre
vecp_tot = np.vstack((vecp,np.zeros((Nbord,Nint))))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(sol_tri[:,0],sol_tri[:,1],vecp_tot[:,0],triangles=tri.simplices)
plt.show()

Il y a quand même un problème avec des dl trop petits, certains modes sont aberrants (et leurs valeurs propres avec), mais les modes propres principaux sont bons, c'est sans doute la qualité du maillage qui provoque cela (il y a des effets de bord).

[A voir en vidéo sur Futura]