Mascota 50 años

Ejercicio de interpolación

  • Desarrollar por grupos
  • integrar soluciones en un mismo gráfico
  • mostrar polinomios y rangos de existencia

# caparazon superior
xi = [113., 117, 134, 153, 169, 184, 194, 203]
fi = [127., 141, 161, 166, 160, 155, 140, 132]
puntos = [[xi,fi]]

# caparazon inferior 1
xi = [113., 123, 130, 149, 182, 197, 208, 211]
fi = [127., 116, 112, 107, 110, 114, 112, 108]
puntos.append([xi,fi])

# caparazon inferior 2
xi = [107., 114, 120, 143, 170, 192, 210]
fi = [124., 114, 108, 99, 99, 106, 107]
puntos.append([xi,fi])

# Patas 01
xi = [110., 116, 120, 121, 134, 143, 153, 168, 173, 177, 181, 188, 194, 201, 205, 210, 214, 218]
fi = [92., 86, 89, 86, 87, 95, 91, 93, 85, 85, 83, 86, 88, 86, 91, 87, 92, 92]
puntos.append([xi,fi])

# Patas 02
xi = [109., 115, 117, 120, 125, 130, 134, 138, 142, 145, 150, 152, 154]
fi = [91., 87, 82, 82, 78, 78, 81, 78, 76, 79, 77, 80, 86]
puntos.append([xi,fi])

# Patas 03
xi = [172., 175, 182, 186, 191, 195, 201, 205, 210, 211, 217]
fi = [86., 79, 79, 77, 82, 78, 81, 79, 82, 85, 90]
puntos.append([xi,fi])

# Patas 04
xi = [113., 118, 122, 127, 132, 136, 140, 144, 149, 152]
fi = [88., 90, 88, 86, 85, 83, 87, 84, 86, 85]
puntos.append([xi,fi])

# Rabito 01
xi = [97., 102, 108, 111, 117, 121]
fi = [120., 113, 108, 105, 102, 102]
puntos.append([xi,fi])

# cabeza 01
xi = [194., 196, 203, 210, 211, 209]
fi = [177., 167, 157, 149, 138, 135]
puntos.append([xi,fi])

# cabeza 02
xi = [195., 199, 207, 214, 219, 224, 229, 234, 239, 242, 244]
fi = [177., 185, 190, 190, 188, 193, 192, 185, 175, 172, 168]
puntos.append([xi,fi])

# cabeza 03
xi = [220., 226, 234, 239, 241]
fi = [164., 162, 163, 163, 163]
puntos.append([xi,fi])

# cabeza 04
xi = [203., 211, 214, 219, 223, 224, 225, 230, 236, 241]
fi = [115., 119, 125, 124, 127, 137, 146, 149, 154, 162]
puntos.append([xi,fi])

# cabeza 05
xi = [208., 212, 215, 219, 221, 225, 228]
fi = [174., 177, 178, 179, 182, 183, 178]
puntos.append([xi,fi])

# cabeza 06
xi = [206., 210, 214, 218, 220, 224, 229, 233, 232]
fi = [179., 182, 182, 181, 174, 171, 170, 175, 180]
puntos.append([xi,fi])

Mascota descansando

Referencia: Burden 9th Ed. Ejercicio 32 p164

La parte superior de una mascota descansando se describe con los puntos presentados en la tabla.  Los puntos se usan para ejercicios de interpolación.

xi = [1,2,5,6,7,8,10,13,17,20,23,24,25,27,27.7,28,29,30]
fi = [3.0,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5,7.0,6.1,5.6,5.8,5.2,4.1,4.3,4.1,3.0]

3.1 Solución de sistema de 3×3

Solución de un sistema de 3×3 como un punto de intersección de planos

Referencia del ejercicio: http://blog.espol.edu.ec/matg1013/1eva_iit2011_t2-sistema-de-ecuaciones/

Considere el siguiente sistema de ecuaciones:

\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}

Se pueden representar como planos en el espacio despejando la variable z para cada ecuación, de la forma:

