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)

 

Ejercicio: 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico

literal a. Planteamiento

A partir de la solución expresión propuesta y simplificada del ejercicio, se reemplaza los valores de R, K, a y g que son constantes. Como se da el valor de t=(15 min)(60 s/min), se unifica la unidad de medida de tiempo a segundos pues g=9.6 m/s².

---

-\frac{4}{3}R h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{Ka}{\pi} t \sqrt{2g}

---

-\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)}

la función para encontrar la raíz se reordena como f(h)=0

f(h) = -\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} +\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)} f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805

Usando la gráfica proporcionada del reservorio semiesférico, el valor de h no puede superar la altura R. Al vaciar el reservorio h=0.

Por lo que el intervalo a usar es [0,3]. Se verifica la existencia de una raíz al con la gráfica y Python.

b. Método de Newton Raphson

El método requiere la derivada f'(h)

f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805 f'(h) = -6\frac{ h^{3/2}}{h} + \frac{h^{5/2}}{h} f'(h) = -6 h^{1/2}+ h^{3/2}

El valor inicial puede ser por ejemplo, el centro del intervalo:

x_0= \frac{a+b}{2} = \frac{0+3}{2} = 1.5

itera=0

f(1.5) = -4(1.5)^{3/2} + 0.4 (1.5)^{5/2} +10.7805= 4.5343 f'(1.5) = -6\frac{ 1.5^{3/2}}{1.5} + \frac{1.5^{5/2}}{1.5} = -5.5113 x_1= x_0 - \frac{f(1.5)}{f'(1.5)}= 1.5-\frac{4.5343}{-5.5113} =2.3227 error_1 =|x1-x0| = |2.3227-1.5| =0.8227

itera=1

f(2.3227) = -4(2.3227)^{3/2} + 0.4 (2.3227)^{5/2} +10.7805= -0.09019 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6043 x_2= 2.3227-\frac{-0.09019}{-5.6043} = 2.3066 error_2 = |2.3066-2.3227| =0.01611

itera=2

f(2.3227) = -4(2.3066)^{3/2} + 0.4 (2.3066)^{5/2} +10.7805= -0.0007104 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6093 x_3= 2.3227-\frac{-0.0007104}{-5.6093} = 2.3066 error_3 = |2.3066-2.3066| =0.000007241

literal c, convergencia

considerando tolerancia de 10-3 (milimétrica), solo serán necesarias 3 iteraciones.

El método converge al mostrarse que los errores disminuyen en cada iteración.

El resultado se encuentra dentro del intervalo y es x = 2.3066

Algoritmo con Python

resultados:

fh:
         ____          ____                   
        /  3          /  5                    
- 4.0*\/  h   + 0.4*\/  h   + 10.7805172327811
dfh:
         ____          ____
        /  3          /  5 
  6.0*\/  h     1.0*\/  h  
- ----------- + -----------
       h             h     
['xi', 'xnuevo', 'tramo']
[[1.5000e+00 2.3227e+00 8.2272e-01]
 [2.3227e+00 2.3066e+00 1.6118e-02]
 [2.3066e+00 2.3066e+00 7.2412e-06]]
raiz en:  2.306612665792577
con error de:  7.241218404896443e-06

>>> f(1.5)
4.534318388683996
>>> df(1.5)
-5.5113519212621505
>>> f(2.3227)
-0.09019959114247378
>>> df(2.3227)
-5.604354799446854
>>> f(2.3066)
7.104683901282272e-05
>>> df(2.3066)
-5.609349350102558

instrucciones:

# 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
h = sym.Symbol('h')
t = 15*60
R = 3 
g = 9.8 
area = 0.01
K = 0.85 
fh = -(4/3)*R*sym.sqrt(h**3)+(2/5)*sym.sqrt(h**5)+(K*area/np.pi)*np.sqrt(2*g)*(t)

# Grafica de f(h)
a = 0
b = 3
muestras = 15

# PROCEDIMIENTO
hi = np.linspace(a,b,muestras)
# para evaluación numérica con numpy
f = sym.lambdify(h,fh)
fi = f(hi)

# derivada con sympy
dfh = sym.diff(fh,h,1)
df = sym.lambdify(h,dfh)

# SALIDA
print('fh:')
sym.pprint(fh)
print('dfh:')
sym.pprint(dfh)

plt.plot(hi,fi)
plt.axhline(0, color='black')
plt.xlabel('h')
plt.ylabel('f')
plt.title(fh)
plt.show()

## Método de Newton-Raphson

# INGRESO
fx  = f
dfx = df
x0 = (a+b)/2 # mitad del intervalo (ejemplo)
tolera = 0.001

# 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)

 

Ejercicio: 3Eva_2023PAOII_T3 Volumen por solido de revolución de un peón

El volumen se calcula a partir de la expresión:

V = \int_{a}^{b} \pi (f(x))^2 dx

