Ejercicio: 1Eva_2024PAOI_T3 Tasas de natalidad en países desarrollados

Los datos de la tabla que se requieren usar para la tabla son los años 1965, 1980, 1995 y 2010.

Promedio Hijos por mujer en Ecuador. INEC 29 Agosto 2019 [4]

tasa	7.32	4.89	3.50	2.79
año	1965	1980	1995	2010
k	0	3	6	9

literal a

Para cuatro datos o muestras el grado del polinomio que se puede plantear es 3. La ecuación se ordena siguiendo el modelo para la matriz de Vandermonde, primer coeficiente es del termino de mayor grado.

p(x) = ax^3+bx^2+cx+d

Por el orden de magnitud del eje x comparado con los valores de las constantes, se usan los valores del vector k que corresponde a los índices de la tabla. El valor del año se obtiene al hacer un cambio de escala que inicia en el año t₀ = 1965 y tamaño de paso Δx = 5 años.

el sistema de ecuaciones usando cada punto será:

p(0) = a(0)^3+b(0)^2+c(0)+d = 7.32 p(3) = a(3)^3+b(3)^2+c(3)+d = 4.89 p(6) = a(6)^3+b(6)^2+c(6)+d = 3.50 p(9) = a(9)^3+b(9)^2+c(9)+d = 2.79

literal b

El sistema de ecuaciones para la Matriz de Vandermonde D y el vector de constantes B en la forma de matriz aumentada D|B:

\left[ \begin{array}{cccc | c} (0)^3 & (0)^2 & (0) & 1 & 7.32 \\ (3)^3 & (3)^2 & (3) &1 &4.89\\ (6)^3 & (6)^2 & (6) &1 & 3.50 \\ (9)^3 & (9)^2 & (9) & 1 & 2.79 \end{array} \right]

literal c

Resolviendo el sistema de ecuaciones con el algoritmo y usando la instrucción np.linalg.solve(A,B)

>>> np.linalg.solve(A,B)
array([-2.22222222e-03,  7.77777778e-02,
       -1.02333333e+00,  7.32000000e+00])

literal d

p(x) = -0.00222 x^3 + 0.07777 x^2 - 1.0233 x + 7.32

Para el ajuste de escala de años, se usa:

año = x*5 +1965

la gráfica usando el algoritmo, muestra que el polinomio pasa por los puntos seleccionados y existe un error entre los puntos que no participan en la construcción del polinomio.

literal e

El año 2020 en x se escala como

2020 = x*5 +1965

x = 11

o revisando la tabla proporcionada en el ejercicio

p(11) = -0.00222 (11)^3 + 0.07777(11)^2 - 1.0233 (11)x + 7.32 = 2.5166

el error se calcula usando el valor medido para el año 2020

error = | 2.51- 2.35| = 0.16

literal f

p(x) = -0.00222 x^3 + 0.07777 x^2 - 1.0233 x + 7.32 = 2.51

para resolverlo se requiere usar cualquiera de los algoritmos de la unidad 2.

>>> import scipy.optimize as opt
>>> opt.bisect(px,11,30,xtol=0.001)
20.312713623046875
>>>

que corresponde al año = 20.31*5+1965 = 2066.55

Sin embargo la interpolación se usa para estimar valores dentro del intervalo de muestras. El concepto de extrapolación corresponde a otro capítulo. Sin embargo la pregunta se la considera para evaluar la comprensión de la unidad 2.

---

Algoritmo con Python

El resultado con el algoritmo se observa como

xi [0 3 6 9]
fi [7.32 4.89 3.5  2.79]
Matriz Vandermonde: 
[[  0.   0.   0.   1.]
 [ 27.   9.   3.   1.]
 [216.  36.   6.   1.]
 [729.  81.   9.   1.]]
los coeficientes del polinomio: 
[-2.22222222e-03  7.77777778e-02 
 -1.02333333e+00  7.32000000e+00]
Polinomio de interpolación: 
-0.00222222224*x**3 + 0.077777778*x**2 - 1.02333333*x + 7.32

 formato pprint
                       3                      2                            
- 0.002222222224*x  + 0.0777777778*x  - 1.023333333*x + 7.32
>>> px(11)
2.5166666666667066

Instrucciones en Python

# 1Eva_2024PAOI_T3 Tasas de natalidad en países desarrollados
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
tasa = [7.32, 6.39, 5.58, 4.89,4.31,3.85,3.50,3.26,3.03,2.79,2.54,2.35]
anio = [1965, 1970, 1975, 1980,1985,1990,1995,2000,2005,2010,2015,2020]
k    = [   0,     1,   2,    3,   4,   5,   6,   7,   8,   9,  10,  11]
cual = [0,3,6,9]
tasamin = 2.1

# PROCEDIMIENTO
fi = np.take(tasa,cual)
xi = np.take(k,cual)

# El polinomio de interpolación
# muestras = tramos+1
muestras = 25

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi)
B = np.array(fi)
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)
   
# SALIDA
print('xi', xi)
print('fi',fi)
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

