s2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Ejercicio: 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Para la ecuación dada con Δx = 1/3, Δt = 0.02, en una revisíón rápida para cumplir la convergencia dt<dx/10, condición que debe verificarse con la expresión obtenida para λ al desarrollar el ejercicio.

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 0 \leq x \leq 2, t>0

literal a. grafica de malla

literal b. Ecuaciones de diferencias divididas a usar

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 \frac{\partial ^2 U}{\partial x^2} = 9 \frac{\partial U}{\partial t} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = 9 \frac{u_{i,j+1}-u_{i,j}}{\Delta t}

se agrupan las constantes,

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

literal d Determine el valor de λ

\lambda = \frac{\Delta t}{9(\Delta x)^2} =\frac{0.02}{9(1/3)^2} = 0.02

valor de λ que es menor que 1/2, por lo que el método converge.

continuando luego con la ecuación general,

\lambda \Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) = u[i,j+1]-u[i,j] \lambda u[i-1,j]-2 \lambda u[i,j] + \lambda u[i+1,j] \Big) = u[i,j+1]-u[i,j]

literal c. Encuentre las ecuaciones considerando las condiciones dadas en el problema.

\lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j] = u[i,j+1]

el punto que no se conoce su valor es u[i,j+1] que es la ecuación buscada.

u[i,j+1] = \lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j]

literal e iteraciones

iteración  i=1, j=0

u[1,1] = \lambda u[0,0]+(1-2 \lambda ) u[1,0] + \lambda u[2,0] u[1,1] =0.02 \cos \Big( \frac{\pi}{2}(0-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}\big(\frac{1}{3}-3\big)\Big) + 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big) \Big) u[1,1] =0.02(0)+(0.96)(-0.5)+0.02(-0.8660)=-0.4973

iteración  i=2, j=0

u[2,1] = \lambda u[1,0]+(1-2 \lambda ) u[2,0] + \lambda u[3,0] u[2,1] = 0.02 \cos \Big( \frac{\pi}{2}(\frac{1}{3}-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}(\frac{2}{3}-3)\Big)+ + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{3}{3}-3\big)\Big) u[2,1] = 0.02 (-0.5) + (0.96 ) (-0.866025) + 0.02 (-1) =-0.8614

iteración  i=3, j=0

u[3,1] = \lambda u[2,0]+(1-2 \lambda ) u[3,0] + \lambda u[4,0] u[3,1] = 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big)\Big)+(1-2 (0.02) ) \cos \Big( \frac{\pi}{2}(1-3)\Big) + + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{4}{3}-3\big)\Big) u[3,1] = 0.02 (-0.866025)+(0.96 ) (-1) + 0.02 (-0,866025) = -0,9946

literal f

la cotas de errores de truncamiento en la ecuación corresponden a segunda derivada O(hx2) y el de primera derivada O(ht), al reemplazar los valores será la suma}

O(hx2) + O(ht) = (1/3)2 + 0.02 = 0,1311

literal g

Resultados usando el algoritmo en Python