xi=[ 0, 3, 5.  , 9.985 , 14.97 , 17.97, 40.04, 43.29, 51.6449, 60]
yi=[15,15,13.25,14.1552, 9.6768,  9.67,  4.64,  4.64,  8.9768, 0.]

Para el intervalo [0,3]

La función es una constante

literal a, b y c

f(x) = 15 V = \int_{0}^{3} \pi (15)^2 dx = \pi (15)^2 x \big|_{0}^{3} = \pi (15)^2 (3-0) = 2120.5750

Para el intervalo [3,5]

literal a

La función es una recta con pendiente negativa

f(x) = -0.875 x + 17.625 V = \int_{3}^{5} \pi (f(x))^2 dx V = \int_{3}^{5} \pi (-0.875 x+17.625)^2 dx

para el integral con cuadratura de Gauss

g(x) = \pi (-0.875*x+17.625)^2

literal b y c

x_a = \frac{5+3}{2} - \frac{5-3}{2}\frac{1}{\sqrt{3}} = 3.4226 x_b = \frac{5+3}{2} + \frac{5-3}{2}\frac{1}{\sqrt{3}} = 4.5773

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{5-3}{2}(g(3.4226) + g(4.5773)) = \frac{5-3}{2} (672.43+582.76) = 1255.1971

Para el intervalo [5, 14.97]

literal a

Del polinomio obtenido en el ejercicio anterior para éste intervalo

