Método Punto Fijo: Calcular profundidad de un tanque

PROBLEMA: El volumen V de un líquido contenido en un tanque horizontal cilíndrico de radio r y longitud L está relacionado con la profundidad del líquido h por
$$V=\left[ { r }^{ 2 }{ cos }^{ -1 }\left( \frac { r-h }{ r }  \right) -\left( r-h\sqrt { 2rh-{ h }^{ 2 } }  \right)  \right] L$$
Determine h para r = 2 m, L = 5 m y V = 8.5 m3.

(Ejercicio 8.8, tomado del libro Métodos Numéricos para Ingenieros de Steven C. Chapra y Raymond P. Canale. 5Ed. McGraw Hill. Mexico. 2007. Página 218)


SOLUCIÓN

Análisis de la situación

En primer lugar, es necesario que sustituyamos los datos que nos dan en el enunciado en la ecuación
$$V=\left[ { r }^{ 2 }{ cos }^{ -1 }\left( \frac { r-h }{ r }  \right) -\left( r-h\sqrt { 2rh-{ h }^{ 2 } }  \right)  \right] L,$$

es decir,
$$V=\left[ { 2 }^{ 2 }{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -\left( 2-h\sqrt { 2*2*h-{ h }^{ 2 } }  \right)  \right] * 5$$
$$8.5=20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -\left( 10-5h\sqrt { 4h-{ h }^{ 2 } }  \right) $$

Como segundo, tenemos que determinar los valores de h para los que se cumple la igualdad.
Una forma alterna de lograr lo mismo es construir una función en una variable (h) a partir de la ecuación

$$8.5=20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -\left( 10-5h\sqrt { 4h-{ h }^{ 2 } }  \right) $$

y determinar todos los interceptos con el eje horizontal. Estos valores serían los valores de h para los cuales también se satisface la ecuación original.

Construyamos la función a partir de la ecuación simplificada:
$$8.5=20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -\left( 10-5h\sqrt { 4h-{ h }^{ 2 } }  \right) $$
$$8.5-20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) +\left( 10-5h\sqrt { 4h-{ h }^{ 2 } }  \right)=0 $$
$$18.5-20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -5h\sqrt { 4h-{ h }^{ 2 } }  =0 $$

Por último, tomemos el miembro izquierdo de la igualdad como la función, es decir,
$$f(h)=18.5-20{ cos }^{ -1 }\left( \frac { 2-h }{ 2 }  \right) -5h\sqrt { 4h-{ h }^{ 2 } }  $$

Listo ya tenemos la función.

Ahora, observemos la gráfica de f(h) para tener idea de cuantos interceptos hay con el eje horizontal.
Por comodidad, cambiaré la variable ‘h’ por ‘x’, es decir, la función quedaría
$$f(x)=18.5-20{ cos }^{ -1 }\left( \frac { 2-x }{ 2 }  \right) -5x\sqrt { 4x-{ x }^{ 2 } }  $$

Ahora sí, veamos su gráfica:

Gráfica 1

Vemos en la Gráfica 1, que hay un único intercepto, y lo podemos ubicar, por ejemplo, en el intervalo [0,10]. En teoría este sería el valor de h que necesitamos encontrar.

Aplicación del método de punto fijo

Antes de usar el código que implementa el método de punto fijo, es necesario que entendamos como se debe proceder a la hora aplicar correctamente el método.

Pasos a seguir:
1. Verificar que la función tenga intercepto con el eje horizontal. En nuestro caso la función f(x) cumple con esto.
2. A partir de f(x)=0, es necesario despejar una x de modo que  nos quede en la forma x=g(x). Veamos la idea con nuestra función:
$$f(x)=18.5-20{ cos }^{ -1 }\left( \frac { 2-x }{ 2 }  \right) -5x\sqrt { 4x-{ x }^{ 2 } }=0  $$

Quitar el término f(x) y despejar a x:
$$18.5-20{ cos }^{ -1 }\left( \frac { 2-x }{ 2 }  \right) -5x\sqrt { 4x-{ x }^{ 2 } }=0  $$
$$-20{ cos }^{ -1 }\left( \frac { 2-x }{ 2 }  \right) = 5x\sqrt { 4x-{ x }^{ 2 } } - 18.5  $$
$${ cos }^{ -1 }\left( \frac { 2-x }{ 2 }  \right) =\frac { 5x\sqrt { 4x-{ x }^{ 2 } } -18.5 }{ -20 } $$
$$\frac { 2-x }{ 2 } =cos\left( \frac { 5x\sqrt { 4x-{ x }^{ 2 } } -18.5 }{ -20 }  \right) $$
$$-x\quad =\quad 2cos\left( \frac { 5x\sqrt { 4x-{ x }^{ 2 } } -18.5 }{ -20 }  \right) -2$$
$$x\quad =\quad 2-2cos\left( \frac { 5x\sqrt { 4x-{ x }^{ 2 } } -18.5 }{ -20 }  \right) $$
Por lo tanto,
$$g(x)=\quad =\quad 2-2cos\left( \frac { 5x\sqrt { 4x-{ x }^{ 2 } } -18.5 }{ -20 }  \right) $$

3. g(x) debe satisfacer las condiciones del Teorema del Punto Fijo:

Teorema del Punto Fijo

Si g(x) es una función continua en [a,b] y los valores resultantes de g(x) están en [a,b] para todo x en [a,b], entonces g(x) tiene por lo menos un punto fijo en [a,b]. Si además, g(x) existe para todo x en (a,b) y |g(x)| ≤ K < 1 para todo x en (a,b), siendo K constante, entonces g(x) tiene un único punto fijo a en [a,b] y la sucesión {Xn}n definida mediante la fórmula de iteración
Xn = g(Xn-1), n=1,2,3,…
converge a a, cualquiera que se X0 que este en [a,b].