Tabla de resultados
[[ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-1.     -0.9946 -0.9893 -0.984  -0.9787 -0.9735 -0.9683 -0.9631 -0.9579  -0.9528]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
T0 = lambda x: np.cos((np.pi/2)*(x-3))
# longitud en x
a = 0.0
b = 2.0
# Constante K
K = 9
# Tamaño de paso
dx = 1/3
dt = 0.02
tramos = int(np.round((b-a)/dx,0))
muestras = tramos + 1
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.linspace(a,b,muestras)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0(xi[1:ultimox])
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

Ejercicio: 2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

La corriente del inductor y(t) para t≥0 se deriva para tener la expresión solo derivadas:

\frac{\delta}{\delta t}y(t) + 2 y(t) + 5 \int_{-\infty}^t y(\tau) \delta \tau = 10 \mu(t)

Para t>0 que es donde transcurre el experimento, el escalón es una constante, se tiene que:

\frac{\delta ^2}{\delta t^2}y(t) + 2 \frac{\delta}{\delta t}y(t) + 5 y(t) = 0

tomando las condiciones iniciales dadas para t=0, y(0)=2, y'(0)=–4

literal a

EL resultadoes perado es el planteamiento del problema. Se reescribe la ecuación con la nomenclatura simplificada y se resordena segun el modelo del método:

y'' = - 2y' - 5 y

luego se sustituye la variable y se convierte a las ecuaciones:

z =y' = f_x(t,y,z) z' = - 2z - 5 y = g_z(t,y,z)

se usa una tabla para llevar el registro de operaciones:

Se plantea las operaciones:

K1y = h * f(ti,yi,zi)
K1z = h * g(ti,yi,zi)

K2y = h * f(ti+h, yi + K1y, zi + K1z)
K2z = h * g(ti+h, yi + K1y, zi + K1z)

yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
ti = ti + h

literal b

El resultado esperado es la aplicación correcta de los valores en las expresiones para al menos tres iteraciones usando h=0.01

itera = 0

K1y = 0.01 y'(0) = 0.01(-4) = -0.04 K1z = 0.01 (- 2z(0) - 5 y(0)) = 0.01(- 2(-4) - 5 (2)) = -0.02 K2y = 0.01 (-4-0.02) = -0.0402 K2z = 0.01 (-2(-4-0.02)-5(2-0.04)) = -0.0176 yi = yi + \frac{K1y+K2y}{2} = 2+\frac{-0.04-0.0402} {2} = 1.9599 zi = zi + \frac{K1z+K2z}{2} = -4 +\frac{-0.02-0.0176}{2} = -4.0188 ti = ti + h = 0+0.01 = 0.01

itera = 1

K1y = 0.01(-4.0188) = -0.040188 K1z = 0.01(- 2(-4.0188) - 5 (1.9599)) = -0.0176 K2y = 0.01 (-4.0188-0.0176) = -0.0403 K2z = 0.01 (-2(-4.0188-0.0176)-5(1.9599-0.040188)) = -0.0152 yi = 1.9599 +\frac{-0.040188-0.0403} {2} = 1.9196 zi = -4.0188 +\frac{-0.0176-0.0152}{2} = -4.0352 ti = ti + h = 0.01+0.01 = 0.02

itera = 2

K1y = 0.01(-4.0352) = -0.040352 K1z = 0.01(- 2(-4.0352) - 5 (1.9196)) = -0.0152 K2y = 0.01 (-4.0352-0.0152) = -0.0405 K2z = 0.01 (-2(-4.0352-0.0152)-5(1.9196-0.040352)) = -0.0129 yi = 1.9196 +\frac{-0.040352-0.0405} {2} =1.8792 zi = -4.0352 +\frac{-0.0152-0.0129}{2} = -4.0494 ti = ti + h = 0.02+0.01 = 0.03

Resultados con el algoritmo en Python

   ti,   yi,    zi,      K1y,    K1z,    K2y,     K2z
[[ 0.00  2.0000 -4.0000  0.0000  0.0000  0.0000   0.0000]
 [ 0.01  1.9599 -4.0188 -0.0400 -0.0200 -0.0402  -0.0176]
 [ 0.02  1.9196 -4.0352 -0.0401 -0.0176 -0.0403  -0.0152]
 [ 0.03  1.8792 -4.0494 -0.0403 -0.0152 -0.0405  -0.0129]
...

Literal c

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

por lo que el error está en el orden de (0.01)3 = 0.000001


Literal d

Se requiere presentar el resultado para el intervalo t entre [0,5]. Siendo el tamaño de paso h=0.01 que es pequeño, se requieren realizar (5-0)/0.01=500 iteraciones, que es más práctico realizarlas usando el algoritmo.

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,7),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -2*z -5*y + 0

t0 = 0
y0 = 2
z0 = -4

h = 0.01
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=4)
print('ti, yi, zi, K1y, K1z, K2y, K2z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y_RK(t)')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('y(t) con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Nota: En el curso TELG1001 Señales y Sistemas, la solución se realiza con Transformadas de Laplace

s2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Ejercicio: 2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Literal a. Integral con Simpson 1/3

Para la ecuación en el intervalo x entre [0,3] aplicando dos veces la fórmula en el intervalo se requieren al menos dos tramos cada Simpson de 1/3. Por lo que la cantidad de tramos es 4 (dos por cada formula, y dos fórmulas) que corresponden a 5 muestras.

El tamaño de paso se calcula como:

h = \frac{b-a}{tramos}=\frac{3-0}{4} = \frac{3}{4} = 0.75

representados en una gráfica como

A = \int_0^3 \frac{e^x \sin(x)}{1+x^2} \delta x

con lo que se define la función del integral f(x)

f(x) = \frac{e^x \sin(x)}{1+x^2}

Con lo que aplicando la fórmula se puede obtener los valores de las muestras:

xi= [0. 0.75   1.5    2.25   3.    ]
fi= [0. 0.9235 1.3755 1.2176 0.2834]

Nota: realizar las expresiones completas para las fórmulas si no adjunta el algoritmo en Python

Aplicando Simpson de 1/3 en cada tramo se tiene:

A_s= \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 0 + 4(0.9235) + 1.3755 \Big) + + \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 1.3755 + 4(1.2176) +0.2834 \Big) A_s = 2.8998

Literal b. Integral con Cuadratura de Gauss

Para usar dos veces el método de Cuadratura de Gauss se usan dos intervalos, con lo que las muestras en x serán:

xj= [0. 1.5 3. ]

se calculan los valores para el tramo [0, 1.5]:

x_{a1} = \frac{0+1.5}{2} - \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 0.3169 x_{b1} = \frac{0+1.5}{2} + \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 1.1830 A_{g1} =\frac{1.5-0}{2} \Big( f(0.3169)+f(1.1830)\Big) = 1.2361

