Ejercicio: 1Eva2017TI_T2 Tanque esférico-volumen
a. Planteamiento del problema
V = \frac{\pi h^{2} (3r-h)}{3}Si r=1 m y V=0.75 m3,
0.75 = \frac{\pi h^{2} (3(1)-h)}{3} 0.75 -\frac{\pi h^{2} (3(1)-h)}{3} = 0 f(h) = 0.75 -\frac{\pi h^{2} (3-h)}{3} = 0.75 -\frac{\pi}{3}(h^{2} (3)-h^{2}h) = 0.75 -\frac{\pi}{3}(3 h^{2} - h^{3}) f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}
b. Intervalo de búsqueda de raíz
El tanque vació tiene h=0 y completamente lleno h= 2r = 2(1) = 2, por lo que el intervalo tiene como extremos:
[0,2]
Verificando que exista cambio de signo en la función:
f(0) = 0.75 -\pi (0)^{2} + \frac{\pi}{3}(0)^{3} = 0.75 f(2) = 0.75 -\pi (2)^{2} + \frac{\pi}{3}(2)^{3}= -3.4387
y verificando al usar la gráfica dentro del intervalo:
Tolerancia
Se indica en el enunciado como 0.01 que es la medición mínima a observar con un flexómetro.
tolera = 0.01
c. Método de Newton-Raphson
d. Método de Punto Fijo
c. Método de Newton-Raphson
El método de Newton-Raphson requiere la derivada de la función:
x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3} f'(h) = -2\pi h + \pi h^{2}Tomando como punto inicial de búsqueda el extremo izquierdo del intervalo, genera una división para cero. Por lo que se mueve un poco a la derecha, algo más cercano a la raiz, viendo la gráfica por ejemplo 0.1
x0 = 0.1
iteración 1
i =0
f(0.1) = 0.75 -\pi (0.1)^{2} + \frac{\pi}{3}(0.1)^{3} =0.7196 f'(0.1) = -2\pi (0.1) + \pi (0.1)^{2} = -0.5969 x_{1} = x_0 -\frac{0.7496 }{-0.0625} = 1.3056 tramo = |x_{1}-x_{0}| = |0.1-1.3056 | = 1.2056iteración 2
i =1
f(1.3056) = 0.75 -\pi (1.3056)^{2} + \frac{\pi}{3}(1.3056)^{3} = -2.2746 f'(1.3056) = -2\pi (1.3056) + \pi (1.3056)^{2} =-2.8481 x_{1} = x_0 -\frac{-2.2746}{-2.8481} = 0.5069 tramo = |x_{2}-x_{1}|=|0.5069-1.3056|=0.7987iteración 3
i =2
f(0.5069) = 0.75 -\pi (0.5069)^{2} + \frac{\pi}{3}(0.5069)^{3} = 0.0789 f'(0.5069) = -2\pi (0.5069) + \pi (0.5069)^{2} =-2.3780 x_{1} = x_0 -\frac{0.0789}{-2.3780} = 0.5401 tramo = |x_{3}-x_{2}| =|0.5401 - 0.5069| = 0.0332Observe que el tramo disminuye en cada iteración , por lo que el método converge, si se siguen haciendo las operaciones se tiene que:
método de Newton-Raphson
i ['xi', 'fi', 'dfi', 'xnuevo', 'tramo']
0 [ 0.1 0.7196 -0.5969 1.3056 1.2056]
1 [ 1.3056 -2.2746 -2.8482 0.507 0.7986]
2 [ 0.507 0.079 -2.378 0.5402 0.0332]
3 [ 5.4019e-01 -1.6689e-03 -2.4774e+00 5.3952e-01 6.7367e-04]
raíz en: 0.5395186666699257Algoritmo en Python
para Método de Newton-Raphson
# 1Eva2017TI_T2 Tanque esférico-volumen
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)
import numpy as np
def newton_raphson(fx,dfx,x0, tolera, iteramax=100,
vertabla=False, precision=4):
'''fx y dfx en forma numérica lambda
xi es el punto inicial de búsqueda
si no converge hasta iteramax iteraciones
la respuesta es NaN (Not a Number)
'''
itera=0
xi = x0
tramo = abs(2*tolera)
if vertabla==True:
print('método de Newton-Raphson')
print('i', ['xi','fi','dfi', 'xnuevo', 'tramo'])
np.set_printoptions(precision)
while (tramo>=tolera):
fi = fx(xi)
dfi = dfx(xi)
xnuevo = xi - fi/dfi
tramo = abs(xnuevo-xi)
if vertabla==True:
print(itera,np.array([xi,fi,dfi,xnuevo,tramo]))
xi = xnuevo
itera = itera + 1
if itera>=iteramax:
xi = np.nan
print('itera: ',itera,
'No converge,se alcanzó el máximo de iteraciones')
return(xi)
# INGRESO
fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
dfx = lambda h: -2*np.pi*h+np.pi*(h**2)
# Parámetros de método
x0 = 0.1
tolera = 0.01
iteramax = 100
# PROCEDIMIENTO
xi= newton_raphson(fx,dfx,x0, tolera, vertabla=True)
# SALIDA
print('raíz en: ', xi)
# GRAFICA
import matplotlib.pyplot as plt
a = 0
b = 2
muestras = 21
xj = np.linspace(a,b,muestras)
fj = fx(xj)
plt.plot(xj,fj, label='f(x)')
plt.plot(xi,0, 'o')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
d. Método de Punto Fijo
Del planteamiento del problema en el literal a, se tiene que:
0.75 = \frac{\pi h^{2} (3(1)-h)}{3}de donde se despeja una h:
\frac{3(0.75)}{\pi (3(1)-h) } = h^{2} h = \sqrt{\frac{3*0.75}{\pi (3-h) }} h = \sqrt{\frac{2.25}{\pi (3-h) }}con lo que se obtienen las expresiones a usar en el método
identidad = h g(h) = \sqrt{\frac{2.25}{\pi (3-h) }}El punto inicial de búsqueda debe encontrarse en el intervalo, se toma el mismo valor que x0 en el método de Newton-Raphson
x0 = 0.10
Iteración 1
x_0= 0.10 g(0.10) = \sqrt{\frac{2.25}{\pi (3-(0.10) }}= 0.4969 tramo= |0.4969-0.10| = 0.3869Iteración 2
x_1= 0.4969 g(0.4969) = \sqrt{\frac{2.25}{\pi (3-(0.4969 ) }}= 0.5349 tramo= |0.5349- 0.4969| = 0.038Iteración 3
x_2 =0.5349 g(0.5349) = \sqrt{\frac{2.25}{\pi (3-(0.5349) }}= 0.5390 tramo= |0.5390 - 0.5349| = 0.0041con lo que se cumple el criterio de tolerancia, y se obtiene la raíz de:
raiz = 0.5394Tabla de resultados, donde se observa que el tramo o error en cada iteración disminuye, por lo que el método converge.
método del Punto Fijo
i ['xi', 'gi', 'tramo']
0 [0.1 0.496955 0.396955]
1 [0.496955 0.534912 0.037956]
2 [0.534912 0.539014 0.004102]
3 [5.390140e-01 5.394631e-01 4.490776e-04]
raíz en: 0.539463113150304
>>> Algoritmo en Python
para Método de Punto-Fijo, recordamos que el método puede ser divergente, por lo que se añade el parámetro iteramax
# 1Eva2017TI_T2 Tanque esférico-volumen
# Algoritmo de punto fijo
import numpy as np
def puntofijo(gx,c,tolera,iteramax=40,vertabla=True, precision=6):
"""
g(x) se obtiene al despejar una x de f(x)
máximo de iteraciones predeterminado: iteramax
si no converge hasta iteramax iteraciones
la respuesta es NaN (Not a Number)
"""
itera = 0 # iteración inicial
tramo = 2*tolera # al menos una iteracion
if vertabla==True:
print('método del Punto Fijo')
print('i', ['xi','gi','tramo'])
np.set_printoptions(precision)
while (tramo>=tolera and itera<=iteramax):
gc = gx(c)
tramo = abs(gc-c)
if vertabla==True:
print(itera,np.array([c,gc,tramo]))
c = gc
itera = itera + 1
respuesta = c
# Valida respuesta
if itera>=iteramax:
respuesta = np.nan
print('itera: ',itera,
'No converge,se alcanzó el máximo de iteraciones')
return(respuesta)
# PROGRAMA ----------------------
# INGRESO
fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
gx = lambda h: np.sqrt(2.25/(np.pi*(3-h)))
c = 0.10 # valor inicial
tolera = 0.001
iteramax = 15
# PROCEDIMIENTO
c = puntofijo(gx,c,tolera,iteramax,vertabla=True)
# SALIDA
print('raíz en: ', c)
# GRAFICA
import matplotlib.pyplot as plt
a = 0.25 # intervalo de gráfica en [a,b]
b = 1
muestras = 21
# calcula los puntos para fx y gx
xj = np.linspace(a,b,muestras)
fj = fx(xj)
gj = gx(xj)
plt.plot(xj,fj, label='f(x)',
linestyle='dashed')
plt.plot(xj,gj, label='g(x)')
plt.plot(xj,xj, label='y=x')
plt.plot(c,c, 'o')
plt.axvline(c, color='magenta',
linestyle='dotted')
plt.axhline(0)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Punto Fijo')
plt.grid()
plt.legend()
plt.show()
Otra forma de probar la convergencia es que |g'(x)|<1 que se observa en la una gráfica adicional, lo que limita aún más el intervalo de búsqueda.
Desarrollo en la siguiente clase.