En otras palabras, para poder aplicar el método de punto fijo sobre g(x), hay que verificar dos aspectos del teorema:

a) Que g(x) sea continua en un intervalo [a,b].
b) Que |g(x)| < 1 para todo x en (a,b).

Apoyémonos en la gráfica de g(x) para verificar estas dos condiciones:

Gráfica 2

Al observar la Gráfica 2, notamos que no hay interrupción o cortes en la gráfica de g(x) sobre el intervalo de [0,4]. Con la mirada puesta en la solución general de nuestro problema, tomemos un intervalo para g(x) que contenga el intercepto de f(x), y en el analicemos las condiciones del teorema de punto fijo.

Hagamos un zoom a la gráfica de f(x) en el intervalo [0,4] para elegir un intervalo para g(x) que nos favorezca el análisis:

Gráfica 3

Observando la Gráfica 3, vemos que el intercepto de f(x) lo podemos ubicar en el intervalo [0,1]. Tomemos este intervalo para analizar las condiciones del teorema del punto fijo.


Análisis de las condiciones del teorema del punto fijo

a) Verificar que g(x) es continua en [0,1].
Comencemos con la función coseno. Sabemos que la función coseno es continua en todos los números reales, por lo tanto si hubiese un problema de continuidad, éste estaría en
$$\sqrt { 4x-{ x }^{ 2 } } $$
Pero sabemos también que esta raíz está definida completamente en el intervalo de [0,4] y que allí no hay discontinuidad.
En conclusión g(x) es continua en [0,1], puesto que [0,1] está contenido en [0,4].

b) |g’(x)| < 1 para todo x en (0,1).
Primero calculemos g’(x). Usando los comandos syms, diff y simplify de Matlab podemos calcular la derivada de g(x):
>> syms x
>> simplify(diff(2-2*cos((5*x*sqrt(4*x-x^2)-18.5)/(-20))))
ans =
-(x*sin((x*(-x*(x - 4))^(1/2))/4 - 37/40)*(x - 3))/(-x*(x - 4))^(1/2)

Luego |g(x)| = abs(-(x*sin((x*(-x*(x - 4))^(1/2))/4 - 37/40)*(x - 3))/(-x*(x - 4))^(1/2))

Ahora grafiquemos a |g(x)|:

Gráfica 4

En la Gráfica 4, podemos observar que ninguno de los valores de |g(x)| supera a 1, ni siquiera a 0.8. Por lo tanto, podemos concluir que |g(x)|<1 para todo x en (0,1).

Bien, hemos visto que la función g(x) satisface las condiciones del teorema de punto fijo y con esto podemos concluir que g(x) tiene un único punto fijo.
Ahora sí, determinemos este punto fijo usando el código en Matlab que implementa el Método de Punto Fijo.

Aplicación del Método del Punto Fijo

Para aplicar este método necesitamos la función g(x) y un punto inicial. Según el teorema del punto fijo éste punto inicial puede ser cualquier valor que este dentro del intervalo [0,1]. Tomemos a 0.1 como punto inicial.

El código en Matlab que implementa el método de punto fijo es:

fun=input('Ingrese la funcion g(x): ','s');
f=inline(fun);
xi=input('Ingrese valor inicial: ');
e=input('Ingrese porcentaje de error: ');

xa=f(xi);
ea=abs((xa-xi)/xa)*100;

while(ea>e)
    xi=xa;
    xa=f(xi);
    ea=abs((xa-xi)/xa)*100;
end
fprintf('La raiz aproximada es:  %12.10f \n',xa);

el cual he llamado puntoFijo.m.

Ejecutemos e ingresemos los datos requeridos:

>> puntoFijo
Ingrese la funcion g(x): 2-2*cos((5*x*sqrt(4*x-x^2)-18.5)/(-20))
Ingrese valor inicial: 0.1
Ingrese porcentaje de error: 10^(-10)
La raiz aproximada es:  0.5301496449

Por lo tanto, el valor de la profundidad que queríamos determinar es h=0.5301496449 m.

Verifiquemos que efectivamente h=0.5301496449 es el valor que estábamos buscando y lo haremos evaluándolo en f(x):
$$f(0.5301496449)=18.5-20{ cos }^{ -1 }\left( \frac { 2-0.5301496449 }{ 2 }  \right) \\-5(0.5301496449)\sqrt { 4(0.5301496449)-{ 0.5301496449 }^{ 2 } } $$
$$f(x)=0$$
Perfecto, esto nos permite concluir que la profundidad que satisface los datos del enunciado del problema es h=0.5301496449 m.

Para aquellos que son fanáticos de Python, les dejo el código del Método de Punto Fijo a continuación:

gx = input('Ingrese la funcion g(x): ')
x0 = float(input('Ingrese el valor inicial: '))
e = float(input('Ingrese porcentaje de error: '))

g = sympify(gx)
x1 = g.subs('x',x0)
em = abs((x1-x0)/x1)*100

while e<em:
    if(g.subs('x',x1)==0):
        em = 0
    else:
        x0 = x1
        x1 = g.subs('x',x0)
        em = abs((x1-x0)/x1)*100

print('La raiz aproximada es: ',x1)

Bueno, hemos llegado al final de esta entrada.

Espero les sirva de mucha ayuda.

Comentarios

Publicar un comentario