<jupyter_start><jupyter_text>Geometría Diferencial DiscretaProyecto final<NAME> & <NAME>Se encontrará una solución numérica al siguiente sistema.Sea $$\Omega=\{(x,y):x^{2}+y^{2}\leq1\}$$Encontrar $f:\Omega\rightarrow\mathbb{R}$ tal que \begin{align*}\Delta f & =-1\\f(x,y) & =10 \quad \text{ para } (x,y)\in\partial\Omega\end{align*}Además comparar gráficamente con la solución exacta$$f(x,y)=\frac{1-x^{2}-y^{2}}{4}+10$$Para ello se importará una malla que sirva como discretización de $\Omega$. Sobre esta se calculará la discretización del laplaciano y utilizando mínimos cuadrados se resolverá finalmente el sistema.---  Librerías e importación de datos<jupyter_code># Librerias
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
import numpy as np
import networkx as nx # Libreria de graficas
import math<jupyter_output><empty_output><jupyter_text>Cargar la primera malla (la más burda)<jupyter_code>a = np.loadtxt("aristas_malla1.txt",delimiter=",",dtype=int) - 1
c = np.loadtxt("caras_malla1.txt",delimiter=",",dtype=int) - 1
# Leyendo geometria de la malla
v = np.loadtxt("coordenadas_malla1.txt",delimiter=",")<jupyter_output><empty_output><jupyter_text>Cargar la segunda malla<jupyter_code>a = np.loadtxt("aristas_malla2.txt",delimiter=",",dtype=int) - 1
c = np.loadtxt("caras_malla2.txt",delimiter=",",dtype=int) - 1
# Leyendo geometria de la malla
v = np.loadtxt("coordenadas_malla2.txt",delimiter=",")<jupyter_output><empty_output><jupyter_text>Cargar la tercera malla (la más fina)<jupyter_code>a = np.loadtxt("aristas_malla3.txt",delimiter=",",dtype=int) - 1
c = np.loadtxt("caras_malla3.txt",delimiter=",",dtype=int) - 1
# Leyendo geometria de la malla
v = np.loadtxt("coordenadas_malla3.txt",delimiter=",")<jupyter_output><empty_output><jupyter_text>Discretización del Laplaciano Matriz de primera derivada exterior ($d_0$)<jupyter_code># Se crea una matriz de ceros
d_0= np.zeros([len(a), len(v)]) 

# Se asignan las orientaciones de las aristas
for i in range(len(a)):
  d_0[i][ a[i][0] ] = -1
  d_0[i][ a[i][1] ] = 1
d_0<jupyter_output><empty_output><jupyter_text>Matriz de la estrella de Hodge para 1 formas ($*_1$)<jupyter_code># Funcion que calcula el circuentro de un triangulo
def circuncentro ( t ): # t es una lista con los numeros de los vertices que forman el triangulo
  # Coordenadas de los puntos
  A=v[t[0]]
  B=v[t[1]]
  C=v[t[2]]
  # Calculo del circuncentro
  D=2*( A[0]*(B[1]-C[1])+B[0]*(C[1]-A[1])+C[0]*(A[1]-B[1]) )
  O_x=( (A[0]*A[0]+A[1]*A[1])*(B[1]-C[1]) + (B[0]*B[0]+B[1]*B[1])*(C[1]-A[1]) + (C[0]*C[0]+C[1]*C[1])*(A[1]-B[1]) )/D
  O_y=( (A[0]*A[0]+A[1]*A[1])*(C[0]-B[0]) + (B[0]*B[0]+B[1]*B[1])*(A[0]-C[0]) + (C[0]*C[0]+C[1]*C[1])*(B[0]-A[0]) )/D
  # Circuncentro
  return np.array([O_x,O_y])

# Matriz de circuncentros para cada triángulo en la matriz c
circuncentros=[]
for i in c:
  circuncentros.append( circuncentro(i) )
# Funcion distancia en R2
def distancia( a , b):
  return math.sqrt( (a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]) )

# Funcion que devuelve los indices de los puntos (en circuncentros[]) que forman la arista dual dada la original
def arista_dual( a_orig ):
  a_set = set(a_orig)
  r = []
  global c
  for i in range(len(c)):
    cara = set(c[i])
    if ( len( a_set & cara )==2 ):
      r.append(i)
  return r

# Funcion que devuelve la longitud de arista dual dada la arista original a_orig
def long_dual( a_orig ):
  # Se obtienen los puntos que definen la arista dual
  p=arista_dual(a_orig)
  if ( len(p) == 0 ):
    print("Error")
  elif ( len(p) == 1 ):
    a=v[a_orig[0]]
    b=v[a_orig[1]]
    medio=[ (a[0]+b[0])/2. , (a[1]+b[1])/2. ]
    return distancia( circuncentros[p[0]] , medio )
  else:
    return distancia( circuncentros[p[0]], circuncentros[p[1]] )