\begin{cases} z=(1+ 2x - 5y)/(9)\\z=(6 -7x-y)/(1)\\z=(-26+3x-7y)/(-1)\end{cases}

Se visualizan los resultados de cada ecuación como los planos mostrados en la gráfica:

La intersección de los planos genera un punto cuyas coordenadas corresponden a la solución del sistema.

Para observar mejor del resultado, ejecute las intrucciones en python propuestas y rote el gráfico resultante con el cursor.


Intrucciones en python

Para generalizar el ingreso de las ecuaciones se usa la forma matricial Ax=B, usando solo los coeficientes.

# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

Para visualizar los planos se requiere cada ecuación como funciones(x,y), por simplicidad se usan formulas en formato lambda. Con la forma matricial se generan las ecuaciones de cada plano usando los coeficientes al despejar la variable z.

# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

Para evaluar las funciones se requieren muestras en cada punto del plano X,Y, semejante a las gráficas en 2D que solo requerian muestra en el eje X.

– Primero se generan las muestras para cada eje en los vectores xi, yi.

a = -5
b = 5
muestras = 21

xi = np.linspace(a,b, muestras)
yi = np.linspace(a,b, muestras)

= Las combinaciones entre las muestras de cada eje se obtienen generando las matrices Xi, Yi que representan la malla de muestras en el plano para evaluar cada punto.

Xi, Yi = np.meshgrid(xi,yi)

– Se evaluan los puntos Xi,Yi en cada ecuación, generando las matrices Zi

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

Se grafica cada ecuación como una «malla de alambre» o wireframe, usando librerias 3D, usando Xi,Yi, Z0, luego con Z1 y Z2

El Valor de la solución se grafica como un punto usando «dispersión» o scatter.

El algoritmo completo se muestra a continuación:

# 1ra Evaluación II Término 2011
# Tema 2. Sistema de ecuaciones 3x3
# Desarrollo conceptual como interseccion de Planos

import numpy as np

# INGRESO Ax=B
A = np.array([[-2, 5, 9],
              [ 7, 1, 1],
              [-3, 7,-1]])

B = np.array([1,6,-26])

a = -5
b = 5
muestras = 21

# PROCEDIMIENTO --------
# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]

xi = np.linspace(a,b, muestras)
yi = np.linspace(a,b, muestras)
Xi, Yi = np.meshgrid(xi,yi)

Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)

# solución al sistema
punto = np.linalg.solve(A,B)

# SALIDA
print('respuesta de A.X=B : ')
print(punto)

# GRAFICA de planos
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

figura = plt.figure()
grafica = figura.add_subplot(111, projection='3d')

grafica.plot_wireframe(Xi,Yi,Z0, color ='blue',   label='Ecuación 1')
grafica.plot_wireframe(Xi,Yi,Z1, color ='green',  label='Ecuación 2')
grafica.plot_wireframe(Xi,Yi,Z2, color ='orange', label='Ecuación 3')

grafica.scatter(punto[0],punto[1],punto[2],
                color = 'red', marker='o', label ='punto')

grafica.set_title('Sistema de ecuaciones 3x3')
grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()
plt.show()

4.4 Videos – acople de aeronaves

Para tener un contexto sobre el cual revisar el tema se tiene los siguientes videos:

1. Acoplamiento de aviones para recarca de combustible . www.AiirSource.com. 30/diciembre/2015.
KC-135 Stratotanker in Action – Aircraft Air Refueling

2. Acoplamiento con estación espacial internacional ISS. RT en español . 2/Julio/2010. 
El carguero ruso Progress M-06M pasó de largo la Estación Espacial Internacional fracasado en su intento de acoplarse

4.1 Normas vs distancias en 3D

Normas de un vector en 3D

La norma de un vector se interpreta como una distancias. Las distancias se referencian al origen, o respecto a otra referencia.

Para el caso de soluciones matriciales, lo que se usa en un método iterativo es la diferencia entre una solucion xi vs xi+1 que se usa como criterio del error.

