Etiqueta: ejercicios resueltos Python

Ejercicio: 1Eva_IIT2007_T3 Interpolación inversa

f(0.50) = 1.648
f(0.65) = 1.915
f( x  ) = 2.117
f(0.80) = 2.225
f(0.95) = 2.5857

Para el algoritmo se intercambian las variables previo a usarlo.
Luego se evalua en el punto buscado, en éste caso: fi=2.117, obteniendo que x es: 0.750321134121361

fi = np.array([0.50 , 0.65 , 0.80,  0.95   ])
x = np.array([1.648, 1.915, 2.225, 2.5857 ])

Polinomio de Lagrange
0.924124055152463*(-3.74531835205992*x + 7.17228464419475)*(x - 2.5857)*(x - 2.225) + 3.12624749298999*(x - 2.5857)*(x - 2.225)*(3.74531835205992*x - 6.17228464419475) - 7.15454716188057*(x - 2.5857)*(x - 1.915)*(1.73310225303293*x - 2.85615251299827) + 3.92689380344011*(x - 2.225)*(x - 1.915)*(1.06643915964594*x - 1.75749173509651)
Expandiendo: 
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px.subs(x,2.117)
0.750321134121361

El algortmo modificado usado es:

# Interpolacion de Lagrange
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# INGRESO , Datos de prueba
fi = np.array([0.50 , 0.65 , 0.80,  0.95   ])
xi = np.array([1.648, 1.915, 2.225, 2.5857 ])

# PROCEDIMIENTO
n = len(xi)
x = sp.Symbol('x')
# Polinomio
polinomio = 0
for i in range(0,n,1):
    # Termino de Lagrange
    termino = 1
    for j  in range(0,n,1):
        if j!=i:
            termino = termino*(x-xi[j])/(xi[i]-xi[j])
    polinomio = polinomio + termino*fi[i]
    
# Expande el polinomio
px = polinomio.expand()

# Salida
print('Polinomio de Lagrange')
print(polinomio)
print('Expandiendo: ')
print(px)

para visualizar el resultado, intercambiando nuevamente los ejes:

>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226
>>> px.subs(x,2.117)
0.750321134121361
>>> px
0.0358848473081501*x**3 - 0.342756582990887*x**2 + 1.44073214117566*x - 1.10404634485226

para revisar con la gráfica, se añaden las líneas,

# GRAFICA
polinomio = sp.lambdify(x,px)
y = np.linspace(1,3,100)
pyi= polinomio(y)
plt.plot(pyi,y)
plt.show()

Ejercicio: 1Eva_IIT2007_T2 Aplicar Gauss-Seidel

Probando solución con Jacobi, enviado como tarea

Desarrrollo Analítico

– Verificar que la matriz es diagonal dominante

No es necesario realizar el pivoteo por filas, ya la matriz tiene la diagonal dominante.

A = [[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 = [ -9.44, 25.27, -48.01, 19.76, -23.63, 62.59]

– revisar el número de condición

cond(A) = || A||.||A^-1||

El número de condición no es «muy alto»,  los valores de la diagonal son los mayores en toda la fila, por lo que el sistema converge.

>>> np.linalg.cond(A)
2.0451853966291011

– realizar iteraciones

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.]

iteración 1

x0 = (-9.44 -0.30(0) -0.15(0) 
      -0.50(0) -0.34(0) -0.84(0))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(0) -0.70(0) 
      -0.90(0) -0.29(0) -0.57(0))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(0) -0.19(0) -0.82(0)
      -0.34(0) -0.37(0))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(0) -0.68(0) -0.86(0)
      -0.53(0) -0.70(0))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(0) -0.30(0) -0.85(0)
      -0.82(0) -0.55(0))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(0) -0.54(0) -0.59(0)
      -0.66(0) -0.31(0))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.23722149  3.9484375  -5.7635054   1.93535749 -3.97142857  6.76648649]
diferencia = X1 - X0 = X1
errado = max(|diferencia|) = 6.76648649

iteración 2

x0 = (-9.44 -0.30(3.9484375) -0.15(-5.7635054) 
      -0.50(1.93535749) -0.34(-3.97142857) -0.84(6.76648649))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(-1.23722149) -0.70(-5.7635054) 
      -0.90(1.93535749) -0.29(-3.97142857) -0.57(6.76648649))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(-1.23722149) -0.19(3.9484375)
      -0.82(1.93535749) -0.34(-3.97142857) -0.37(6.76648649))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(-1.23722149) -0.68(3.9484375)
      -0.86(-5.7635054) -0.53(-3.97142857) -0.70(6.76648649))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(-1.23722149) -0.30(3.9484375) 
      -0.85(-5.7635054) -0.82(1.93535749) -0.55(6.76648649))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(-1.23722149) -0.54(3.9484375) 
      -0.59(-5.7635054) -0.66(1.93535749) -0.31(-3.97142857))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.97395113  3.95743644 -6.05925771  1.96068604 -4.09171178  6.95612152]
