Tema 4. Aproxime la solución de la EDP elíptica:
\frac{\partial ^2 U}{\partial x^2} + \frac{\partial ^2 U}{\partial y^2} = \frac{x}{y} + \frac{y}{x}1 < x < 2
1 < y < 2
use h = k = 0.5
Tema 3. El área de la sección transversal As (m2) de un lago, a cierta profundidad, se calcula a partir del volumen utilizando la diferenciación:
A_s(Z) = -\frac{\delta V}{\delta z} (Z)Donde V = volumen (m3) y z = profundidad (m), se mide a partir de la superficie en dirección del fondo.
La concentración promedio de una sustancia que varía con la profundidad \overline{c} (g/m3) se obtiene por integración:
\overline{c} = \frac{\int_0^{Z_t} c(Z) A_s(Z) \delta Z}{\int_0^{Z_t}A_s(z) \delta Z}Donde Zt es la profundidad total (m).
Determine la concentración promedio con base en los siguientes datos:
| z (m) | 0 | 4 | 8 | 12 | 16 |
| V (106 m3) | 9.82 | 5.11 | 1.96 | 0.393 | 0.000 |
| c (g/m3) | 10.2 | 8.5 | 7.4 | 5.2 | 4.1 |
zi = [0. , 4 , 8 , 12 , 16] vi = [9.82, 5.11, 1.96, 0.393, 0.] ci = [10.2, 8.5 , 7.4 , 5.2 , 4.1]
Tema 2. Use un método de Runge Kutta para sistemas y aproxime la solución de la siguiente EDO de orden superior
y'''+ 2y'' - y'- 2y = e^t0 ≤ t ≤ 1
y(0) = 1
y'(0) = 0
y''(0) = 0
con h = 0.25
Tema 1. La razón de crecimiento específico g de una levadura que produce un antibiótico es una función de la concentración del alimento c,
g = \frac{2c}{4+0.8c + c^2 +0.2 c^3}
Como se ilustra en la figura, el crecimiento parte de cero a muy bajas concentraciones debido a la limitación de alimento.
También parte de cero en altas concentraciones debido a los efectos de toxicidad.
a) Encuentre el valor de c para el cual el crecimiento es un máximo.
b) Evalúe la función g del problema 1 para c=0,1,2,3, y encuentre el trazador cúbico natural para aproximar el máximo de g, encuentre el error.
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
Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2
\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)por ser los tamaños de paso iguales en ambos ejes, se simplifica la ecuación a usar:
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,
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 2i=1, j=1
i = 2, j =1
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]] >>>
# 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()
Tema 3. (30 puntos) Use el método de diferencias progresivas para aproximar la solución de la siguiente ecuación diferencial parcial parabólica:
\frac{\partial u}{\partial t} - \frac{1}{\pi ^2} \frac{\partial ^2 u}{\partial x^2} =00 ≤ x ≤ 1, t>0
condiciones de borde: u(0,t) = u(1,t) = 0, t>0,
condiciones iniciales u(x,0) = cos(π(x-0.5)), 0 ≤ x ≤ 1
a) Use dx = 0.2 y dt = 0.01. Realize 3 pasos en el tiempo.
b) Estime el error.
c) Calcule la temperatura promedio para t=0 como el área bajo la curva, mediante el método de Simpson. Repita para t=0.01 y calcule el porcentaje que disminuye.
Rúbrica: Construir la malla hasta 5 puntos, plantear la ecuación en el nodo i,j hasta 5 puntos, calcular el estimado de u(i,j) hasta la tercera fila hasta 5 puntos, calcular la temperatura media estimada hasta 5 puntos.
Tema 2. (40 puntos)
a) Usando un polinomio de grado dos obtenga una fórmula central para la primera derivada y otra para la segunda derivada (la tabla tiene al menos 3 nodos, xi-1, xi, xi+1 )
b) Use el método de diferencias finitas para aproximar la solución al siguiente problema con valor de frontera:
y'' = -3y'+2y+2x+30 ≤ x ≤ 1
y(0) = 2
y(1) = 1
use h = 0.25
c) Estime el error
Rúbrica: Plantear un polinomio hasta 5 puntos, deducir la fórmula de la primera derivada hasta 5 puntos, deducir la fórmula de la segunda derivada hasta 5 puntos, plantear el error en las fórmulas hasta 5 puntos. Indicar los nodos en el intervalo hasta 5 puntos, plantear la ecuación en el nodo i hasta 5 puntos, resolver el sistema hasta 5 puntos y estimar el error hasta 5 puntos.
Tema 1. (30 puntos) El movimiento de un sistema acoplado masa-resorte está descrito por la ecuación diferencial ordinaria que sigue:
Donde:
x = el desplazamiento desde la posición de equilibrio (m),
t = tiempo (s),
m = 20 kg masa,
c = 5 (N s/m) coeficiente de amortiguamiento (sub_amortiguado) y
k = 20 (N/m) constante del resorte
La velocidad inicial es cero y el desplazamiento inicial es 1 m.
a) Resuelva esta ecuación con un método numérico para 0<= t <= 15 s, (solo planteo)
b) Realice 3 iteraciones con h=0.1 s
c) Estime el error acumulado en la tercera iteración.
Rúbrica: Plantear el sistema 5 hasta puntos, Plantear el modelo del método numérico hasta 10 puntos, Realizar 3 iteraciones hasta 10 puntos y estimar el error hasta 5 puntos.
Ejercicio: 1Eva_IT2017_T4 Componentes eléctricos
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 x0, x1, x2. 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.
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.
Tema 4. Bono: Deduzca la fórmula de la derivada de y(x1) en función de y(x0), y(x1) y y(x2) considerando h constante.