# Se crea una matriz de ceros
estrella1= np.zeros([len(a), len(a)]) ### CHECAR SI SE DEBE QUEDAR COMO MATRIZ DE ENTEROS  d_0= np.zeros([len(a), len(v)], dtype = int)

# Se asignan las orientaciones de las aristas
for i in range(len(a)):
  p=v[a[i][0]] 
  q=v[a[i][1]]
  estrella1[i][i] =(long_dual(a[i]))/(distancia(p, q))

estrella1<jupyter_output><empty_output><jupyter_text>Matriz de derivada exterior para 1 formas ($d_1$)<jupyter_code># Esta es otra manera de calcular d_1

d_1= np.zeros([len(v), len(a)]) 

# Se asignan las orientaciones de las aristas
for i in range(len(v)):
  for j in range(len(a)):
    if ( i==a[j][0] ):
      d_1[i][j]=1
    elif ( i==a[j][1] ):
      d_1[i][j]=-1<jupyter_output><empty_output><jupyter_text>$$d_1^d = -d_0^T$$<jupyter_code>d_1= -np.transpose(d_0)
d_1<jupyter_output><empty_output><jupyter_text>Estrella de Hodge para 2 formas ($*_2$)<jupyter_code># Arreglo vs numeración de los vertices
vs = np.array([i for i in range(len(v))])

# Arreglo vc tal que vc[i] son las caras a las que pertenece el vertice i
vc = [[] for i in vs]
for i in range(len(c)):
    vc[c[i][0]].append(c[i])
    vc[c[i][1]].append(c[i])
    vc[c[i][2]].append(c[i])<jupyter_output><empty_output><jupyter_text>Cálculo del área de la celda dual<jupyter_code>def semiAreaDual(i, ci):
    # Funcion para calcular la aportacion del triangulo ci al area de [v_i]*
    # i es el indice del vertice (dual)
    l = []
    if ci[0] != i:
        l.append(distancia(v[i], v[ci[0]]))
    if ci[1] != i:
        l.append(distancia(v[i], v[ci[1]]))
    if ci[2] != i:    
        l.append(distancia(v[i], v[ci[2]]))
    l = np.array(l)
    l = l/2.0
    
    cci = circuncentro(ci) 
    
    z = distancia(cci, v[i])
    
    area = (l[0]/2.0)*(z**2 - l[0]**2)**0.5+(l[1]/2.0)*(z**2 - l[1]**2)**0.5 # esta formula se deduce sacando el area de los triangulos
    # (v[i], cci, p1), (v[i], cci, p2), donde p1 y p2 son los puntos medios de los segmentos del triangulo ci que contienen a v[i]
    # y utilizando teorema de pitagoras para calcular distancia(p1, cci) y distancia(p2, cci)
    
    return area
def areaDual(i):
    area = 0.0
    for j in range(len(vc[i])):
        area += semiAreaDual(i, vc[i][j])
    
    return area
n = len(vs)
estrella2 = np.zeros((n,n))
for i in range(n):
    estrella2[i][i] = areaDual(i)**(-1)
estrella2<jupyter_output><empty_output><jupyter_text>Matriz del Laplaciano ($\Delta_0$)<jupyter_code>L=np.matmul(estrella2, np.matmul(d_1, np.matmul(estrella1, d_0)))
L<jupyter_output><empty_output><jupyter_text>Solución Podemos agregar las condiciones de frontera $$f(x,y)=10, \quad \forall \quad (x,y)\in \partial \Omega,$$ a nuestro sistema de ecuaciones como sigue. Sea $n$ el número de vértices de la malla asociada a $\Omega$ y $k$ el número de vértices en la frontera. Sea $A$ la matriz de $(n+k)\times n$ tal que las primeras $n\times n$ entradas corresponden a la matriz $[\Delta]$ asociada a la discretización del Laplaciano y la entrada $A_{(n+i)\times j} = 1$ si el vértice $j$ es el $i$-ésimo vértice en la frontera, y $0$ en las demás entradas. Sea $q\in \mathbb{R}^{n+k}$ (como vector columna) tal que las primeras $n$ entradas son iguales a $-1$ y las últimas $k$ entradas son iguales a $10$. De esta manera, al solucionar la ecuación$$ Af = q, \quad f\in \mathbb{R}^n,$$ obtendremos un vector $f$ que representa la función $f_0$ evaluada en los vértices de la malla tal que $\Delta f_0 = -1$ y que cumple las condiciones de frontera.Para solucionar el sistema $Af=q$, dado que $A$ tiene más renglones que columnas, en principio puede que no haya solución (aunque sabemos que sí existe). Quisiéramos de todos modos encontrar $f$ que minimice $||Af-q||^2$, y en el mejor de los casos dicha norma será 0, que es cuando sí hay solución. Tenemos así un problema de mínimos cuadrados lineales el cual sabemos que se soluciona con el sistema $A^TAf = A^Tq$.<jupyter_code>def minimosCuadrados(A, q):
    At = np.transpose(A)
    AtA = At @ A
    Atq = At @ q
    
    f = np.linalg.solve(AtA, Atq)
    return f