diferencia = X1 - X0 = [-0.73672964,  0.00899894, -0.29575231,  0.02532855, -0.12028321, 0.18963504]
errado = max(|diferencia|) = 0.736729635697

iteración 3

x0 = (-9.44 -0.30(3.95743644) -0.15(-6.05925771) 
      -0.50(1.96068604) -0.34(-4.09171178) -0.84(6.95612152))/7.63 
   = -9.44/7.63 = -1.23722149
x1 = (25.27 -0.38(-1.97395113) -0.70(-6.05925771) 
      -0.90(1.96068604) -0.29(-4.09171178) -0.57(6.95612152))/6.40 
   = 25.27/6.40 = 3.9484375
x2 = (-48.01 -0.83(-1.97395113) -0.19(3.95743644) 
      -0.82(1.96068604) -0.34(-4.09171178) -0.37(6.95612152))/8.33 
   = -48.01/8.33 = -5.7635054
x3 = (19.76 -0.50(-1.97395113) -0.68(3.95743644) 
     -0.86(-6.05925771) -0.53(-4.09171178) -0.70(6.95612152))/10.21 
   = 19.76/10.21 = 1.93535749
x4 = (-23.63 - 0.71(-1.97395113) -0.30(3.95743644) 
      -0.85(-6.05925771) -0.82(1.96068604) -0.55(6.95612152))/5.95 
   = -23.63/5.95 = -3.97142857
x5 = (62.59 - 0.43(-1.97395113) -0.54(3.95743644) 
      -0.59(-6.059257710) -0.66(1.96068604) -0.31(-4.09171178))/9.25 
   = 62.59/9.25 = 6.76648649
X1 = [-1.98566781  4.0185268  -5.9920623   2.01431955 -3.98302285  7.01093224]
diferencia = X1 - X0 = [-0.01171668,  0.06109037,  0.0671954 ,  0.05363351,  0.10868893, 0.05481072]
errado = max(|diferencia|) = 0.108688931048

Desarrollo numérico con Python

Se verifica el resultado obtenido realizando A.Xi y comparando con el vecto B
en la tabla se usa el signo de errado para la gráfica.

X0:  [ 0.  0.  0.  0.  0.  0.]
Xi:  [-1.23722149  3.9484375  -5.7635054   1.93535749 -3.97142857  6.76648649]
errado:  6.76648648649
Xi:  [-1.97395113  3.95743644 -6.05925771  1.96068604 -4.09171178  6.95612152]
errado:  -0.736729635697
Xi:  [-1.98566781  4.0185268  -5.9920623   2.01431955 -3.98302285  7.01093224]
errado:  0.108688931048
Xi:  [-2.00378293  3.99452422 -6.00443878  1.99576482 -4.0067623   6.9961552 ]
errado:  -0.0240025781157
Xi:  [-1.99869529  4.00195452 -5.99863447  2.00153846 -3.99769932  7.00130746]
errado:  0.00906298532211
Xi:  [-2.00045097  3.99933614 -6.00047801  1.99948184 -4.00078219  6.99955127]
errado:  -0.00308287453133
Xi:  [-1.99984629  4.00022733 -5.99983706  2.00017793 -3.99973154  7.00015339]
errado:  0.0010506528253
Xi:  [-2.00005265  3.9999222  -6.00005579  1.99993915 -4.00009178  6.9999475 ]
errado:  -0.0003602430897
Xi:  [-1.99998199  4.00002662 -5.99998091  2.00002082 -3.99996859  7.00001796]
errado:  0.000123195668206
Xi:  [-2.00000616  3.99999089 -6.00000653  1.99999287 -4.00001075  6.99999385]
errado:  -4.21623029112e-05

respuesta de A.X=B : 
[-2.00000616  3.99999089 -6.00000653  1.99999287 -4.00001075  6.99999385]
iteraciones:  10

A.Xi:  [ -9.44006312  25.26992175 -48.01007303  19.75990236 -23.63008584
  62.58992368]
   B: [ -9.44 25.27 -48.01 19.76 -23.63 62.59]

el gráfico de los errores vs iteraciones es:

---

Aplicando Gauss-Seidel

la tabla de aproximaciones sucesivas para el vector X es:

Tabla de iteraciones con AX=B: 
[[ 0.          0.          0.          0.          0.          0.        ]
 [-1.23722149  4.02189753 -5.73196479  2.21089228 -3.51242077  6.91477852]
 [-2.0322951   3.92844105 -6.03202587  1.98957766 -3.97864919  7.00774962]
 [-1.9976784   4.00317297 -6.00049341  1.99807691 -4.00081785  6.99990294]
 [-1.99994191  4.00036665 -5.99978715  2.00000392 -4.00004739  6.99996363]
 [-2.00001274  3.99998231 -5.99999516  2.00000635 -3.99999579  7.00000072]
 [-2.00000008  3.99999833 -6.00000078  1.99999991 -3.99999985  7.00000015]
 [-1.99999994  4.00000007 -6.00000001  1.99999997 -4.00000002  7.        ]]
>>> 

que se obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.
Tarea: X=TX+C

# Algoritmo Gauss-Seidel
# solución de matrices
# métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
#      Rodriguez 5.2 p.162
import numpy as np

def gauss_seidel_tabla(A,B,tolera,X = [0], itermax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    if (len(X)==n):
        X = np.zeros(n, dtype=float)
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    tabla = [np.copy(X)]
        
    itera = 0
    while (errado>tolera or itera>itermax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(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
        tabla.append(np.copy(X))

    tabla = np.array(tabla)
        
    # No converge
    if (itera>itermax):
        X=0
    return(tabla)

# 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

# PROCEDIMIENTO
n = len(A)
X = np.zeros(n, dtype=float)
respuesta = gauss_seidel_tabla(A,B, tolera, X)

# SALIDA
print('Tabla de iteraciones con AX=B: ')
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:

>>> np.linalg.norm(A, np.inf)
13.479999999999999

Tarea: incluir la norma infinito para T

Ejercicio: 1Eva_IIT2007_T1 Distribución binomial acumulada

Dado F=0.4, dado que n=5 y k=1

Con lo que la fórmula para el objetivo se convierte en:

F = \Bigg( \begin{array}{c} 5 \\ 0 \end{array} \Bigg) p ^0 (1-p)^{(5-0)} + \Bigg( \begin{array}{c} 5 \\ 1 \end{array} \Bigg) p ^1 (1-p)^{(5-1)} = 0.4

para buscar la raíz:

f(p) = (1-p)^{5} + 5p (1-p)^{4} - 0.4

Los valores de las combinatorias se calculan como:

>>> import scipy.special as sts
>>> sts.comb(5,0,repetition=False)
1.0
>>> sts.comb(5,1,repetition=False)
5.0
>>> 

como referencia se revisa la gráfica para f(p)

# 1ra Eval II Término 2007
# Tema 1. Distribución binomial acumulada
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sts

fp = lambda p: (1-p)**5 + 5*p*((1-p)**4) - 0.4

a = 0
b = 1
pasos = 100

# PROCEDIMIENTO
xi = np.linspace(a,b,pasos+1)
p_i = fp(xi)

# SALIDA
plt.plot(xi,p_i)
plt.axhline(0)
plt.show()

f(p) = (1-p)^5 + 5p(1-p)^4 - 0.4 = (1-p)^4 (1 - p + 5p) - 0.4 = (1-p)^4 (1 + 4p) - 0.4 = (1-p)^2 (1-p)^2 (1 + 4p) - 0.4 = (1-2p+p^2) (1-2p+p^2) (1 + 4p) - 0.4 = (1 - 4p + 6p^2 - 4p^3 +p^4 ) (1 + 4p) - 0.4 = 1 - 10p^2 + 20p^3 + 15p^4 + 4p^5 - 0.4 = 0.6 - 10p^2 + 20p^3 + 15p^4 + 4p^5

y su derivada:

f'(p) = - 20p + 60p^2 + 60p^3 +20p^4
con lo que se puede desarrollar Newton-Raphson…

Verificando el polinomio a usar con python:

>>> import sympy as sp
>>> p = sp.Symbol('p')
>>> poli = (1-p)**5 + 5*p*((1-p)**4) - 0.4
>>> pol = poli.expand()
>>> pol
4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
>>> pol.diff(p,1)
20*p**4 - 60*p**3 + 60*p**2 - 20*p

Tarea: con lo que puede continuar con el procedimiento para encontrar la raiz…

verificar con raiz: 0.3649852264049102

---

Newton-Raphson

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
import numpy as np

def newtonraphson(funcionx, fxderiva, c, tolera):
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = c - funcionx(c)/fxderiva(c)
        tramo = abs(xnuevo-c)
        c = xnuevo
    return(c)


# PROGRAMA #######################
import matplotlib.pyplot as plt

funcionx = lambda p: 4*p**5 - 15*p**4 + 20*p**3 - 10*p**2 + 0.6
fxderiva = lambda p: 20*p**4 - 60*p**3 + 60*p**2 - 20*p

# INGRESO
c = 0.2
tolera = 0.0000001

# PROCEDIMIENTO
raiz = newtonraphson(funcionx, fxderiva, c, tolera)

# SALIDA
print('raiz en: ', raiz)

el resultado es:

raiz en:  0.3649852264049102
>>> funcionx(0.3649852264049102)
1.4099832412739488e-14
>>>