la derivada se puede obtener usando la librería Sympy, con las instrucciones:
# simplificando el polinomioimport sympy as sym
x = sym.Symbol('x')
p = 2*x**3 - 5*x**2 + 3*x -0.1
Q2 = p/(x-respuesta)
dQ2 = sym.diff(Q2,x,1)
print('Q2')
print(Q2)
print('dQ2')
print(dQ2)
Usando los resultados de Q2(x) y actualizando las expresiones para f(x) y f'(x) se obtienen los siguientes resultados, que verifican las raíces encontradas en el literal a.
Se usa la tabla para plantear los polinomios usando los datos de cada fila. Se evalúa el polinomio con el valor de K=25 (en miles) para completar los datos de la columna K=25
L\K (miles $)
10
20
25
30
40
0
0
0
0
0
10
11.0000
13.0813
14.4768
15.5563
20
18.4997
22.0000
24.3470
26.1626
25
Y
30
25.0746
29.8189
33.0000
35.4608
Para la fila L=0 no es necesario realizar un polinomio, pues se observa que para 0 trabajadores la producción es nula, no existe producción p(K)=0 y pK(25)=0.
Para la fila L=10, se disponen de 4 puntos, por lo que se plantea un polinomio de grado 3:
Se repite el ejercicio del literal a, con el mismo algoritmo para no rescribir todo por las filas por columnas. En el algoritmo se intercambian las filas por columnas y se transpone la matriz.
# INGRESO
...# Se intercambian los valores de fila columna para repetir el algoritmo# obteniendo Y(L=25,k)
M = np.transpose(M)
kj = [0, 10, 20, 30]
li = [10, 20, 30, 40]
l_busca = 25 # trabajadores
k_busca = 25 # inversión en miles de $
Se ajustan las etiquetas de las gráficas y se obtiene el polinomio a usar, que es la producción Y(k) para un capital K dado:
Para el algoritmo de Newton- Raphson, el ejercicio se plantea igualando el polinomio al valor buscado de producción Y=30. Se escribe luego la expresión en como f(x)=0 y f'(x):.
Considerando como x la cantidad de refrigeradoras, y cantidad de cocinas, las ecuaciones a plantear son:
{200x+100y=10002x+1.05y=10.4
La forma matricial del ejercicio se convierte a:
(20021001.05)(xy)=(100010.4)
la matriz aumentada es
(20021001.05100010.4)
Para el pivoteo parcial por filas, dado que el mayor valor de la primera columna se encuentra en la diagonal, no se requiere y la matriz aumentada se mantiene igual.
Para el proceso de eliminación hacia adelante, se incia con el pivote=200
Matriz aumentada
[[ 200. 100. 1000. ]
[ 2. 1.05 10.4 ]]
Pivoteo parcial:
Pivoteo por filas NO requerido
Elimina hacia adelante:
factor: 0.01 para fila: 1
[[2.e+02 1.e+02 1.e+03]
[0.e+00 5.e-02 4.e-01]]
solución:
[1. 8.]
>>>
literal b
Matriz aumentada
[[ 200. 100. 1000. ]
[ 2. 1.1 10.4]]
Pivoteo parcial:
Pivoteo por filas NO requerido
Elimina hacia adelante:
fila 0 pivote: 200
factor: 0.01 para fila: 1
[[2.e+02 1.e+02 1.e+03]
[0.e+00 1.e-01 4.e-01]]
solución:
[3. 4.]
>>>
literal c
Observación: el pequeño cambio de volumen de la cocina no es consistente con los resultados.
El asunto es que la forma de la refrigeradora o cocina no se adapta al volumen disponible, pues son objetos rígidos. Por lo que el sistema de ecuaciones estaría mal planteado.
Las ecuaciones tendrían sentido si esta diseñando el mejor «tamaño» para que entren la mayor cantidad dentro de un container, sin embargo los tamaños de las refrigeradoras y cocinas se encuentran estandarizados.
Revisamos el número de condición, que resulta ser del orden de 104, lo que confirma que el sistema está mal condicionado.
# 1Eva_IIIT2007_T1 Container: Cocinas y Refrigeradorasimport numpy as np
defpivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
if vertabla==True:
print('Matriz aumentada')
print(AB)
print('Pivoteo parcial:')
# Pivoteo por filas AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print(AB)
return(AB)
defgauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
''' Gauss elimina hacia adelante
tarea: verificar términos cero
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
if vertabla==True:
print('Elimina hacia adelante:')
for i inrange(0,n,1):
pivote = AB[i,i]
adelante = i+1
if vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(adelante,n,1):
if (np.abs(pivote)>=casicero):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - factor*AB[i,:]
if vertabla==True:
print(' factor: ',factor,' para fila: ',k)
else:
print(' pivote:', pivote,'en fila:',i,
'genera division para cero')
if vertabla==True:
print(AB)
return(AB)
defgauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
''' Gauss sustituye hacia atras
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Sustitución hacia atras
X = np.zeros(n,dtype=float)
ultfila = n-1
ultcolumna = m-1
for i inrange(ultfila,0-1,-1):
suma = 0
for j inrange(i+1,ultcolumna,1):
suma = suma + AB[i,j]*X[j]
X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
return(X)
# INGRESO
A = [[200, 100 ],
[ 2, 1.05]]
B = [1000, 10.4]
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
AB = gauss_eliminaAdelante(AB,vertabla=True)
X = gauss_sustituyeAtras(AB,vertabla=True)
# SALIDAprint('solución: ')
print(X)
Dado que no se establece en el enunciado el vector inicial, se usará el vector cero. La tolerancia requerida es 10-5
X0 = [ 0. 0. 0. 0. 0. 0.] = [x0,x1,x2,x3,x4,x5]
La tabla obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.
Tarea: X=TX+C
Instrucciones en Python
# 1Eva_IIT2007_T2 Aplicar Gauss-Seidel# Algoritmo Gauss-Seidelimport numpy as np
defpivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
M = np.copy(AB)
# Pivoteo por filas AB
tamano = np.shape(M)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(M[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(M[i,:])
M[i,:] = M[dondemax+i,:]
M[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(pivoteado, 'intercambiar: ',i,dondemax)
if vertabla==Trueand pivoteado==0:
print('Pivoteo por filas NO requerido')
return(M)
defgauss_seidel(A,B,tolera,X0, iteramax=100, vertabla=False, precision=5):
''' Método de Gauss Seidel, tolerancia, vector inicial X0
para mostrar iteraciones: vertabla=True
'''
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
diferencia = 2*tolera*np.ones(n, dtype=float)
errado = np.max(diferencia)
X = np.copy(X0)
itera = 0
if vertabla==True:
print('Iteraciones Gauss-Seidel')
print('itera,[X],[errado]')
np.set_printoptions(precision)
print(itera, X, np.array([errado]))
while (errado>tolera and itera<iteramax):
for i inrange(0,n,1):
xi = B[i]
for j inrange(0,m,1):
if (i!=j):
xi = xi-A[i,j]*X[j]
xi = xi/A[i,i]
diferencia[i] = np.abs(xi-X[i])
X[i] = xi
errado = np.max(diferencia)
itera = itera + 1
if vertabla==True:
print(itera, X, np.array([errado]))
if (itera>iteramax): # No converge
X = itera
print('iteramax superado, No converge')
return(X)
# Programa de prueba ######## INGRESO
A = np.array([[7.63, 0.30, 0.15, 0.50, 0.34, 0.84],
[0.38, 6.40, 0.70, 0.90, 0.29, 0.57],
[0.83, 0.19, 8.33, 0.82, 0.34, 0.37],
[0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
[0.71, 0.30, 0.85, 0.82, 5.95, 0.55],
[0.43, 0.54, 0.59, 0.66, 0.31, 9.25]])
B = np.array([-9.44,25.27,-48.01,19.76,-23.63,62.59])
tolera = 0.00001
X = np.zeros(len(A), dtype=float)
# PROCEDIMIENTO
n = len(A)
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
AB = pivoteafila(A,B, vertabla=True)
A = AB[:,:n]
B = AB[:,n]
respuesta = gauss_seidel(A,B, tolera, X, vertabla=True)
# SALIDAprint('Matriz aumentada y pivoteada:')
print(AB)
print('Vector Xi: ')
print(respuesta)
En el caso de la norma infinito, para la matriz A, se puede usar el algoritmo desarrollado en clase.
Como valor para verificar su algoritmo, se obtuvo:
![Cobb Douglas Y[L,K]](https://blog.espol.edu.ec/analisisnumerico/files/2017/11/CobbDouglas01.png)
