Dado el vector inicial X0= [0, 0, 0], y una tolerancia de 0.0001, e desarrollan al menos 3 iteraciones:
El error se verifica para éste ejercicio como el mayor valor de la diferencia entre iteraciones consecutivas. Analogía al video del acople entre las aeronaves. Si una coordenada aún no es la correcta …. menor a lo tolerado, pues NO hay acople…
La gráfica muestra las coordenadas de las aproximaciones, observe el espiral que se va cerrando desde el punto inicial en verde al punto final en rojo.
Se debe controlar el número de iteraciones para verificar la convergencia con iteramax.
NO es necesario almacenar los valores de los puntos, solo el último valor y el contador itera permite determinar si el sistema converge.
# 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante# MATG1052 Métodos Numéricos# propuesta: edelros@espol.edu.ec,import numpy as np
defjacobi(A,B,X0, tolera, iteramax=100, vertabla=False):
''' Método de Jacobi, tolerancia, vector inicial X0
para mostrar iteraciones: vertabla=True
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
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)
tabla = [np.copy(X0)]
itera = 0
if vertabla==True:
print('Iteraciones Jacobi')
print('itera,[X],errado')
print(itera, xnuevo, errado)
whilenot(errado<=tolera or itera>iteramax):
for i inrange(0,n,1):
nuevo = B[i]
for j inrange(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)
tabla = np.concatenate((tabla,[X]),axis = 0)
itera = itera + 1
if vertabla==True:
print(itera, xnuevo, errado)
# No convergeif (itera>iteramax):
X=itera
print('iteramax superado, No converge')
return(X,tabla)
defpivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas, entrega matriz aumentada AB
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+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print('AB')
return(AB)
# PROGRAMA Búsqueda de solucion# INGRESO --------# INGRESO
A = [[-2, 5, 9],
[ 7, 1, 1],
[-3, 7,-1]]
B = [1, 6,-26]
X0 = [0, 0, 0]
tolera = 0.0001
iteramax = 100
# PROCEDIMIENTO --------
AB = pivoteafila(A,B,vertabla=True)
# separa matriz aumentada en A y B
[n,m] = np.shape(AB)
A = AB[:,:m-1]
B = AB[:,m-1]
[X, puntos] = jacobi(A,B,X0,tolera,vertabla=True)
iterado = len(puntos)
# SALIDA --------print('respuesta de A.X=B : ')
print(X)
print('iterado, incluyendo X0: ', iterado)
# GRAFICA DE ITERACIONESimport matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
figura = plt.figure() # Una hoja para dibujar
grafica = figura.add_subplot(111,projection = '3d')
# Puntos de la iteración
pxi = puntos[:,0]
pyi = puntos[:,1]
pzi = puntos[:,2]
grafica.plot(pxi,pyi,pzi, color = 'magenta',
label = 'puntos[i]')
grafica.scatter(pxi[0],pyi[0],pzi[0], color = 'green', marker='o')
grafica.scatter(pxi[iterado-1],pyi[iterado-1],pzi[iterado-1],
color = 'red', marker='o')
grafica.set_title('Puntos de iteración')
grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()
plt.show()
Gráfica de planos
Dado que el sistema de ecuaciones es de tres incógnitas, la solución se puede interpretar como la intersección de los planos formados por cada una de las ecuaciones en el espacio X,Y,Z.
se comenta la última instrucción del algoritmo anterior, #plt.show(), y se añaden las siguientes líneas al algoritmo anterior para obtener la gráfica:
# Tema 2. Sistema de ecuaciones 3x3# Concepto como interseccion de Planos##import matplotlib.pyplot as plt##from mpl_toolkits.mplot3d import Axes3D
ax = -5 # Intervalo X
bx = 5
ay = ax-2 # Intervalo Y
by = bx+2
muestras = 11
# 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(ax,bx, muestras)
yi = np.linspace(ay,by, 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)
# SALIDAprint('respuesta de A.X=B : ')
print(punto)
# Interseccion entre ecuación 1 y 2# PlanoXZ, extremo inferior de y
Aa = np.copy(A[0:2,[0,2]])
Ba = np.copy(B[0:2])
Ba = Ba-ay*A[0:2,1]
pta = np.linalg.solve(Aa,Ba)
pa = np.array([ay])
pxza = np.array([pta[0],ay,pta[1]])
# PlanoXZ, extremo superior de y
Bb = np.copy(B[0:2])
Bb = Bb-by*A[0:2,1]
ptb = np.linalg.solve(Aa,Bb)
pb = np.array([by])
pxzb = np.array([ptb[0],by,ptb[1]])
# GRAFICA de planos
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111, projection='3d')
graf3D.plot_wireframe(Xi,Yi,Z0,
color ='blue',
label='Z1')
graf3D.plot_wireframe(Xi,Yi,Z1,
color ='green',
label='Z2')
graf3D.plot_wireframe(Xi,Yi,Z2,
color ='orange',
label='Z3')
# recta intersección planos 1 y 2
graf3D.plot([pxza[0],pxzb[0]],
[pxza[1],pxzb[1]],
[pxza[2],pxzb[2]],
label='Sol 1y2',
color = 'violet',
linewidth = 4)
# Punto solución del sistema 3x3
graf3D.scatter(punto[0],punto[1],punto[2],
color = 'red',
marker='o',
label ='punto',
linewidth = 6)
graf3D.set_title('Sistema de ecuaciones 3x3')
graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_xlim(ax, bx)
graf3D.set_xlim(ay, by)
graf3D.legend()
graf3D.view_init(45, 45)
# rotacion de ejesfor angulo inrange(45, 360+45, 5 ):
graf3D.view_init(45, angulo)
plt.draw()
plt.pause(.001)
plt.show()
Tarea: Realice las modificaciones necesarias para usar el algoritmo de Gauss-Seidel. Luego compare resultados.
Revisión de resultados
Si se usan las ecuaciones sin pivorear y cambiar a diagonal dominante, como fueron presentadas en el enunciado, el algoritmo Jacobi NO converge, el error aumenta, y muy rápido como para observar la espiral hacia afuera en una gráfica.
con lo que se pueden obtener cada precio unitario en función de k.
Como variante, se continua siguiendo el procedimieno de Gauss, dejando como tarea el uso de Gauss-Jordan
En el literal c se indica que el valor de k es 65, con lo que se requiere sustituir en la solución el valor de K para encontrar los precios unitarios.
⎣⎡532395184176535⎦⎤
Se encuentra que:
el vector solución X es:
[[-0.18181818]
[ 5.18181818]
[ 2.36363636]]
Lo que muestra que debe existir un error en el planteamiento del enunciado, considerando que los precios NO deberían ser negativos como sucede con el primer precio unitario de la respuesta.
que es lo que suponemos ser trata de corregir en el literal d, al indicar que se cambie en la matriz el valor de 5 por 5.1. Los resultados en éste caso son más coherentes con el enunciado. Todas las soluciones son positivas.
El error relativo de los precios encontrados entre las ecuaciones planteadas es:
diferencia = [0.335968-0.181818,
3.886693-5.181818,
3.626482-2.363636]
= [0.154150, -1.295125, 1.262845]
error(dolares) = max|diferencia| = 1.295125
Por las magnitudes de los precios, el error se aprecia
usando el error relativo referenciado
al mayor valor de la nueva solución
error relativo = 1.295125/3.886693 = 0.333220
es decir de aproximadamente 33%
Para revisar otra causa del error se analiza el número de condición de la matriz:
El número de condición resulta lejano a 1, por lo que para éste problema:
pequeños cambios en la matriz de entrada producen grandes cambios en los resultados.
por ejemplo: un 0.1/5= 0.02 que es un 2% de variación en la entrada produce un cambio del 33% en el resultado.
# 1Eva_IT2011_T3_MN Precios unitarios en factura, k# Método de Gauss-Jordan# Solución a Sistemas de Ecuaciones# de la forma A.X=Bimport 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+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print(AB)
return(AB)
defgauss_eliminaAdelante(AB, vertabla=False,lu=False,casicero = 1e-15):
''' Gauss elimina hacia adelante, a partir de,
matriz aumentada y pivoteada.
Para respuesta en forma A=L.U usar lu=True entrega[AB,L,U]
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
L = np.identity(n,dtype=float) # Inicializa Lif 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,:]
L[k,i] = factor # llena Lif vertabla==True:
print(' factor: ',factor,' para fila: ',k)
else:
print(' pivote:', pivote,'en fila:',i,
'genera division para cero')
respuesta = AB
if vertabla==True:
print(AB)
if lu==True:
U = AB[:,:n-1]
respuesta = [AB,L,U]
return(respuesta)
defgauss_eliminaAtras(AB, vertabla=False, precision=5, casicero = 1e-15):
''' Gauss-Jordan elimina hacia atras
Requiere la matriz triangular inferior
Tarea: Verificar que sea triangular inferior
'''
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
ultfila = n-1
ultcolumna = m-1
if vertabla==True:
print('Elimina hacia Atras:')
for i inrange(ultfila,0-1,-1):
pivote = AB[i,i]
atras = i-1 # arriba de la fila iif vertabla==True:
print(' fila',i,'pivote: ', pivote)
for k inrange(atras,0-1,-1):
if (np.abs(AB[k,i])>=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')
AB[i,:] = AB[i,:]/AB[i,i] # diagonal a unos
X = np.copy(AB[:,ultcolumna])
if vertabla==True:
print(AB)
return(X)
# PROGRAMA ------------------------# INGRESO
A = np.array([[2,5,4],
[3,9,8],
[5,3,1]])
B = np.array([[35],
[65],
[17]])
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
AB = gauss_eliminaAdelante(AB,vertabla=True)
X = gauss_eliminaAtras(AB,vertabla=True)
# SALIDAprint('solución X: ')
print(X)
Se observa la gráfica de la función que muestra el valor máximo en el intervalo entre [2, 5] , alrededor de 3
Para encontrar el valor del máximo, se requiere usar la derivada respecto al tiempo y encontrar el cruce por cero.
dtdC=−te−t/3+e−t/3
Adicionalmente, para el literal b en la gráfica se incluye como meta donde C(t) disminuye a 1/4 o 0.25. Lo que estaría en el intervalo entre [10, 14].
Instrucciones en Python para observar la función:
# 1ra Evaluación I Término 2011# Tema 1. Fondo de Inversionimport numpy as np
import matplotlib.pyplot as plt
ft = lambda t: t*np.exp(-t/3)
a=0
b=20
tolera = 0.0001
muestras = 101
meta = 0.25
# PROCEDIMIENTO# Observar la función entre [a,b]
ti = np.linspace(a,b,muestras)
fti = ft(ti)
# Salida# Gráfica
plt.plot(ti,fti, label='f(t)')
plt.axhline(meta, color = 'r')
plt.axhline(0, color = 'k')
plt.legend()
plt.show()
literal a
Para encontrar el máximo se puede determinar la derivada con Sympy.
Para la derivada se usa la forma simbólica de la función, que se convierte a forma numérica lambda para evaluarla de forma más fácil y obtener la gráfica.
# Literal a) usando derivada simbólicaimport sympy as sp
x = sp.Symbol('x')
fxs = x*sp.exp(-x/3)
dfxs = fxs.diff(x,1)
# convierte la expresión a lambda
dfxn = sp.utilities.lambdify(x,dfxs,'numpy')
dfxni = dfxn(ti)
print('derivada de la función: ')
print(dfxs)
# Gráfica de la derivada.
plt.plot(ti,dfxni, label='df(t)/dt')
plt.axhline(0, color = 'k')
plt.legend()
plt.show()
derivada de la función:
-x*exp(-x/3)/3 + exp(-x/3)
Se busca la raíz con algún método, por ejemplo bisección, siendo f(x) = 0
# 1Eva_IT2011_T1_MN Fondo de Inversión# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
defbiseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
'''
Algoritmo de Bisección
Los valores de [a,b] son seleccionados
desde la gráfica de la función
error = tolera
'''
fa = fx(a)
fb = fx(b)
tramo = np.abs(b-a)
itera = 0
cambia = np.sign(fa)*np.sign(fb)
if cambia<0: # existe cambio de signo f(a) vs f(b)if vertabla==True:
print('método de Bisección')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
np.set_printoptions(precision)
while (tramo>=tolera and itera<=iteramax):
c = (a+b)/2
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if vertabla==True:
print(itera,np.array([a,c,b]),np.array([fa,fc,fb]))
if (cambia<0):
b = c
fb = fc
else:
a = c
fa = fc
tramo = np.abs(b-a)
if vertabla==True:
print(' ',tramo)
itera = itera + 1
respuesta = c
# Valida respuestaif (itera>=iteramax):
respuesta = np.nan
else:
print(' No existe cambio de signo entre f(a) y f(b)')
print(' f(a) =',fa,', f(b) =',fb)
respuesta=np.nan
return(respuesta)
# INGRESO
fx = lambda t: -t*np.exp(-t/3)/3 + np.exp(-t/3)
a = 2
b = 5
tolera = 0.0001
# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDAprint('raíz en: ', respuesta)
Determine el monto de la inversión inicial A necesaria para que el máximo sea igual a un millón de dólares. Como el máximo se encuentra en t=3, se tiene que:
C(t)=Ate−t/31=A(3)e−3/3A=(3)e−11=0.906093
considerando que las unidades se encuentran en millones.
literal b
Para el literal b, se busca la raíz usando el método de Newton-Raphson como se indica en el enunciado.
En la función nueva se usa el valor de A encontrado, y la meta establecida.
Se obtiene la misa derivada del problema anterior multiplicada por A, por ser solo un factor que multiplica a la función original. El valor de meta es una constante, que se convierte en cero al derivar.