Si se observa como el error entre vectores x1 y xi+1se tendría:

x1     =  [ 1  2  3]
x2     =  [ 2  4 -1]
errado =  [-1 -2  4]

se genera con simples instrucciones de python:

# error como la distancia entre dos puntos
# para el caso en 3D
import numpy as np

x0 = np.array([0,0, 0])
x1 = np.array([1,2, 3])
x2 = np.array([2,4,-1])

errado = x1-x2

print('x1 = ', x1)
print('x2 = ', x2)
print('errado = ', errado)

Sin embargo, para visualizar el concepto se obtiene una gráfica:

La norma de x1 se interpreta como la distancia al origen, es decir:

||x|| = \Big[ \sum_{i=0}^{n} x_i ^2 \Big] ^{1/2} ||x_1|| = \sqrt{x^2+y^2+z^2}
errado =  [-1 -2  4]
||errado|| =  4.58257569495584


Observar en el contexto de los temas en video:

La final, si debe existir acoplamiento, es necesario calcular la raiz cuadrada de los cuadrados?

¿ o solo toman en cuenta las diferencias entre las coordenadas?

¿cuál de las formas tiene menos operaciones?


Norma infinito

||x|| = max\Big[ |x_i| \Big]

Extender el proceso a matrices


# error como la distancia entre dos puntos
# para el caso en 3D
import numpy as np

x0 = np.array([0,0, 0])
x1 = np.array([1,2, 3])
x2 = np.array([2,4,-1])

errado = x1-x2

erradonorma = np.linalg.norm(errado)

print('x1 = ', x1)
print('x2 = ', x2)
print('errado = ', errado)
print('||errado|| = ', erradonorma)

# Grafica
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
figura = plt.figure()
ejes = figura.add_subplot(111,projection = '3d')

punto = x0
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'blue',
             marker='o', label = 'x0')

punto = x1
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'red',
             marker='o', label = 'x1')

punto = x2
[x, y , z] = punto
ejes.scatter(x,y,z, c = 'green',
             marker='o', label = 'x2')