f(x) = -0.1083 x^2+1.8047 x + 6.9341 V = \int_{5}^{14.97} \pi (f(x))^2 dx V = \int_{5}^{14.97} \pi (-0.1083 x^2+1.8047 x + 6.9341)^2 dx [ g(x) = \pi (-0.1083 x^2+1.8047 x + 6.9341)^2

literal b y c

x_a = \frac{14.97+5}{2} - \frac{14.97-5}{2}\frac{1}{\sqrt{3}} = 7.1069 x_b = \frac{14.97+5}{2} + \frac{14.97-5}{2}\frac{1}{\sqrt{3}} = 12.8630

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{14.97-5}{2}(g(7.1069) + g(12.8630)) = \frac{14.97-5}{2} (641.5176+469.8124) = 5539.9805

literal d

El volumen total es la suma de los resultados de cada una de las secciones:

V_{total} = 2120.5750 + 1255.1971 + 5539.9805 + ...

Tarea: continuar con los otros intervalos para obtener el volumen total.

literal e

Resultados con el algoritmo

para f(x):
[xa,xb]= [0.6339745962155612, 2.366025403784439]
[f(xa),f(xb)]= [915.4418590078262, 760.106947830205]
Volumenfx:  2513.323210257047

para f(x):
[xa,xb]= [3.4226497308103743, 4.577350269189626]
[f(xa),f(xb)]= [672.4334354361756, 582.7637293668464]
Volumenfx:  1255.1971648030221

para f(x):
[xa,xb]= [7.106908908089714, 12.863091091910285]
[f(xa),f(xb)]= [641.5176069162055, 469.8124936184865]
Volumenfx:  5539.98055116544

Instrucciones en Python

# 3Eva_2023PAOII_T3 Volumen por solido de revolución de un peón
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO
# intervalos eje x, funciones
#a=0 ; b=3 ; f = lambda x: np.pi*(15)**2 + 0*x
#a=3 ; b=5 ; f = lambda x: np.pi*(-0.875*x+17.625)**2
# con el polinomio del tema 2
a=5 ; b=14.97 ; f = lambda x: np.pi*(-0.1083*x**2+1.8047*x + 6.9341)**2
# con el circulo del tema 1
#a=5 ; b=14.97 ; f = lambda x: np.pi*(np.sqrt(6.7**2-(x-8.6)**2) + 7.6)**2

xmuestras = 31

# PROCEDIMIENTO
# Volumen para f(x) con Cuadratura de Gauss
xa = (b+a)/2 + (b-a)/2*(-1/np.sqrt(3)) 
xb = (b+a)/2 + (b-a)/2*(1/np.sqrt(3))
Volumenfx = (b-a)/2*(f(xa)+f(xb))
xab = [xa,xb]
fab = [f(xa),f(xb)]

# SALIDA
print("para f(x):")
print("[xa,xb]=", xab)
print("[f(xa),f(xb)]=", fab)
print("Volumenfx: ",Volumenfx)
print()

Añadir para graficar en 2D

# grafica 2D para el area de corte en eje y,x --------
# muestras en x
xi = np.linspace(a, b, xmuestras)
fi = f(xi)
f0 = np.zeros(xmuestras,dtype=float)

# grafica 2D
figura2D = plt.figure()
grafica2D = figura2D.add_subplot(111)
grafica2D.plot(xi,fi,color='blue',label='f(x)')
grafica2D.fill_between(xi,fi,f0,color='lightblue')
grafica2D.grid()
grafica2D.set_title('Area para sólido de revolución')
grafica2D.set_xlabel('x')
grafica2D.set_ylabel('f(x)')

Añadir para graficar en 3D

# Para grafica 3D -----------------------
# angulo w de rotación
w_a = 0
w_b = 2*np.pi
w_muestras = 31

# grafica 3D muestras en x y angulo w
wi = np.linspace(w_a, w_b, w_muestras)
X, W = np.meshgrid(xi, wi)
# proyeccion en cada eje 
Yf = f(xi)*np.cos(W)
Zf = f(xi)*np.sin(W)

# grafica 3D
figura3D = plt.figure()
grafica = figura3D.add_subplot(111, projection='3d')

grafica.plot_surface(X, Yf, Zf,
color='blue', label='f(x)',
alpha=0.6, rstride=6, cstride=12)

grafica.set_title('Sólido de revolución')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
# grafica.legend()
eleva = 30
rota = -45
deltaw = 5
grafica.view_init(eleva, rota)

# rotacion de ejes
for angulo in range(rota, 360+rota, deltaw ):
grafica.view_init(eleva, angulo)
plt.draw()
plt.pause(.001)
plt.show()

.

Ejercicio: 3Eva_2023PAOII_T2 perfil de un peón

literal a y b

Para el planteamiento de los polinomios, es necesario observar la gráfica proporcionada para el ejercicio y la correspondiente tabla de datos:

xi	0	3	5	9.985	14.97	17.97	40.04	43.29	51.6456	60
yi	15	15	13.25	14.155	9.676	9.676	4.64	4.64	8.976	0

Para el intervalo [0,3], El perfil se describe con una constante, por lo tanto:

p_1(x) = 15

Para el intervalo [3,5] el perfil se describe como una recta con pendiente negativa.

p_2(x) =a_0 + a_1 x a_1 = \frac{ 13.25-15}{5-3} = -0.875 p_2(5) = 13.25 =a_0 -0.875 (5) a_0 = 13.25 +0.875 (5) = 17.625 p_2(x) = -0.875 x + 17.625

Para el con los puntos xi = [5, 9.985, 14.97] el perfil se describe con tres puntos, por lo que el polinomio de interpolación puede ser de grado 2. Usando el método de Lagrange:

p_3(x) =\frac{(x-9.985)( x-14.97)}{(5-9.985)( 5-14.97)} 13.25 + +\frac{(x-5)( x-14.97)}{(9.985-5)( 9.985-14.97)} 14.155 + +\frac{(x-5)( x-9.985)}{(14.97-5)( 14.97-9.985)} 9.676 +

simplificando la expresión con Python y Sympy:

P_3(x) = -0.1083 x^2 + 1.8047 x + 6.9341

literal c

El error para los polinomios acorde a la gráfica proporcionada es cero para los tramos [0,3] y [3,5], dado que el primero es una constante y el segundo es una recta con pendiente.

El error de mayor magnitud se presenta con la interpolación entre el círculo de la ecuación del tema 1 y el polinomio de grado 2 en el intervalo [5, 14.97]. Tomando como ejemplo el punto intermedio del tramo derecho en:

x_k = \frac{14.97-5}{4} + 9.985= 12.4775 p_3(12.4775) = -0.1083(12.4775)^2 + 1.8047(12.4775) + 6.9341 = 12.5889 f(12.4775) = \sqrt{6.7^2-(12.4775-8.6)^2} + 7.6 = 13.0639 errado = |p_3(12.4775) - f(12.4775)| = 0.4750

literal d

Se puede mejorar la aproximación del polinomio aumentando el grado y número de puntos a usar dentro del intervalo. El resultado se debería acercar mas al perfil superior del círculo de referencia del tema 1.

Resultados para el intervalo [5, 14.97]

    valores de fi:  [13.25   14.1552  9.6768]
divisores en L(i):  [ 49.70045  -24.850225  49.70045 ]

Polinomio de Lagrange, expresiones
0.266597183727713*(x - 14.97)*(x - 9.985) 
- 0.569620596996607*(x - 14.97)*(x - 5.0) 
+ 0.194702462452553*(x - 9.985)*(x - 5.0)

Polinomio de Lagrange: 
-0.108320950816341*x**2 + 1.80477420224566*x + 6.93415275918024

Instrucciones Python

# 3Eva_2023PAOII_T2 perfil de un peón
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
xj=[ 0, 3, 5.  , 9.985 , 14.97 , 17.97, 40.04, 43.29, 51.6449, 60]
yj=[15,15,13.25,14.1552, 9.6768,  9.67,  4.64,  4.64,  8.9768, 0.]


#  Datos de prueba
ia = 2
ib = 4+1
xi = np.array(xj[ia:ib])
fi = np.array(yj[ia:ib])

g = lambda x: np.sqrt(6.7**2-(x-8.6)**2) + 7.6
muestras = 21

# PROCEDIMIENTO
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación 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)
gi = g(pxi)

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.plot(pxi,gi, label = 'g(x)')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.grid()
plt.show()