se calculan los valores para el tramo [1.5, 3]:

x_{a2} = \frac{1.5+3}{2} - \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 1.8169 x_{b2} = \frac{1.5+3}{2} + \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 2.6830 A_{g2} =\frac{3-1.5}{2} \Big( f(1.8169)+f(2.6830)\Big) = 1.6329

El total del integral para el intervalo  [0,3]

A_g = A_{g1} + A_{g2} = 2.8691

Al comparar los resultados entre los métodos del literal a y b

errado = |A_s - A_g| = 2.8998 - 2.8691 = 0.0307

Instrucciones integradas en Python

# 2Eva_2022PAOI_T1
# Comparar integrales numéricos Simpson
#   y Cuadratura de Gauss
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: (np.exp(x)*np.sin(x))/(1+x**2)
a = 0
b = 3

# PROCEDIMIENTO
# Aplicando Simpson 1/3
tramos = 4
muestras = tramos+1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
hs = xi[1]-xi[0]
As = (1/3)*(3/4)*(fi[0]+4*fi[1]+fi[2])
As = As + (1/3)*(3/4)*(fi[2]+4*fi[3]+fi[4])
erradoS = 2*(hs**5)

# Aplicando Cuadratura de Gauss
tramosG = 2
muestrasG = tramosG+1
xj = np.linspace(a,b,muestrasG)
hg = xj[1]-xj[0]

xa1 = (xj[0]+xj[1])/2 - (1/np.sqrt(3))*(xj[1]-xj[0])/2
xb1 = (xj[0]+xj[1])/2 + (1/np.sqrt(3))*(xj[1]-xj[0])/2
Ag1 = (hg/2)*(fx(xa1)+fx(xb1))

xa2 = (xj[1]+xj[2])/2 - (1/np.sqrt(3))*(xj[2]-xj[1])/2
xb2 = (xj[1]+xj[2])/2 + (1/np.sqrt(3))*(xj[2]-xj[1])/2
Ag2 = (hg/2)*(fx(xa2)+fx(xb2))

Ag = Ag1 + Ag2

# error entre métodos
errado = np.abs(As-Ag)

# SALIDA
print('xi=',xi)
print('fi=',fi)
print('A Simpson =', As)
print('error Truncamiento Simpson 2*(h^5):',
      erradoS)
print('Cuadratura de Gauss xa,xb por tramos:',
      [xa1,xb1,xa2,xb2])
print('  fx(xa),fx(xb) por tramos:',
      [fx(xa1),fx(xb1),fx(xa2),fx(xb2)])
print('  integral de cada tramo:', [Ag1,Ag2])
print('A Gauss =', Ag)
print('errado entre Simpson y Gauss',errado)

# Grafica con mejor resolucion
xk = np.linspace(a,b,5*tramos+1)
fk = fx(xk)

plt.plot(xk,fk)
plt.plot(xi,fi,'o')
for unx in xi:
    plt.axvline(unx,color='red',
                linestyle='dotted')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('(np.exp(x)*np.sin(x))/(1+x**2)')
plt.grid()
plt.show()

s1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Ejercicio: 1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Día del mes 1 8 15 22
Contagios 1 5.6 27 43.5

a) Realice el planteamiento del sistema de ecuaciones que se usaría usando el método de interpolación polinómica.

El modelo de polinomio de grado máximo que se puede obtener es grado 3:

p_3(t) = a_0 + a_1 t + a_2 t^2 + a_3 t^3

por lo que usando los valores de los puntos dados en la tabla:

p_3(1) = a_0 + a_1 (1) + a_2 (1)^2 + a_3 (1)^3 = 1 p_3(8) = a_0 + a_1 (8) + a_2 (8)^2 + a_3 (8)^3 = 5.6 p_3(15) = a_0 + a_1 (15) + a_2 (15)^2 + a_3 (15)^3 = 27 p_3(22) = a_0 + a_1 (22) + a_2 (22)^2 + a_3 (22)^3 = 43.5

b) Realice el planteamiento del sistema de ecuaciones en su forma matricial y muestre la matriz aumentada.

\begin{pmatrix} 1 & 1 & 1^2 & 1^2\\ 1 & 8 & 8^2 & 8^3 \\ 1 & 15 & 15^2 & 15^3 \\ 1 & 22 & 22^2 & 22^3 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} \begin{pmatrix} 1 \\ 56 \\ 27 \\43.5 \end{pmatrix}

matriz aumentada,

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 8 & 8^2 & 8^3 & 56 \\ 1 & 15 & 15^2 & 15^3 & 27\\ 1 & 22 & 22^2 & 22^3 & 43.5\end{pmatrix}

c) Desarrolle el pivoteo parcial por filas, indicando las operaciones realizadas en éste proceso

pivoteo parcial por filas

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 22 & 22^2 & 22^3 & 43.5 \\ 1 & 15 & 15^2 & 15^3 & 27 \\ 1 & 8 & 8^2 & 8^3 & 56\end{pmatrix}