linea = np.concatenate(([x0],[x1]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||x1||')

linea = np.concatenate(([x0],[x2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||x2||')

linea = np.concatenate(([x1],[x2]),axis = 0)
x = linea[:,0]
y = linea[:,1]
z = linea[:,2]
ejes.plot(x,y,z, label = '||error||')

ejes.set_xlabel('eje x')
ejes.set_ylabel('eje y')
ejes.set_zlabel('eje z')
ejes.legend()
plt.show()

1.4.2 Taylor-Polinomio con gráfico animado

Referencia: Burden, Capítulo 1.1 Ejemplo 3.  p11, pdf 21. Chapra, 4.1 p80, pdf104.   Taylor Series (Wikipedia)

El ejercicio se presenta como un complemento para la sección 1.4 para animar la gráfica. Solo como referencia a lo realizado para exponer el tema:

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy

import numpy as np
import sympy as sym

def politaylor(funcionx,x0,n):
    i = 0
    polinomio = 0
    while (i <= n):
        derivada   = funcionx.diff(x,i)
        derivadax0 = derivada.subs(x,x0)
        divisor   = np.math.factorial(i)
        terminoi  = (derivadax0/divisor)*(x-x0)**i
        polinomio = polinomio + terminoi
        i = i + 1
    return(polinomio)

# PROGRAMA  -------------
# Capitulo 1 Ejemplo 2, Burden p11, pdf 21
# INGRESO
x = sym.Symbol('x')
funcionx = sym.cos(x) 

x0 = 0          
n  = 10         # Grado polinomio Taylor
a  = -5 ; b = 5  # x entre [a,b]
muestras = 51

# PROCEDIMIENTO
# tabla polinomios
px_tabla = []
for grado in range(0,n,1):
    polinomio = politaylor(funcionx,x0,grado)
    px_tabla.append(polinomio)

# SALIDA
print('grado :  polinomio')
for grado in range(0,n,1):
    px = px_tabla[grado]
    print(str(grado)+ ' : '+str(px))
    
    print('polinomio: ')
    sym.pprint(px)
    print()


# GRAFICA - TABLA polinomios ------
xi = np.linspace(a,b,muestras)
# Forma lambda, simplifica evaluación
fx = sym.utilities.lambdify(x,funcionx,'numpy')
fi = fx(xi)

# lineas de cada grado de polinomio
px_lineas = np.zeros(shape=(n,muestras), dtype =float)
for grado in range(0,n,1):
    polinomio = px_tabla[grado]
    px = sym.utilities.lambdify(x,polinomio,'numpy')
    px_lineas[grado] = px(xi)


# SALIDA
# GRAFICA CON ANIMACION ------------
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Parametros de trama/foto
narchivo = 'Taylor01' # nombre archivo
retardo = 700   # milisegundos entre tramas
tramas = len(px_lineas)
ymax = 2*np.max(np.abs(fi))

# GRAFICA figura
figura, ejes = plt.subplots()

# Función Base
fx_linea, = ejes.plot(xi,fi,'r')

# Polinomios de tablapoli grado = 0
px_unalinea, = ejes.plot(xi, px_lineas[0],
                         '-.', label='grado: 0')

# Configura gráfica
plt.xlim([a,b])
plt.ylim([-ymax,ymax])
plt.axhline(0, color='k')  # Linea horizontal en cero
plt.title('Polinomio Taylor: '+'f(x) = ' + str(funcionx))
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# cuadros de texto en gráfico
txt_x = (b+a)/2
txt_y = ymax*(1-0.1)
texto_poli = ejes.text(txt_x, txt_y*(1),
                      'p(x):',
                      horizontalalignment='center')
texto_grado = ejes.text(txt_x, txt_y*(1-0.1),
                        'grado:',
                        horizontalalignment='center')

# Nueva Trama
def unatrama(i,xi,pxi):
    
    # actualiza cada linea
    px_unalinea.set_xdata(xi)
    px_unalinea.set_ydata(pxi[i])
    etiquetap = 'p'+str(i)+'(x) = '+str(px_tabla[i])
    px_unalinea.set_label(etiquetap)
    
    # actualiza texto
    texto_poli.set_text(etiquetap)
    texto_grado.set_text('Grado: '+str(i))
    
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numerocolor = i%10
        lineacolor = 'C'+str(numerocolor)
    px_unalinea.set_color(lineacolor)
    
    return (px_unalinea, texto_poli, texto_grado)

# Limpia Trama anterior
def limpiatrama():
    
    px_unalinea.set_ydata(np.ma.array(xi, mask=True))
    px_unalinea.set_label('')
    
    texto_poli.set_text('')
    texto_grado.set_text('')
    
    return (px_unalinea,texto_poli, texto_grado)

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs = (xi,px_lineas),
                              init_func = limpiatrama,
                              interval = retardo,
                              blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_GIFanimado.gif', writer='imagemagick')
# ani.save(narchivo+'_video.mp4')
plt.draw()
plt.show()

9.3 EDP hiperbólicas

Referencia:  Chapra PT8.1 p.860 pdf.884,  Rodriguez 10.4 p.435

Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.

\frac{\delta ^2 u}{\delta t^2}=c^2\frac{\delta ^2 u}{\delta x^2}

Los extremos de la cuerda de longitud unitaria están sujetos.

u[x,t] , 0<x<1, t≥0
u(0,t) = 0 , t≥0
u(1,t) = 0 , t≥0

Al inicio la cuerda está estirada por su punto central:

u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}

Se suelta la cuerda, con velocidad cero desde la posición inicial:

\frac{\delta u(x,0)}{\delta t} = 0

Semejante al procedimiento para EDP parabólicas y elípticas, se discretiza con diferencias finitas divididas:

\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

con error del orden O(\Delta x)^2 + O(\Delta t)^2.
se reagrupa de la forma:

u_{i,j+1}-2u_{i,j}+u_{i,j-1} = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \big( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \big)

el término constante, por facilidad se simplifica con el valor de 1

\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1

si c = 2 y Δx = 0.2, se deduce que Δt = 0.1

que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:

u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}

en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:

u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}

para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero

\frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0 u_{i,-1}=u_{i,1}

que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:

j=0

u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0} u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0} 2 u_{i,1} = u_{i+1,0}+u_{i-1,0} u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}

que se aplica solo cuando j = 0

que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:

Algoritmo en Python:

# Ecuaciones Diferenciales Parciales
# Hiperbólica. Método explicito
import numpy as np

def cuerdainicio(xi):
    n = len(xi)
    y = np.zeros(n,dtype=float)
    for i in range(0,n,1):
        if (xi[i]<=0.5):
            y[i]=-0.5*xi[i]
        else:
            y[i]=0.5*(xi[i]-1)
    return(y)

# INGRESO
x0 = 0
xn = 1 # Longitud de cuerda
y0 = 0
yn = 0 # Puntos de amarre
t0 = 0
c = 2
# discretiza
tramosx = 16
tramost = 32
dx = (xn-x0)/tramosx 
dt = dx/c

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
tj = np.arange(0,tramost*dt+dt,dt)
n = len(xi)
m = len(tj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0] = cuerdainicio(xi)
u[0,:] = y0
u[n-1,:] = yn
# Aplicando condición inicial
j = 0
for i in range(1,n-1,1):
    u[i,j+1] = (u[i+1,j]+u[i-1,j])/2
# Para los otros puntos
for j in range(1,m-1,1):
    for i in range(1,n-1,1):
        u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j]

# SALIDA
np.set_printoptions(precision=2)
print('xi =')
print(xi)
print('tj =')
print(tj)
print('matriz u =')
print(u)

con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.

# GRAFICA
import matplotlib.pyplot as plt
for j in range(0,m,1):
    y = u[:,j]
    plt.plot(xi,y)
plt.title('EDP hiperbólica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# **** GRÁFICO CON ANIMACION ***********
import matplotlib.animation as animation

# Inicializa parametros de trama/foto
retardo = 70   # milisegundos entre tramas
tramas = m
maximoy = np.max(np.abs(u))
figura, ejes = plt.subplots()
plt.xlim([x0,xn])
plt.ylim([-maximoy,maximoy])

# lineas de función y polinomio en gráfico
linea_poli, = ejes.plot(xi,u[:,0], '-')
plt.axhline(0, color='k')  # Eje en x=0
plt.title('EDP hiperbólica')
# plt.legend()
# txt_x = (x0+xn)/2
txt_y = maximoy*(1-0.1)
texto = ejes.text(x0,txt_y,'tiempo:',horizontalalignment='left')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# Nueva Trama
def unatrama(i,xi,u):
    # actualiza cada linea
    linea_poli.set_ydata(u[:,i])
    linea_poli.set_xdata(xi)
    linea_poli.set_label('tiempo linea: '+str(i))
    texto.set_text('tiempo['+str(i)+']')
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    linea_poli.set_color(lineacolor)
    return linea_poli, texto

# Limpia Trama anterior
def limpiatrama():
    linea_poli.set_ydata(np.ma.array(xi, mask=True))
    linea_poli.set_label('')
    texto.set_text('')
    return linea_poli, texto

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs=(xi,u),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo video y GIFAnimado :
# ani.save('EDP_hiperbólica.mp4')
ani.save('EDP_hiperbolica.gif', writer='imagemagick')
plt.draw()
plt.show()

con el algoritmo se cambian las condiciones iniciales de la cuerda y se observan los resultados.
Realizar con otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas.

9.2.2 EDP Elípticas método implícito

con el resultado desarrollado en EDP elípticas para:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{ \delta y^2} = 0

y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

se puede plantear que:

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

con lo que para el método implícito, se plantea un sistema de ecuaciones para determinar los valores en cada punto desconocido.

j=1, i =1

u_{2,1}-4u_{1,1}+u_{0,1} + u_{1,2} +u_{1,0} = 0 u_{2,1}-4u_{1,1}+Ta + u_{1,2} +Tc= 0 -4u_{1,1}+u_{2,1}+u_{1,2} = -(Tc+Ta)

j=1, i =2

u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +u_{2,0} = 0 u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +Tc = 0 u_{1,1}-4u_{2,1}+u_{3,1}+ u_{2,2}= -Tc

j=1, i=3

u_{4,1}-4u_{3,1}+u_{2,1} + u_{3,2} +u_{3,0} = 0 Tb-4u_{3,1}+u_{2,1} + u_{3,2} +Tc = 0 u_{2,1} -4u_{3,1} + u_{3,2} = -(Tc+Tb)

j=2, i=1

u_{2,2}-4u_{1,2}+u_{0,2} + u_{1,3} +u_{1,1} = 0 u_{2,2}-4u_{1,2}+Ta + u_{1,3} +u_{1,1} = 0 -4u_{1,2}+u_{2,2}+u_{1,1}+u_{1,3} = -Ta

j = 2, i = 2

u_{1,2}-4u_{2,2}+u_{3,2} + u_{2,3} +u_{2,1} = 0

j = 2, i = 3

u_{4,2}-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 Tb-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 u_{2,2} -4u_{3,2}+ u_{3,3} +u_{3,1} = -Tb

j=3, i = 1

u_{2,3}-4u_{1,3}+u_{0,3} + u_{1,4} +u_{1,2} = 0 u_{2,3}-4u_{1,3}+Ta + Td +u_{1,2} = 0 -4u_{1,3}+u_{2,3}+u_{1,2} = -(Td+Ta)

j=3, i = 2

u_{3,3}-4u_{2,3}+u_{1,3} + u_{2,4} +u_{2,2} = 0 u_{3,3}-4u_{2,3}+u_{1,3} + Td +u_{2,2} = 0 +u_{1,3} -4u_{2,3}+u_{3,3} +u_{2,2} = -Td

j=3, i=3

u_{4,3}-4u_{3,3}+u_{2,3} + u_{3,4} +u_{3,2} = 0 Tb-4u_{3,3}+u_{2,3} + Td +u_{3,2} = 0 u_{2,3}-4u_{3,3}+u_{3,2} = -(Td+Tb)

con las ecuaciones se arma una matriz:

A = np.array([
    [-4, 1, 0, 1, 0, 0, 0, 0, 0],
    [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
    [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
    [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
    [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
    [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
    [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
    [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
    [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
    ])
B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
                  -Ta,   0,    -Tb,
              -(Td+Ta),-Td,-(Td+Tb)])

que al resolver el sistema de ecuaciones se obtiene:

>>> Xu
array([ 56.43,  55.71,  56.43,  60.  ,  60.  ,  60.  ,  63.57,  64.29,
        63.57])

ingresando los resultados a la matriz u:

xi=
[ 0.   0.5  1.   1.5  2. ]
yj=
[ 0.    0.38  0.75  1.12  1.5 ]
matriz u
[[ 60.    60.    60.    60.    60.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 50.    55.71  60.    64.29  70.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 60.    60.    60.    60.    60.  ]]
>>>

Algoritmo usado para resolver el problema:

# Ecuaciones Diferenciales Parciales
# Elipticas. Método implícito
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 60
Tb = 60
Tc = 50
Td = 70
# dimensiones de la placa
x0 = 0
xn = 2
y0 = 0
yn = 1.5
# discretiza, supone dx=dy
tramosx = 4
tramosy = 4
dx = (xn-x0)/tramosx 
dy = (yn-y0)/tramosy 
maxitera = 100
tolera = 0.0001

A = np.array([
    [-4, 1, 0, 1, 0, 0, 0, 0, 0],
    [ 1,-4, 1, 0, 1, 0, 0, 0, 0],
    [ 0, 1,-4, 0, 0, 1, 0, 0, 0],
    [ 1, 0, 0,-4, 1, 0, 1, 0, 0],
    [ 0, 1, 0, 1,-4, 1, 0, 1, 0],
    [ 0, 0, 1, 0, 1,-4, 0, 0, 1],
    [ 0, 0, 0, 1, 0, 0,-4, 1, 0],
    [ 0, 0, 0, 0, 1, 0, 1,-4, 1],
    [ 0, 0, 0, 0, 0, 1, 0, 1,-4],
    ])
B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb),
              -Ta,0,-Tb,
              -(Td+Ta),-Td,-(Td+Tb)])


# PROCEDIMIENTO
# Resuelve sistema ecuaciones
Xu = np.linalg.solve(A,B)
[nx,mx] = np.shape(A)

xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0] = Tc
u[:,m-1] = Td
u[0,:] = Ta
u[n-1,:] = Tb
u[1:1+3,1] = Xu[0:0+3]
u[1:1+3,2] = Xu[3:3+3]
u[1:1+3,3] = Xu[6:6+3]

# SALIDA
np.set_printoptions(precision=2)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

La gráfica de resultados se obtiene de forma semejante al ejercicio con método iterativo.

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()
ax = Axes3D(figura)
ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)
plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Se podría estandarizar un poco más el proceso para que sea realizado por el algoritmo y sea más sencillo generar la matriz con más puntos. Tarea.

9.2.1 EDP elípticas método iterativo

con el resultado desarrollado en EDP elípticas para:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{ \delta y^2} = 0

y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1

se puede plantear que:

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

que reordenando para un punto central desconocido se convierte a:

u_{i,j} = \frac{1}{4} \big[ u_{i+1,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} \big]

con lo que se interpreta que cada punto central es el resultado del promedio de los puntos alrededor del rombo formado en la malla.

El cálculo numerico se puede realizar de forma iterativa haciendo varias pasadas en la matriz, promediando cada punto. Para revisar las iteraciones se controla la convergencia junto con un máximo de iteraciones.

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 60
Tb = 60
Tc = 50
Td = 70
# dimensiones de la placa
x0 = 0
xn = 2
y0 = 0
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25 
dy = 0.25 
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz u
u = np.zeros(shape=(n,m),dtype = float)
# valores en fronteras
u[0,:] = Ta
u[n-1,:] = Tb
u[:,0] = Tc
u[:,m-1] = Td

# valor inicial de iteración
promedio = (Ta+Tb+Tc+Td)/4
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
while not(itera>=maxitera and converge==1):
    itera = itera +1
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
np.set_printoptions(precision=2)
print('converge = ', converge)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

Un ejemplo de resultados:

converge = 1
xi=
[ 0.    0.25  0.5   0.75  1.    1.25  1.5   1.75  2.  ]
yi=
[ 0.    0.25  0.5   0.75  1.    1.25  1.5 ]
matriz u
[[ 50.    60.    60.    60.    60.    60.    70.  ]
 [ 50.    55.6   58.23  60.    61.77  64.4   70.  ]
 [ 50.    54.15  57.34  60.    62.66  65.85  70.  ]
 [ 50.    53.67  56.97  60.    63.03  66.33  70.  ]
 [ 50.    53.55  56.87  60.    63.13  66.45  70.  ]
 [ 50.    53.67  56.97  60.    63.03  66.33  70.  ]
 [ 50.    54.15  57.34  60.    62.66  65.85  70.  ]
 [ 50.    55.6   58.23  60.    61.77  64.4   70.  ]
 [ 50.    60.    60.    60.    60.    60.    70.  ]]
>>>

Cuyos valores se interpretan mejor en una gráfica, en este caso 3D:

La gráfica de resultados requiere ajuste de ejes, pues el índice de filas es el eje x, y las columnas es el eje y. La matriz graficada se obtiene como la transpuesta de u

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()
ax = Axes3D(figura)
ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)
plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()