# Grafica
plt.plot(anio,tasa,'o')
plt.plot(xi*5+anio[0],fi,'o',color='red', label='[xi,fi]')
plt.plot(xin*5+anio[0],yin,color='red', label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.axhline(0)
plt.axhline(2.1,linestyle='-.')
plt.xlabel('anio')
plt.ylabel('tasa')
plt.title(polinomio)
plt.show()

 

Ejercicio: 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada

literal a

El planteamiento de las ecuaciones según se cita en la solución iterativa es el promedio de los nodos adjacentes:

a = \frac{100+40+b+c}{4} b = \frac{40+20+a+d}{4} c = \frac{100+80+a+d}{4} d = \frac{80+20+c+b}{4}

reordenando las ecuaciones para su forma maticial:

4a -b -c = 100+40 4b -a -d= 40+20 4c -a -d= 100+80 4d -c -b= 80+20

literal b

la forma matricial del sistema de ecuaciones es:

\begin{bmatrix} 4 & -1 &-1 & 0 \\ -1 & 4 & 0 &-1 \\ -1 & 0 & 4 &-1 \\ 0 & -1 &-1 & 4 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 140 \\ 60 \\ 180 \\ 100 \end{bmatrix}

matriz aumentada:

\left[ \begin{array} {cccc|c} 4 & -1 &-1 & 0 & 140\\ -1 & 4 & 0 &-1 & 60 \\ -1 & 0 & 4 &-1 & 180 \\ 0 & -1 &-1 & 4 & 100 \end{array}\right]

al revisar el procedimiento de pivoteo parcial por filas se muestra que ya se encuentra pivoteada la matriz. Por lo que no es necesario realizar cambios en las filas.

literal c

Las expresiones a partir de la matriz aumentada y pivoteada para usar en un método iterativo como Gauss-Seidel son:

a = 0.25(b+c+140) b = 0.25*(a+d+60) c = 0.25(a+d+180) d = 0.25(c+b+100)

el vector inicial debe considerar que las temperaturas serán medias entre los  bordes de la placa

X₀ = [ (100+40)/2, (40+20)/2, (100+80)/2, (80+20)/2 ]

X₀ = [ 70, 30, 90, 50 ]

literal d iteraciones

itera = 0

a = 0.25(30+90+140)= 65

b = 0.25*(65+50+60) = 43.75

c = 0.25(65+50+180) = 73.75

d = 0.25(73.75+43.75+100) = 54.375

error = max|[65-70, 43.75-30, 73.75-90, 54.375-50]|

error = max|[ -5, 13.75, -16.25, 4.375| = 16.25

X₁ = [ 65, 43.75, 73.75, 54.375 ]

itera = 1

a = 0.25(43.75+73.75+140)= 64.375

b = 0.25*(64.375+54.375+60) = 44.687

c = 0.25(64.375+54.375+180) = 74.687

d = 0.25(74.687+44.687+100) = 54.843

error = max|[64.375-65, 44.687-43.75, 74.687-73.75, 54.843-54.375]|

error = max|[-0.625, 0.937, 0.937, 0.468| = 0.937

X₂ = [ 64.375, 44.687, 74.687, 54.843 ]

itera = 2

a = 0.25(44.687+74.687+140)= 64.843

b = 0.25*(64.843 +54.843+60) = 44.921

c = 0.25(64.843+54.843+180) = 74.921

d = 0.25(74.921+44.921+100) = 54.960

error = max|[64.843-64.375, 44.921-44.687, 74.921-74.687, 54.960-54.843]|

error = max|[0.468, 0.234, 0.234, 0.117| = 0.468

X₃ = [ 64.843, 44.921, 74.921, 54.960 ]

literal e

El método converge pues los errores en cada iteración disminuyen. Los valores obtenidos en el resultado son acorde a las temperaturas esperadas en los nodos.

Algoritmos con Python

Los resultados obtenidos son:

X 0 = [65.    43.75  73.75  54.375]
X 1 = [64.375   44.6875  74.6875  54.84375]
X 2 = [64.84375   44.921875  74.921875  54.9609375]
X 3 = [64.9609375  44.98046875 74.98046875 54.99023438]
X 4 = [64.99023438 44.99511719 74.99511719 54.99755859]
X 5 = [64.99755859 44.9987793  74.9987793  54.99938965]
X 6 = [64.99938965 44.99969482 74.99969482 54.99984741]
X 7 = [64.99984741 44.99992371 74.99992371 54.99996185]
X 8 = [64.99996185 44.99998093 74.99998093 54.99999046]
X 9 = [64.99999046 44.99999523 74.99999523 54.99999762]
respuesta X: 
[[64.99999046]
 [44.99999523]
 [74.99999523]
 [54.99999762]]
verificar A.X=B: 
[[139.99997139]
 [ 59.99999285]
 [179.99999285]
 [100.        ]]
>>>

el algoritmo usado es:

# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada
# Método de Gauss-Seidel
import numpy as np
# INGRESO
A = np.array([[4,-1,-1,0],
              [-1,4,0,-1],
              [-1,0,4,-1],
              [0,-1,-1,4]],dtype=float)
B = np.array([140,60,180,100],dtype=float)
X0  = np.array([70,30,90,50],dtype=float)

tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  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('X',itera,'=',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)