d) Usando el método directo de Gauss-Jordan, muestre las expresiones necesarias para el algoritmo.

eliminación hacia adelante

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1-1 & 22-1 & 22^2-1^2 & 22^3 -1^3& 43.5 - 1\\ 1-1 & 15-1 & 15^2 -1^2& 15^3 -1^3& 27 -1\\ 1-1 & 8-1 & 8^2 -1^2& 8^3 -1^3& 56-1\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14 & 224 & 3376& 26\\ 0 & 7 & 63& 511 & 55\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14-\frac{14}{21} 21 & 224-\frac{14}{21}483& 3376 -\frac{14}{21}10647& 26-\frac{14}{21}42.5\\ 0 & 7-\frac{7}{21}21 & 63-\frac{7}{21}483& 511-\frac{7}{21}10647 & 55-\frac{7}{21}42.5\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98 & -3038 & 40.83 \end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98-\frac{98}{-98}98 & -3038 -\frac{98}{-98}3722& 40.83 -\frac{98}{-98}(-2.33)\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & 0 & 686 & -7.23\end{pmatrix}

realizando el proceso de eliminación hacia atrás, semejante al método anterior se obtiene

\begin{pmatrix} 1 & 0 & 0 & 0 & 2.98\\ 0 & 1 & 0 & 0& -2.39 \\ 0 & 0 & 1 & 0 & 0.424 \\ 0 & 0 & 0 &1 & -0.0105\end{pmatrix}

con lo que el vector resultado es:

X= [2.98, -2.39, 0.42, -0.0105]

El polinomio de interpolación resultante es:

p(t)= 2.98 -2.39 t + 0.42 t^2 -0.0105 t^3

e) Para el día 19 se encuentra que el valor correspondiente a contagios es de 37%. Estime el error presentado del modelo para ese día.

p(19)= 2.98 -2.39 (19) + 0.42 (19)^2 -0.0105 (19)^3 = 38.42 error = |38.42-37| = 1.42

f) Desarrolle el ejercicio usando otro método para encontrar el polinomio de interpolación.

usando diferencias finitas

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[  0.    1.    1.    4.6  16.8 -21.7   0. ]
 [  1.    8.    5.6  21.4  -4.9   0.    0. ]
 [  2.   15.   27.   16.5   0.    0.    0. ]
 [  3.   22.   43.5   0.    0.    0.    0. ]]

polinomio:

p(t) = 1+\frac{4.6}{1! (7)}(t-1) + + \frac{16.8}{2!(7^2)}(t-1)(t-8) + +\frac{-21.7}{3!(7^3}(t-1)(t-8)(t-15)

Algoritmo en Python

Para literal f

# Polinomio interpolación
# Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([1,8,15,22],dtype=float)
fi = np.array([1,5.6,27,43.5],dtype=float)

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()

s1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos

Ejercicio: 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos

a) Realice el planteamiento del sistema de ecuaciones que permita determinar la cantidad máxima de pacientes de cada grupo que podrían ser atendidos usando todos los productos disponibles. Una vez planteadas las ecuaciones, se le indica que la capacidad de atención para emergencia sea fija en K = 10 pacientes (variable libre).

Producto\ Paciente Maternidad Pos – operatorio Covid_19 emergencia Suministro diario
Producto A 0.2 0.1 1.7 0.25 135
Producto B 0.5 2 0.05 0.4 320
Producto C 1.5 0.2 0.75 1.4 410

Sistema de ecuaciones usando la columna de emergencias como variable libre y valor constante para la variable igual a K=10

0.2 x_1 + 0.1 x_2 +1.7 x_3 = 135-0.25 K 0.5 x_1 + 2 x_2 +0.05 x_3 = 320-0.4 K 1.5 x_1 + 0.2 x_2 +0.75 x_3 = 410-1.4 K

b) Muestre los pasos detallados para la matriz aumentada y pivoteo parcial por filas.

matriz aumentada

\begin{pmatrix} 0.2 & 0.1 & 1.7 & 135-0.25*10 \\ 0.5 & 2 &0.05 & 320-0.4*10 \\ 1.5 & 0.2 & 0.75 &410-1.4*10 \end{pmatrix}

pivoteo parcial por filas

\begin{pmatrix} 1.5 & 0.2 & 0.75 & 396 \\ 0.5 & 2 &0.05 & 316 \\ 0.2 & 0.1 & 1.7 & 132.5 \end{pmatrix}

c) Desarrolle al menos 3 iteraciones para el método requerido, con expresiones completas.
1.5 x_1 + 0.2 x_2 +0.75 x_3 = 396

0.5 x_1 + 2 x_2 +0.05 x_3 = 316 0.2 x_1 + 0.1 x_2 +1.7 x_3 = 132.5

ecuaciones a usar

