IMT1001-2023: Introduccion a la Ingenieria Matematica
Modulo Analisis Numerico
Prof. Manuel A. Sánchez
In [214]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import pandas as pd
from IPython.display import display, HTML
display(HTML("""<style>.output {display: flex;align-items: center;text-align: center;}</style>"""))
Metodo de diferencias finitas para la ecuacion de Poisson¶
\begin{equation} -\Delta u = f, \quad \text{en } \Omega\subseteq \mathbb R^{d}, \qquad u=g,\quad \text{sobre } \partial \Omega \end{equation}Poisson en dimension $d=1$¶
\begin{equation} - \frac{u_{j-1}-2u_{j} + u_{j+1}}{h^{2}} = f(x_{j}), \quad 1\leq j \leq N-1, \qquad u_0=\alpha,\,u_N = \beta \end{equation}Ejemplo 1:¶
$$ f(x) = 4\pi^{2}\sin(2\pi (x-b)),\quad \alpha = \sin(-2\pi b ),\,\beta = \sin(2\pi(1-b)) $$In [247]:
def example1_1d():
a = 0.0
b = 1.0
c = 0
u = lambda x: np.sin(2.0*np.pi*(x-c))
alpha = np.sin(2.0*np.pi*(a-c)); beta = np.sin(2.0*np.pi*(b-c))
f = lambda x: 4.0*(np.pi)**2*np.sin(2.0*np.pi*(x-c))
return u, (a,b), (alpha, beta), f
In [260]:
def solve_fd1d(domain, bdata, f, N=10):
x = np.linspace(domain[0],domain[1],N+1)
h = (domain[1]-domain[0])/(N)
A = np.diag(2.0*np.ones(N-1)) + np.diag(-1.0*np.ones(N-2),1)+ np.diag(-1.0*np.ones(N-2),-1)
F = f(x[1:N]); F[0]+=bdata[0]/h**2; F[-1]+=bdata[1]/h**2
uh = np.zeros(N+1)
uh[0] = bdata[0]
uh[N] = bdata[1]
uh[1:N] = np.linalg.solve(A,h**2*F)
return x, uh, A
In [294]:
uexact, domain, bdata, f = example1_1d()
x, uh, A1d = solve_fd1d(domain, bdata, f, N=10)
fig, ax = plt.subplots(1,2, figsize=(14,6))
ax[1].spy(A1d, marker='o')
ax[0].plot(x, uh,color='C01', marker='o', linestyle='None', label=r'$u_h$')
xplot = np.linspace(domain[0], domain[1], 100)
ax[0].plot(xplot, uexact(xplot), label=r'$u$')
ax[0].legend()
plt.show()
Ejemplo 1:¶
$$ f(x) = \begin{cases} 100 & \text{si } 0.2\leq x\leq 0.5 \\ -50 & \text{si } 0.6\leq x\leq 0.8 \\ 0 & \text{en otro caso} \end{cases},\quad \alpha = 0,\,\beta = 1 $$In [250]:
def example2_1d():
a = 0.0
b = 1.0
alpha = 0.0; beta = 1.0
def f(x):
if x<=0.5 and x>=0.2:
return 100.0
elif x<=0.8 and x>=0.6:
return -50.0
else:
return 0
def fvec(x):
return np.array([f(xj) for xj in x])
xplot = np.linspace(-1, 2, 100)
return None, (a,b), (alpha, beta), fvec
In [251]:
_, domain, bdata, f = example2_1d()
x, uh, A1d = solve_fd1d(domain, bdata, f, N=600)
fig, ax = plt.subplots(1,2, figsize=(14,6))
xplot = np.linspace(domain[0], domain[1], 100)
ax[1].plot(xplot, f(xplot),color='C00', marker='None', linestyle='-', label=r'$f$')
ax[0].plot(x, uh,color='C01', marker='None', linestyle='--', label=r'$u_h$')
ax[0].legend()
ax[1].legend()
plt.show()
Test de convergencia¶
In [295]:
uexact, domain, bdata, f = example1_1d()
NN = [2**l for l in range(3,12)]; hh = [(domain[1]-domain[0])/(N) for N in NN]
Eh = np.zeros(len(NN), dtype=np.float64); Emax = np.zeros(len(NN), dtype=np.float64)
for n, N in enumerate(NN):
x, uh, _ = solve_fd1d(domain, bdata, f, N=N)
Eh[n] = np.sqrt( hh[n]*(uexact(x)-uh).dot(uexact(x)-uh))
Emax[n] = np.max( np.abs(uexact(x)-uh) )
r = [ np.log(Eh[l+1]/Eh[l])/np.log(hh[l+1]/hh[l]) for l in range(len(Eh)-1)]; r.insert(0,'-')
rmax = [ np.log(Emax[l+1]/Emax[l])/np.log(hh[l+1]/hh[l]) for l in range(len(Emax)-1)]; rmax.insert(0,'-')
df = pd.DataFrame({'h':hh, 'L2 error': Eh, 'rate':r, 'Max error': Emax, 'rate max':rmax})
df
Out[295]:
En dimension $d=2$¶
\begin{equation} - \frac{u_{i+1,j}+u_{i,j+1}-4u_{i,j} + u_{i-1,j}+u_{i,j-1} }{h^{2}} = f_{ij}, \quad 1\leq i,j,\leq N-1 \end{equation}y \begin{equation} u_{i,j} = g(x_{i,j}), \quad \{i,j=0 \wedge i,j=n \} \end{equation}
Ejemplo 1:¶
\begin{align*} f(x,y) & = 8\pi^{2}\sin(2\pi (x+y-c)),\\ g_{b}(x,y) & = \sin(2\pi (x-c) ), \quad \text{en } y = 0 \\ g_{r}(x,y) & = \sin(2\pi(1+y-c)), \quad \text{en } x = 1 \\ g_{t}(x,y) & = \sin(2\pi (x+1-c) ), \quad \text{en } y = 1 \\ g_{l}(x,y) & = \sin(2\pi(y-c)), \quad \text{en } x = 0 \end{align*}In [274]:
def example1_2d():
a = 0.0; b = 1.0
c = 0.1
u = lambda x,y : np.sin(2.0*np.pi*(x+y-c))
f = lambda x,y : 8.0*(np.pi)**2*np.sin(2.0*np.pi*(x+y-c))
def gboundary():
gb = lambda x,y: np.sin(2.0*np.pi*(x+0.0-c))
gr = lambda x,y: np.sin(2.0*np.pi*(1.0+y-c))
gt = lambda x,y: np.sin(2.0*np.pi*(x-0.0-c))
gl = lambda x,y: np.sin(2.0*np.pi*(0+y-c))
g = {'bottom':gb, 'right':gr, 'top': gt, 'left':gl}
return g
return u, (a,b,a,b), gboundary(), f
In [275]:
def rectangulargrid2d(domain,N):
x = np.linspace(domain[0],domain[1],N+1)
y = np.linspace(domain[2],domain[3],N+1)
X, Y = np.meshgrid(x,y)
Nodes = np.zeros( ((N+1)**2,2))
Nodes[:,0] = X.reshape((N+1)**2)
Nodes[:,1] = Y.reshape((N+1)**2)
# boundary nodes
bottom= [*range(0, N+1)]
left = [*range(N+1, (N+1)*N, N+1)]
right = [*range(2*N+1, (N+1)*N, N+1)]
top = [*range((N+1)*N, (N+1)**2)]
allnodes = bottom+left+top+right
boundary = {'bottom':bottom, 'right':right, 'top':top, 'left':left, 'all':allnodes}
return Nodes, boundary, X, Y
In [288]:
uexact, domain, g, f = example1_2d()
X, Y, uh, A2d, Nodes = solve_fd2d(domain, g, f, N=40)
In [289]:
fig = plt.figure(figsize =(14, 9))
ax = fig.add_subplot(111, projection='3d', )
Z = uh.reshape(X.shape)
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, rstride=1, cstride=1)
# Customize the z axis.
ax.set_zlim(-1.01, 1.01)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
fig.colorbar(surf, shrink=0.75, aspect=15)
plt.show()
Test de convergencia¶
In [292]:
uexact, domain, g, f = example1_2d()
NN = [2**l for l in range(3,8)]
hh = [(domain[1]-domain[0])/(N) for N in NN]
Eh = np.zeros(len(NN), dtype=np.float64)
Emax = np.zeros(len(NN), dtype=np.float64)
for n, N in enumerate(NN):
X, Y, uh, _, Nodes = solve_fd2d(domain, g, f, N=N)
Eh[n] = np.sqrt( hh[n]**2*(uexact(Nodes[:,0], Nodes[:,1])-uh).dot(uexact(Nodes[:,0], Nodes[:,1])-uh))
Emax[n] = np.max( np.abs(uexact(Nodes[:,0], Nodes[:,1])-uh) )
r = [ np.log(Eh[l+1]/Eh[l])/np.log(hh[l+1]/hh[l]) for l in range(len(Eh)-1)]; r.insert(0,'-')
rmax = [ np.log(Emax[l+1]/Emax[l])/np.log(hh[l+1]/hh[l]) for l in range(len(Emax)-1)]; rmax.insert(0,'-')
df = pd.DataFrame({'N nodes': (N+1)**2, 'h':hh, 'L2 error': Eh, 'rate':r, 'Max error': Emax, 'rate max':rmax})
df
Out[292]:
In [ ]: