Autor: Edison Del Rosario

Ejercicio: 3Eva_IT2017_T4 EDP elíptica, placa desplazada

La ecuación del problema en forma contínua:

\frac{\delta ^2 U}{\delta x^2} + \frac{\delta ^2 U}{\delta y^2} = \frac{x}{y} + \frac{y}{x}

1 <  x < 2
1 <  y < 2

Se convierte a la versión discreta usando diferencias divididas centradas

---

\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = \frac{x_i}{y_j} + \frac{y_j}{x_i}

---

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy²

\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1

por ser los tamaños de paso iguales en ambos ejes, se simplifica la ecuación a usar:

---

u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

Por simplicidad se usará el método iterativo en el ejercicio, por lo que se despeja la ecuación del centro del rombo formado por los puntos,

---

4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] -\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

u[i,j] = \frac{1}{4}\Big( u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] -\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)\Big)

---

Iteraciones:

Se utiliza una matriz de ceros para la iteración inicial. En el ejercicio se muestran cálculos para 3 nodos, el resto se realiza con el algoritmo en Python.

Para varias iteraciones se usa Δx =Δy = 1/4 = 0.25

y las ecuaciones para los valores en las fronteras o bordes de la placa

U(x,1)= x \ln (x), U(x,2) = x \ln (4x^{2}),1 \lt x \lt 2 U(1,y)= y \ln(y), U(2,y) = 2y \ln (2y), 1 \lt x \lt 2

i=1, j=1

---

u[1,1] = \frac{1}{4}\Big( u[0,1]+u[2,1] + + u[1,0]+u[1,2] -(0.25)^2\Big( \frac{1.25}{1.25} + \frac{1.25}{1.25}\Big)\Big)

---

u[1,1] = \frac{1}{4}\Big(1.25 \ln (1.25)+0 + + 1.25 \ln(1.25) + 0 -(0.25)^2\Big( \frac{1.25}{1.25} + \frac{1.25}{1.25}\Big)\Big)

---

i = 2, j =1

---

u[2,1] = \frac{1}{4}\Big( u[1,1]+u[3,1] + + u[2,0]+u[2,2] -(0.25)^2\Big( \frac{1.5}{1.5} + \frac{1.5}{1.5}\Big)\Big)

---

---

Método iterativo

usando el método iterativo se obtiene los siguientes resultados:

iteraciones:  15
error entre iteraciones:  6.772297286980838e-05
solución para u: 
[[0.         0.27892944 0.60819766 0.97932763 1.38629436]
 [0.27892944 0.69781162 1.1792239  1.7127402  2.29072683]
 [0.60819766 1.1792239  1.8252746  2.53384036 3.29583687]
 [0.97932763 1.7127402  2.53384036 3.42800537 4.38467039]
 [1.38629436 2.29072683 3.29583687 4.38467039 5.54517744]]
>>> 

Algoritmo en Python

# 3Eva_IT2017_T4 EDP elíptica, placa desplazada
# método iterativo
import numpy as np

# INGRESO
# longitud en x
a = 1
b = 2
# longitud en y
c = 1
d = 2
# tamaño de paso
dx = 0.25
dy = 0.25
# funciones en los bordes de la placa
abajo     = lambda x,y: x*np.log(x)
arriba    = lambda x,y: x*np.log(4*(x**2))
izquierda = lambda x,y: y*np.log(y)
derecha   = lambda x,y: 2*y*np.log(2*y)
# función de la ecuación
fxy = lambda x,y: x/y + y/x

# control de iteraciones
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
# tamaño de la matriz
n = int((b-a)/dx)+1
m = int((d-c)/dy)+1
# vectores con valore de ejes
xi = np.linspace(a,b,n)
yj = np.linspace(c,d,m)
# matriz de puntos muestra
u = np.zeros(shape=(n,m),dtype=float)

# valores en los bordes
u[:,0]   = abajo(xi,yj[0])
u[:,m-1] = arriba(xi,yj[m-1])
u[0,:]   = izquierda(xi[0],yj)
u[n-1,:] = derecha(xi[n-1],yj)

# valores interiores
# para menos iteraciones
mitadx = int(n/2)
mitady = int(m/2)
promedio = (u[mitadx,0]+u[mitadx,m-1]+u[0,mitady]+u[n-1,mitady])/4
u[1:n-1,1:m-1] = promedio

# método iterativo
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
    itera = itera +1
    # copia u para calcular errores entre iteraciones
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            # usar fórmula desarrollada para algoritmo
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(dy**2)*fxy(xi[i],yj[j]))/4 
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
print('iteraciones: ',itera)
print('error entre iteraciones: ',erroru)
print('solución para u: ')
print(u)

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

# matrices de ejes para la gráfica 3D
X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()

grafica = Axes3D(figura)
grafica.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()

Ejercicio: 1Eva_IT2017_T4 Componentes eléctricos

Desarrollo Analítico

Solo puede usar las cantidades disponibles de material indicadas, por lo que las cantidades desconocidas de producción por componente se convierten en las incógnitas x₀, x₁, x₂. Se usa el índice cero por compatibilidad con las instrucciones Python.

	Material 1	Material 2	Material 3
Componente 1	5 x0	9 x0	3 x0
Componente 2	7 x1	7 x1	16 x1
Componente 3	9 x2	3 x2	4 x2
Total	945	987	1049

Se plantean las fórmulas a resolver:

5 x0 +  7 x1 + 9 x2 = 945
9 x0 +  7 x1 + 3 x2 = 987
3 x0 + 16 x1 + 4 x2 = 1049

Se reescriben en la forma matricial AX=B

\begin{bmatrix} 5 & 7 & 9\\ 9 & 7 & 3 \\ 3& 16 & 4\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 945 \\ 987 \\ 1049 \end{bmatrix}

Se reordena, pivoteando por filas para tener la matriz diagonalmente dominante:

\begin{bmatrix} 9 & 7 & 3\\ 3 & 16 & 4 \\ 5& 7 & 9\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 987 \\ 1049 \\ 945 \end{bmatrix}

Se determina el número de condición de la matriz, por facilidad con Python:

numero de condicion: 4.396316324708121

Obtenido con:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]])

B = np.array([987.,1049.,945.])
# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

# SALIDA
print('numero de condicion:', ncond)

Recordando que la matriz debe ser tipo real, se añade un punto a los dígitos.

El número de condición es cercano a 1, por lo que el sistema si debería converger a la solución.

Desarrollo con Python

La forma AX=B permite usar los algoritmos desarrollados, obteniendo la solución. Se verifica el resultado al realizar la multiplicación de A con el vector respuesta, debe ser el vector B con un error menor al tolerado.

respuesta con Jacobi
[[62.99996585]
 [44.99998037]
 [34.9999594 ]]
verificando:
[[ 986.99943346]
 [1048.99942111]
 [ 944.99932646]]

Si interpreta el resultado, se debe obtener solo la parte entera [63,45,35] pues las unidades producidas son números enteros.

instrucciones:

# 1Eva_IT2017_T4 Componentes eléctricos
import numpy as np

def jacobi(A,B,tolera,X0,iteramax=100, vertabla=False):
    ''' Método de Jacobi, tolerancia, vector inicial X0
        para mostrar iteraciones: tabla=1
    '''
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)
    xnuevo = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('itera,[X0],errado')
        print(itera, xnuevo, errado)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        X = np.copy(xnuevo)
        itera = itera + 1
        if vertabla==True:
            print(itera, xnuevo, errado)
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)


# INGRESO
A = np.array([[9., 7.,3.],
              [3.,16.,4.],
              [5., 7.,9.]],dtype=float)

B = np.array([987.,1049.,945.],dtype=float)
tolera = 1e-4
X0 = [0.,0.,0.]

# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

respuesta = jacobi(A,B,tolera,X0,vertabla=True)

verifica = np.dot(A,respuesta)
verificaErrado = np.max(np.abs(B-np.transpose(verifica)))

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificar A.X:')
print(verifica)
print('max|A.X - B|')
print(verificaErrado)

al ejecutar el algoritmo se determina que se requieren 83 iteraciones para cumplir con con el valor de error tolerado.