x_1 = \frac{1}{1.5}(396 - 0.2 x_2 - 0.75 x_3) x_2 = \frac{1}{2}(316 - 0.5 x_1 - 0.05 x_3) x_3 = \frac{1}{1.7}(132.5 - 0.2 x_1 - 0.1 x_2)

Dado que el valor de la variable libre se establecía con K=10, el vector inicial podría ser el doble de este valor

X_0 = [20,20,20]

itera=1

X_0 = [20,20,20] x_1 = \frac{1}{1.5}(396 - 0.2(20) - 0.75(20)) =251.33 x_2 = \frac{1}{2}(316 - 0.5(20) - 0.05 (20)) = 152.5 x_3 = \frac{1}{1.7}(132.5 - 0.2 (20) - 0.1 (20)) =74.41 diferencia = [231.33, 132.5, 54.41] errado = max |[231.33, 132.5, 54.41]| = 231.33

itera=2

X_1 = [251.33, 152.5, 74.41] x_1 = \frac{1}{1.5}(396 - 0.2(152.5) - 0.75(74.41)) =206.46 x_2 = \frac{1}{2}(316 - 0.5(251.33) - 0.05 (74.41)) = 96.30 x_3 = \frac{1}{1.7}(132.5 - 0.2 (251.33) - 0.1 (152.5)) =39.40 diferencia = [-44.86, -59.27, -211.92] errado = max |[-44.86, -59.27, -211.92]| = 211.92

itera=3

X_2 = [231.85, 105.39, 48.16] x_1 = \frac{1}{1.5}(396 - 0.2(105.39) - 0.75(48.16)) =231.85 x_2 = \frac{1}{2}(316 - 0.5(231.85) - 0.05 (48.16)) = 105.39 x_3 = \frac{1}{1.7}(132.5 - 0.2 ( 231.85) - 0.1 (105.39)) =48.16 diferencia = [25.39, 121.09, -158.29] errado = max |[25.39, 121.09, -158.29]| = 158.29 X_3 = [231.85, 105.39, 48.16]

d) Realice las observaciones necesarias sobre los errores entre iteraciones y la convergencia.

El error entre iteraciones disminuye, el método converge a:

respuesta X: 
[[228.22]
 [ 99.81]
 [ 44.92]]

e) Si se decide no atender a los pacientes del grupo emergencias, ¿Qué aumento individual de cada una de otros grupos de pacientes podría soportarse con la cantidad diaria de alimento disponible? (use el algoritmo.py).

Se interpreta como K=0 y usando el algoritmo se obtiene:

respuesta X: 
[[237.23]
 [ 99.54]
 [ 45.64]]

el aumento de pacientes entre grupos es

diferencia = [9.01, -0.27, 0.72]

en términos de pacientes, como número entero, solo se gana 9 pacientes para el primer grupo. Se descarta la parte decimal si los pacientes no se cuentan en miles, por lo que las otras variaciones se interpretan como cero.

Algoritmo en Python

Datos tomados luego de pivoteo parcial por filas

Resultados:

0 [20. 20. 20.]
0 [251.33333333 152.5         74.11764706]
1 [206.60784314  93.31372549  39.10784314]
2 [232.00424837 105.37034314  47.85121107]
3 [226.02501538  98.80265763  44.15418589]
4 [228.74921937 100.38989151  45.24395951]
5 [227.99270138  99.68159617  44.83009822]
6 [228.2940714   99.8810722   44.96076477]
7 [228.20214132  99.80246303  44.91357559]
8 [228.23621714  99.82662528  44.92901496]
9 [228.22527582  99.81772034  44.92358473]
10 [228.22917825  99.82059143  44.92539577]
11 [228.22788993  99.81957054  44.92476777]
12 [228.22834004  99.81990832  44.92497939]
13 [228.2281892   99.8197905   44.92490656]
14 [228.22824132  99.81983004  44.92493124]
numero de condicion: 2.166985328561448
respuesta con Jacobi
[[228.22824132]
 [ 99.81983004]
 [ 44.92493124]]
verificando:
[[396.00002641]
 [316.00002729]
 [132.00001438]]

Algoritmo en Python

# 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos
import numpy as np

def jacobi(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    xnuevo = np.copy(X)

    itera = 0
    print(itera, X)
    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)
        print(itera, xnuevo)
        X = np.copy(xnuevo)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
    return(X)


# INGRESO
A = np.array([[1.5, 0.2, 0.75],
              [0.5, 2.0, 0.05],
              [0.2, 0.1, 1.7 ]],
             dtype=float)

B = np.array([396.,316.,132.],
             dtype=float)
tolera = 1e-4

X = np.array([20.,20.,20.],
             dtype=float)

# PROCEDIMIENTO

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

respuesta = jacobi(A,B,tolera,X)

verifica = np.dot(A,respuesta)

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificando:')
print(verifica)

s1Eva_2022PAOI_T1 Impacto en trayectoria del drone

Ejercicio: 1Eva_2022PAOI_T1 Impacto en trayectoria del drone