def frontierCondition(L, q, frontier):
    # Recibe la discretizacion del Laplaciano L, el vector independiente q, y los indices de los vertices de la malla en la frontera de Omega (la bola de radio 1)
    # Agrega filas a la matriz para considerar las condiciones de frontera y las entradas correspondientes a q
    n = L.shape[0]
    k = len(frontier)
    Ladd = np.zeros([k,n])
    for i in range(k):
        Ladd[i][frontier[i]] = 1.0
    
    L = np.vstack((L, Ladd))
    
    q = np.hstack((q, [10 for i in range(len(frontier))]))
    
    return L, q<jupyter_output><empty_output><jupyter_text>Solución del problema planteado Grafo asociado a la malla<jupyter_code>G = nx.Graph()
G.add_nodes_from(vs)
G.add_edges_from(a)

# Funcion para encontrar los puntos que son frontera
def frontera():
  global G
  f=[]
  for i in G.nodes:
    if (G.degree(i)<=4):
      f.append(i)
  return f

# Vectores necesarios para resolver el problema 
q = np.array([-1 for i in range(len(v))])
frontier = frontera()
# Se aumenta la matriz laplaciana y el vector q
A, qa = frontierCondition(L, q, frontier)

# Se resuelve el sistema
f = minimosCuadrados(A , qa)
f<jupyter_output><empty_output><jupyter_text>Resultados Comparación entre solucion numérica y teóricaSe calculará el error de aproximación entre la solución teórica y la solución encontrada con la discretización.<jupyter_code># Funcion teorica que resuelve el sistema
def solucionTeorica(p):
  return (1-p[0]*p[0]-p[1]*p[1])/4 + 10

# Diferencia entre solucion numerica y teorica
def diferencia(f, v):
  error=0
  for i in range(len(f)):
    error+=abs(f[i]-solucionTeorica(v[i]))
  return error

error=diferencia(f,v)
print("Error de aproximacion: " + str(error) )<jupyter_output>Error de aproximacion: 73.76677068821262<jupyter_text>Solución teórica<jupyter_code>sol_teo = np.array([solucionTeorica(p) for p in v])
sol_teo<jupyter_output><empty_output><jupyter_text>Multiplicación de la matriz Laplaciana multiplicada por el vector que representa a la solución teórica $\Delta f$.<jupyter_code>L @ sol_teo<jupyter_output><empty_output><jupyter_text>Multiplicación de la matriz del Laplaciano por la solución encontrada<jupyter_code>L@f<jupyter_output><empty_output><jupyter_text>GraficaciónPrimero, esta es una visualización de la discretización utilizada para aproximar la función.<jupyter_code>plt.scatter(v[:,0],v[:,1], s=1)<jupyter_output><empty_output><jupyter_text>Ahora se muestra gráficamente la diferencia de la solución teórica con la solución aproximada obtenida.<jupyter_code># Se crea la malla de graficación
X = np.arange(-1, 1, 0.01)
Y = np.arange(-1, 1, 0.01)
X, Y = np.meshgrid(X, Y)

# Se grafica la solución teórica
Z = (1-X**2-Y**2)/4+10
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.set_zlim(9,11)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.viridis, alpha=0.4)

# Se grafican los puntos de la solución aproximada
ax.scatter(v[:,0] , v[:,1] , f,  color='red')

# Gráfica
plt.show()<jupyter_output><empty_output><jupyter_text>La diferencia proyectando sobre el eje  $xz$. En rojo la solución teórica, en azul la aproximada<jupyter_code>X=v[:,0]
x=np.arange(-1, 1, 0.01)
plt.axis([-1,1,8,12])
plt.plot(x, (1-x**2)/4+10, 'r', X, f, 'bd')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15267792954708842, 'fraction_numerical': 0.024442846872753415, 'mean_word_length': 4.754395036194416, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 82, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '530586', 'license': 'MIT', 'n_tokens_mistral': 4749, 'n_tokens_neox': 4629, 'n_words': 1190, 'path': 'EcuacionPoisson/EcuacionPoisson.ipynb', 'repo_name': 'miguelruor/DiscreteDifferencialGeometry', 'size': '11146'}