Desarrollo analítico

a) Realice el planteamiento del problema usando inicialmente las trayectorias en el eje x, donde para el intervalo de operación del misil antidrone, se observa más de un impacto.

x1(t) = x2(t)

f(t) = cos(t) –  sin(0.75 t) =0

y1(t) = y2(t)

sin(2 t) =kt

k = \frac{sin(2 t)}{t}

b) Usando el método de Newton-Raphson encuentre el valor de t en el cual se pretende realizar el impacto al drone. Realice al menos 3 iteraciones de forma analítica, use tolerancia de 10-4,

f(t) = cos(t) - sin(0.75 t) f'(t) = - sin(t) - 0.75 cos(0.75 t)

Como punto inicial para encontrar la raíz de f(t) podría ser t0=4 para el punto marcado en rojo. Para el método de Newton-Raphson se tiene que

t_{i+1} = t_i - \frac{f(t_i)}{f'(t_i)} error = |t_{i+1} - t_i|

iteración 1 t0=4

t_1 = 4 - \frac{cos(4) - sin(0.75*4)}{- sin(4) - 0.75 cos(0.75*4)} t_1 = 4 - \frac{-0.7947}{1.4992} = 4.5300 tramo = |4.5300-4| = 0.53

iteración 2 t1=4.53

t_{2} = 4.53 - \frac{cos(4.53) - sin(0.75*4.53)}{- sin(4.53) - 0.75 cos(0.75*4.53)} t_2 = 4.53 - \frac{-0.0717}{1.7089} = 4.4880 tramo= |4.4880-4.53| = 0.042

iteración 3 t2=4.4880

t_{3} = 4.4880 - \frac{cos(4.4880) - sin(0.75*4.4880)}{ - sin(4.4880) - 0.75 cos(0.75*4.4880)} t_3 = 4.4880 - \frac{-0.0000179}{1.7061} = 4.4879 tramo= |4.4879-4.4880| = 0.0001

c) Realice el análisis de la convergencia del método.

El error disminuye, el método converge. La raiz se encuentra en t=4.487989505154422

d) Con el resultado de t anterior, determine el valor de la constante k para la expresión de y2(t) que asegura el impacto contra el drone.

y_1(t) = y_2(t) sin(2 t) =kt k = \frac{sin(2 t)}{t} = \frac{sin(2*4.487989505154422)}{4.487989505154422} = 0.096714

Desarrollo con Algoritmo

Resultados

['xi', 'xnuevo', 'tramo']
[[4.0000e+00 4.5301e+00 5.3009e-01]
 [4.5301e+00 4.4880e+00 4.2071e-02]
 [4.4880e+00 4.4880e+00 3.0277e-05]]
raiz en:  4.487989505154422
con error de:  3.0276981949128867e-05

Algoritmo en Python

# Método de Newton-Raphson
import numpy as np

# INGRESO
fx  = lambda t: np.cos(1*t) - np.sin(0.75*t)
dfx = lambda t: -np.sin(t) - np.cos(0.75*t)*0.75

x0 = 4
tolera = 0.0001

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

s2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas

Ejercicio: 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas

literal a

La expresión dada en el enunciado para EDO, se reordena para definir la funcion a usar con Runge-Kutta:

\frac{\delta y(t)}{\delta t} + \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2 = 0 \frac{\delta y(t)}{\delta t} = - \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2

siendo h = 0.5,  con y(0) = 0.15 m y d= 0.01 m ajustando las unidades de medida.

\frac{\delta y(t)}{\delta t} = - \frac{0.01^2}{4}\sqrt{2 (9.8) \text{ }y(t)}\Bigg[\frac{tan (\pi/4)}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t} = - (1.1068e-4) \sqrt{ y(t)}\Bigg[\frac{1}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t}= - (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2}

literal b

se inicia el cálculo del siguiente punto de la tabla

i t y
0 0 0.15
1 0.5 0.1490
2 1 0.1480
3 1.5 0.1471

i = 0

K_1 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2} \Bigg) K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.15}}{0.15^2}\Bigg) = -9.5258e-04 K_2 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)+K_1}}{(y(t)+K_1)^2} \Bigg) K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.15+-9.5258e-04}}{(0.15-9.5258e-04)^2} \Bigg) K_2 = -9.6173e-04 y_1 = y_0 + \frac{K_1 + K_2}{2} y_1 = 0.15 + \frac{-9.5258e-04 -9.6173e-04}{2} = 0.149

i = 1

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.149}}{0.149^2}\Bigg) =-9.6177e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.149-9.7120e-04}}{(0.149-9.7120e-04)^2} \Bigg) K_2 = -9.7116e-04 y_2 = y_1 + \frac{-9.6177e-04 + -9.7116e-04}{2} = 0.1480

i = 2

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.1480}}{0.1480^2}\Bigg) = -9.7120e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.1480-9.7120e-04}}{(0.1480-9.7120e-04)^2} \Bigg) K_2= -9.8084e-04 y_3 = y_2 + \frac{-9.7120e-04 + -9.8084e-04}{2} = 0.1471

literal c

Resultados usando Algoritmo, se encuentra que el embudo se vacia entre 3.15 y 3.20 segundos

 [ t , y , K1 , K2 ]
[[ 0.0000e+00  1.5000e-01  0.0000e+00  0.0000e+00]
 [ 5.0000e-01  1.4904e-01 -9.5258e-04 -9.6173e-04]
 [ 1.0000e+00  1.4808e-01 -9.6177e-04 -9.7116e-04]
 [ 1.5000e+00  1.4710e-01 -9.7120e-04 -9.8084e-04]
 [ 2.0000e+00  1.4611e-01 -9.8088e-04 -9.9078e-04]
 [ 2.5000e+00  1.4512e-01 -9.9083e-04 -1.0010e-03]
...
[ 3.1000e+01  2.8617e-02 -7.5631e-03 -1.0583e-02]
 [ 3.1500e+01  1.0620e-02 -1.1431e-02 -2.4563e-02]
 [ 3.2000e+01         nan -5.0566e-02         nan]
 [ 3.2500e+01         nan         nan         nan]

Instrucciones Python

# 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0,K1,K2]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,K1,K2]
    return(estimado)

# INGRESO
d = 0.01
theta = np.pi/4
g = 9.8
d1y = lambda t,y: -(d**2)/4*np.sqrt(2*g*y)*(np.tan(theta)/y)**2

t0 = 0
y0 = 0.15
h  = 0.5
muestras = 70

# PROCEDIMIENTO
tabla = rungekutta2(d1y,t0,y0,h,muestras)

# SALIDA
np.set_printoptions(precision=4)
print('[ t , y , K1 , K2 ]')
print(tabla)

s2Eva_2021PAOII_T1 Promedio de precipitación de lluvia en área

Ejercicio: 2Eva_2021PAOII_T1 Promedio de precipitación de lluvia en área

Los datos de la tabla corresponden a las medidas de precipitación lluviosa en cada cuadrícula. El área observada es de 300×250, los tramos horizontales 6 y los verticales 4.

Los tamaños de paso son:
h = 300/6 = 50
k = 250/4 = 62.5

f(xi,yj)
i \ j 1 2 3 4 5 6
1 0.02 0.36 0.82 0.65 1.7 1.52
2 3.15 3.57 6.25 5 3.88 1.8
3 0.98 0.98 2.4 1.83 0.04 0.01
4 0.4 0.04 0.03 0.03 0.01 0.08

Para obtener el integral doble, usando la forma compuesta de Simpson, primero desarrolla por filas.

i=1

I_1 = \frac{3}{8}(50) ( 0.02+3(0.36)+3(0.82)+0.65) + \frac{1}{3}(50) (0.65+4(1.7)+1.52) = 228.43

i=2

I_2 = \frac{3}{8}(50) ( 3.15+3(3.57)+3(6.25)+5) + \frac{1}{3}(50) (5+4(3.88)+1.8) = 1077.18

i=3

I_3 = \frac{3}{8}(50) ( 0.98+3(0.98)+3(2.4)+1.83) + \frac{1}{3}(50) (1.83+4(0.04)+0.01) = 276.14

i=4

I_4 = \frac{3}{8}(50) ( 0.4+3(0.04)+3(0.03)+0.03) + \frac{1}{3}(50) (0.03+4(0.01)+0.08) = 14.5

con los resultados, luego se desarrolla por columnas:

columnas = [ 228.43,  1077.18,  276.14,  14.5 ]

I_{total} = \frac{3}{8}(62.5) ( 228.43+3(1077.18)+3(276.14)+14.5) I_{total} = 100850.09

para el promedio se divide el integral total para el área de rectángulo.

Precipitacion_{promedio} = \frac{100850.09}{(300)(250)} = 1.34

Solución con algoritmo

Resultados usando algoritmos con Python:

integral por fila: [ 228.4375     1077.1875      276.14583333   14.5       ]
integral total: 100850.09765625
promedio precipitación:  1.34466796875
>>>

Instrucciones con Python

# 2Eva_2021PAOII_T1 Promedio de precipitación de 
#  lluvia en área
import numpy as np

def simpson_compuesto(vector,h):
    m = len(vector)
    suma = 0
    j = 0
    while (j+3)<m:  # Simpson 3/8
        f = vector[j:j+4]
        s = (3*h/8)*(f[0]+3*f[1]+3*f[2]+f[3])
        suma = suma + s
        j = j + 3
    while (j+2)<m: # Simpson 1/3
        f = vector[j:j+3]
        s = (h/3)*(f[0]+4*f[1]+f[2])
        suma = suma + s
        j = j + 2
    while (j+1)<m: # Trapecio
        f = vector[j:j+2]
        s = (h/2)*(f[0]+f[1])
        suma = suma + s
        j = j + 1
    return(suma)
    

# INGRESO
A = [[0.02, 0.36, 0.82, 0.65, 1.7 , 1.52],
     [3.15, 3.57, 6.25, 5.  , 3.88, 1.8 ],
     [0.98, 0.98, 2.4 , 1.83, 0.04, 0.01],
     [0.4 , 0.04, 0.03, 0.03, 0.01, 0.08]]
base = 300
altura = 250

h = base/6
k = altura/4

# PROCEDIMIENTO
A = np.array(A)
[n,m] = np.shape(A)
columna = np.zeros(n,dtype=float)
for i in range(0,n,1):
    vector = A[i]
    integralfila = simpson_compuesto(vector,h)
    columna[i] = integralfila

total = simpson_compuesto(columna,k)
promedio = total/(base*altura)
    
# SALIDA
print('integral por fila:', columna)
print('integral total:', total)
print('promedio precipitación: ', promedio)

s1Eva_2021PAOII_T3 Nutrientes en asociación de cultivos

Ejercicio: 1Eva_2021PAOII_T3 Nutrientes en asociación de cultivos

a. Realice el planteamiento del sistema de ecuaciones, presente en la forma Ax=B.

Las variables del ejercicio corresponden al consumo de nutrientes:

a = consumo de nutrientes planta de plátano
b = consumo de nutrientes planta de café
c = consumo de nutrientes planta de cacao

que generan las siguientes expresiones de consumo segun la tabla del ejercicio:

\begin{cases} 40a+110b+310c = 750 \\ 400a+15b+25c = 445 \\ 200a+560b+310c = 10 \end{cases}

que en forma A.x=B:

\begin{bmatrix} 40 && 110 && 310 \\ 400 && 15 && 25 \\ 200 && 560 && 310 \end{bmatrix}.\begin{bmatrix} a \\ b \\ c \end{bmatrix}=\begin{bmatrix} 750 \\ 445 \\ 10 \end{bmatrix}

b. De ser necesario, realice operaciones con la matriz aumentada para mejorar la convergencia con un método iterativo.

\begin{bmatrix} 40 && 110 && 310 && 750\\ 400 && 15 && 25 && 445\\ 200 && 560 && 310 &&10\end{bmatrix}

aplicando pivoteo parcial por filas

\begin{bmatrix} 400 && 15 && 25 && 445 \\ 200 && 560 && 310 &&10\\ 40 && 110 && 310 && 750\end{bmatrix}

c. En el contexto del problema, proponga un vector inicial y tolerancia.

el vector inicial no se usa con ceros, un suelo sin nutrientes no es útil para un cultivo, el vector de unos puede considerarse una opción. Otra opción es consultar con los agricultores para disponer de valores iniciales.

[1,1,1]

la toreleancia podría ser de 0.001 si consideramos que los nutrientes estan dados en kilogramos, la tolerancia sería en gramos.

d. Realice 3 iteraciones con el método de Gauss-Seidel y estime el error (papel y lápiz)

las ecuaciones para el método iterativo serán

\begin{bmatrix} 400 && 15 && 25 && 445 \\ 200 && 560 && 310 &&10\\ 40 && 110 && 310 && 750\end{bmatrix} a = \frac{445- 15 b - 25 c}{400} b = \frac{10- 200 a - 310 c}{560} c = \frac{750- 40 a -110 c}{400}

iteración 1

a = \frac{445- 15 (1) - 25 (1)}{400} = 1.0125 b = \frac{10- 200 (1.0125) - 310 (1)}{560}= -0.8973 c = \frac{750- 40 (1.0125) -110 (-0.8973)}{400} = 2.6071 error_a = |1.0125-1| = 0.0125 error_b = |-0.89736-1| = 1.8973 error_c =|2.6071-1| = 1.6071

error_iteración = max| [0.0125, 1.8973, 1.6071]| =1.8973

iteración 2 y 3 se deja como tarea:

Resultados:

respuesta X: 
[[ 0.99999468]
 [-1.99993467]
 [ 2.9999775 ]]
verificar A.X=B: 
[[444.99829063]
 [ 10.02854772]
 [750.        ]]
>>> 

e. Describa y justifique su observación sobre la convergencia del método y estime una descripción de los resultados.
El resultado muestra que el consumo de la planta b (café) es negativo, puede interpretarse como un error de datos por revisar o se iterpreta como si la planta aporta nutrientes al suelo.
Instrucciones en Python

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# INGRESO
A = np.array([[40,110,310],
              [400,15,25],
              [200,560,310]],dtype=float)
B = np.array([[750],[445],[10]],dtype=float)

X0  = np.array([1.,1.,1.],dtype=float)

tolera = 0.001
iteramax = 50

# PROCEDIMIENTO

# Matriz aumentada
AB  = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)
    
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]

A = np.copy(AB[:n,:m])
B = np.copy(AB[:,m])

#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
        print(diferencia,X)
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)