Autor: Edison Del Rosario

Ejercicio: 2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

\frac{\delta}{\delta t}x(t) = b x(t) - d (x(t))^2 \frac{\delta}{\delta t}y(t) = b y(t) - d (y(t))^2 +r b (x(t)-y(t))

literal a

simplificando la nomenclatura

x' = b x - d x^2 y' = b y - d y^2 +r b (x-y)

sustituyendo constantes, y considerando x(0)=1 ; y(0)=0.01 ; h=0.5

x' = 0.02 x - 0.015 x^2 y' = 0.02 y - 0.015 y^2 +0.1(0.02) (x-y)

el planteamiento de Runge Kutta se hace junto a la primera iteración, además de encontrarse en las instrucciones con Python.

literal b

Se describen 3 iteraciones usando los resultados de la tabla con el algoritmo, para mostrar la comprensión del algoritmo.

t = 0
K1x = 0.5 (0.02 (1) – 0.015 (1)² = 0.0025
K1y = 0.5(0.02 (0.01) – 0.015 (0.01)² +0.1(0.02) (1-0.01)= 0.001089

K2x = 0.5 (0.02 (1+0.0025) – 0.015 (1+0.0025)²= 0.00248
K2y = 0.5(0.02 (0.01+0.00108) – 0.015 (0.01+0.00108)²+0.1(0.02) ((1+0.0025)-(0.01+0.00108)) = 0.001101

x₁ = 1 + (1/2)(0.0025+0.00248) = 1.0025
y₁ = 0.01 + (1/2)(0.001089+0.001101) = 0.01109
t₁ = 0 + 0.5 =0.5

t=0.5
K1x = 0.5 (0.02 (1.0025) – 0.015 (1.0025)² = 0.002487
K1y = 0.5(0.02 (0.01109) – 0.015 (0.01109)² +0.1(0.02) (1.0025-0.01109)= 0.001101

K2x = 0.5 (0.02 (1.0025+ 0.002487) – 0.015 (1.0025+ 0.002487)²= 0.002474
K2y = 0.5(0.02 (0.01109+0.001101) – 0.015 (0.01109+0.001101)²+0.1(0.02) ((1.0025+ 0.002487)-(0.01109+0.001101)) = 0.001113

x₂ = 1.0025 + (1/2)(0.002487+0.002474) = 1.0050
y₂ = 0.01109 + (1/2)(0.001101+0.001113) = 0.01220
t₂ = 0.5 + 0.5 = 1

t=1
K1x = 0.5 (0.02 (1.0050) – 0.015 (1.0050)² = 0.002474
K1y = 0.5(0.02 (0.01220) – 0.015 (0.01220)² +0.1(0.02) (1.0050-0.01220)= 0.001113

K2x = 0.5 (0.02 (1.0050+ 0.002474) – 0.015 (1.0050+ 0.002474)²= 0.002462
K2y = 0.5(0.02 (0.01220+0.001113) – 0.015 (0.01220+0.001113)²+0.1(0.02) ((1.0050+ 0.002474)-(0.01220+0.001113)) = 0.001126

x₃ = 1.0050 + (1/2)(0.002474+0.002462) = 1.0074
y₃ = 0.01220 + (1/2)(0.001113+0.001126) = 0.01332
t₃ = 1 + 0.5 = 1.5

Resultado con el algoritmo

Para obtener los datos de las iteraciones, primero se ejecuta el algoritmo para pocas iteraciones.
Para la pregunta sobre 200 años, se incrementa las iteraciones a 2 por año y las condiciones iniciales, es decir 401 muestras.

 [ ti, xi, yi]
 [ ti, xi, yi]
[[0.0000e+00 1.0000e+00 1.0000e-02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
 [5.0000e-01 1.0025e+00 1.1095e-02 2.5000e-03 1.0892e-03 2.4875e-03 1.1014e-03]
 [1.0000e+00 1.0050e+00 1.2203e-02 2.4875e-03 1.1014e-03 2.4749e-03 1.1136e-03]
 [1.5000e+00 1.0074e+00 1.3323e-02 2.4749e-03 1.1137e-03 2.4623e-03 1.1260e-03]
 [2.0000e+00 1.0099e+00 1.4455e-02 2.4624e-03 1.1260e-03 2.4497e-03 1.1384e-03]
 [2.5000e+00 1.0123e+00 1.5600e-02 2.4498e-03 1.1384e-03 2.4371e-03 1.1509e-03]
 [3.0000e+00 1.0148e+00 1.6757e-02 2.4371e-03 1.1509e-03 2.4245e-03 1.1634e-03]
 [3.5000e+00 1.0172e+00 1.7926e-02 2.4245e-03 1.1635e-03 2.4118e-03 1.1761e-03]
 [4.0000e+00 1.0196e+00 1.9109e-02 2.4118e-03 1.1761e-03 2.3991e-03 1.1888e-03]
...
 [1.9950e+02 1.3252e+00 1.1561e+00 ... 1.7202e-03 8.1217e-05 1.7059e-03]
 [2.0000e+02 1.3252e+00 1.1578e+00 ... 1.7060e-03 8.0418e-05 1.6918e-03]
 [2.0050e+02 1.3253e+00 1.1595e+00 ... 1.6919e-03 7.9628e-05 1.6778e-03]]
>>> 

Observación: La población identificada como protestante, continua creciendo, mientras que la proporción de «conformistas» se reduce según los parámetros indicados en el ejercicio. Los valores de natalidad y defunción cambian con el tiempo mucho más en años por otras variables, por lo que se deben realizar ajustes si se pretende extender el modelo.

Instrucciones en Python

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,7),dtype=float)
    tabla[0] = [t0,x0,y0,0,0,0,0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,K1x,K1y,K2x,K2y]
    tabla = np.array(tabla)
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
b = 0.02
d = 0.015
r = 0.1

# Ecuaciones
f = lambda t,x,y : (b-d*x)*x
g = lambda t,x,y : (b-d*y)*y + r*b*(x-y)

# Condiciones iniciales
t0 = 0
x0 = 1
y0 = 0.01

# parámetros del algoritmo
h = 0.5
muestras = 401

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi, K1x, K1y, K2x, K2y]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt

plt.plot(ti,xi, label='xi poblacion')
plt.plot(ti,yi, label='yi protestante')

plt.title('población y protestantes')
plt.xlabel('t años')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('población y protestantes [xi,yi]')
plt.xlabel('x población')
plt.ylabel('y protestantes')
plt.grid()
plt.show()

Ejercicio: 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

literal a

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt v = 1800 \ln\Big(\frac{160000}{160000-2500t}\Big) - 9.8t

Seleccionando el método de Simpson de 3/8, se requieren al menos 3 tramos o segmentos para usarlo, que generan 4 muestras. El vector de tiempo se obtiene como:

v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3
h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)

siendo los vectores:

ti = [ 0. 10. 20. 30.]
vi = [ 0. 207.81826623 478.44820899 844.54060574]

la aplicación del método de Simpson de 3/8 es:

I = \frac{3}{8}(10) \Bigg(1800 \ln\Big(\frac{160000}{160000-2500(0)}\Big) - 9.8(0) +3(1800 \ln\Big(\frac{160000}{160000-2500(10)}\Big) - 9.8(10)) +3(1800 \ln\Big(\frac{160000}{160000-2500(20)}\Big) - 9.8(20)) +1800 \ln\Big(\frac{160000}{160000-2500(30)}\Big) - 9.8(30) \Bigg) = I = \frac{3}{8}(10) \Big(v(0)+3(v(10))+3(v(20))+v(30) \Big) I = \frac{3}{8}(10) \Big(0+3(207.81)+3(478.44)+844.54 \Big) I = 10887.52

literal b

para el primer segmento se usa t entre [0,10]

x_a = \frac{0+10}{2} + \frac{1}{\sqrt{3}}\frac{10-0}{2} = 7.88 x_b = \frac{0+10}{2} - \frac{1}{\sqrt{3}}\frac{10-0}{2} = 2.11 I = \frac{10-0}{2}\Big(v(7.88)+v(2.11)\Big)=995.79

para el 2do segmento se usa t entre [10,20]

x_a = \frac{10+20}{2} + \frac{1}{\sqrt{3}}\frac{20-10}{2} = 17.88 x_b = \frac{10+20}{2} - \frac{1}{\sqrt{3}}\frac{20-10}{2} = 12.11 I = \frac{20-10}{2}\Big(v(17.88)+v(12.11)\Big) =3368.42

para el 3er segmento se usa t entre [20,30]

x_a = \frac{20+30}{2} + \frac{1}{\sqrt{3}}\frac{30-20}{2} = 27.88 x_b = \frac{20+30}{2} - \frac{1}{\sqrt{3}}\frac{30-20}{2} = 22.11 I = \frac{30-20}{2}\Big(v(27.88)+v(22.11)\Big) = 6515.23 Altura = 995.79+ 3368.42 + 6515.23 = 10879.44

literal c

el error es la diferencia entre los métodos
error_entre = |10887.52-10879.44| = 8.079

Resultados con algoritmo

Método de Simpon 3/8
ti
[ 0. 10. 20. 30.]
vi
[ 0. 207.81826623 478.44820899 844.54060574]
Altura con Simpson 3/8 : 10887.52511781406
segmento Cuad_Gauss :    [995.792, 3368.421, 6515.231]
Altura Cuadratura Gauss: 10879.445437288954
diferencia s3/8 y Cuad_Gauss: 8.079680525106596
>>>

Instrucciones en Python

# 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos
import numpy as np

# INGRESO
v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3

# PROCEDIMIENTO literal a
def integrasimpson38_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 3/8
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-3)):
        h  = xi[i+1]-xi[i]
        h1 = (xi[i+2]-xi[i+1])
        h2 = (xi[i+3]-xi[i+2])
        dh = abs(h-h1)+abs(h-h2)
        if dh<tolera:# tramos iguales
            unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
            unS38 = (3/8)*h*unS38
            suma = suma + unS38
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 3
    if (i+1)<n: # incompleto, tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma +str(n-(i+1))+' tramos'
    return(suma)

h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)
altura = integrasimpson38_fi(ti,vi)

# SALIDA
print('Método de Simpon 3/8')
print('ti')
print(ti)
print('vi')
print(vi)
print('Altura con Simpson 3/8 :',altura)

# PROCEDIMIENTO literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

area = 0
area_i =[]
for i in range(0,tramos,1):
    deltaA = integraCuadGauss2p(v,ti[i],ti[i+1])
    area = area + deltaA
    area_i.append(deltaA)
# SALIDA
print('segmento Cuad_Gauss :   ', area_i)
print('Altura Cuadratura Gauss:', area)

print('diferencia s3/8 y Cuad_Gauss:',altura-area)

import matplotlib.pyplot as plt
plt.plot(ti,vi)
plt.plot(ti,vi,'o')
plt.title('v(t)')
plt.xlabel('t (s)')
plt.ylabel('v (m/s)')
plt.grid()
plt.show()

Ejercicio: 1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

La variable independiente para la trayectoria es tiempo, con datos en el vector de ti.

ti  = [0, 1, 2, 3, 4]
xti = [2, 1, 3, 4, 2]
yti = [0, 1, 5, 1, 0]

Considerando que los puntos marcan posiciones por donde debe pasar el dron y se define la trayectoria, se usarán todos los puntos. Cada polinomio será de grado 4 al incluir los 5 puntos disponibles para cada eje.

Nota: podría usar polinomios de menor grado, siempre que considere que se debe completar la trayectoria y regresar al punto de salida.

px(t) = 2\frac{(t-1)(t-2)(t-3)(t-4)}{(0-1)(0-2)(0-3)(0-4)} + 1 \frac{(t-0)(t-2)(t-3)(t-4)}{(1-0)(1-2)(1-3)(1-4)} + 3 \frac{(t-0)(t-1)(t-3)(t-4)}{(2-0)(2-1)(2-3)(2-4)} + 4 \frac{(t-0)(t-1)(t-2)(t-4)}{(3-0)(3-1)(3-2)(3-4)} + 2 \frac{(t-0)(t-1)(t-2)(t-3)}{(4-0)(4-1)(4-2)(4-3)}

simplificando con el algoritmo:

px(t) = \frac{1}{12}t^4 - \frac{7}{6}t^3 + \frac{53}{12}t^2 - \frac{13}{3}t + 2

Realizando lo mismo con el algoritmo para polinomio de Lagrange se obtiene:

py(t) = \frac{11}{12}t^4 - \frac{22}{3}t^3 + \frac{205}{12}t^2 - \frac{29}{3}t

se muestra la gráfica de trayectorias por cada eje vs tiempo

Observaciones: La trayectoria usada tiene el mismo punto de salida como de retorno. La trayectoria presenta lóbulos que podrían ser reducidos y minimizar uso de recursos como bateria. Considere usar trazadores cúbicos y observe la misma gráfica de trayectorias x(t) vs y(t).

Resultado con el algoritmo:

Polinomio de Lagrange x: 
x**4/12 - 7*x**3/6 + 53*x**2/12 - 13*x/3 + 2
Polinomio de Lagrange y: 
11*x**4/12 - 22*x**3/3 + 205*x**2/12 - 29*x/3

Algoritmo en Python

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
ti  = [0,1,2,3,4]
xti = [2,1,3,4,2]
yti = [0,1,5,1,0]

# PROCEDIMIENTO
x = sym.Symbol('x')

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.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*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# para ejex
polinomiox = interpola_lagrange(ti,xti)
polisimplex = polinomiox.expand()
px = sym.lambdify(x,polisimplex)

# para ejey
polinomioy = interpola_lagrange(ti,yti)
polisimpley = polinomioy.expand()
py = sym.lambdify(x,polisimpley)

# Puntos para la gráfica
muestras = 101
a = np.min(ti)
b = np.max(ti)
ti = np.linspace(a,b,muestras)
pxi = px(ti)
pyi = py(ti)

# SALIDA
print('Polinomio de Lagrange x: ')
print(polisimplex)
print('Polinomio de Lagrange y: ')
print(polisimpley)

# Gráfica
figura, enplano = plt.subplots()
plt.scatter(xti,yti, color='red')
plt.plot(pxi,pyi)
plt.ylabel('y(t)')
plt.xlabel('x(t)')
plt.title('trayectoria 2D')
plt.grid()

figura, entiempo = plt.subplots()
plt.plot(ti,pxi, label = 'px')
plt.plot(ti,pyi, label = 'py')
plt.legend()
plt.title('posicion en tiempo')
plt.xlabel('t')
plt.ylabel('p(t)')
plt.grid()

plt.show()

Ejercicios: 1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Se requiere determinar la distribución de cupos en base a los costos relativos al promedio por estudiante para docencia, infraestructura y servicios mostrados en la tabla.

Costo referencial /carrera	Mecatrónica	Computación	Civil	Matemáticas
Docencia	1.5	0.9	0.6	0.7
Infraestructura	0.8	1.4	0.4	0.5
Servicios	0.45	0.55	1.1	0.5

Con los datos del total de recursos relativos al promedio por estudiante disponibles son docencia 271, infraestructura 250 y servicios 230.

1.5 a + 0.9 b + 0.6 c + 0.7 d = 271 0.8 a + 1.4 b + 0.4 c + 0.5 d = 250 0.45 a + 0.55 b + 1.1 c + 0.5 d = 230

se indica que en carreras como matemáticas de baja demanda, se establece el cupo de 10,

1.5 a + 0.9 b + 0.6 c + 0.7 (10) = 271 0.8 a + 1.4 b + 0.4 c + 0.5(10) = 250 0.45 a + 0.55 b + 1.1 c + 0.5(10) = 230

el sistema se convierte en:

1.5 a + 0.9 b + 0.6 c = 271 - 0.7 (10) 0.8 a + 1.4 b + 0.4 c = 250 - 0.5(10) 0.45 a + 0.55 b + 1.1 c = 230 - 0.5(10)

Para usar un método iterativo se convierte a matriz aumentada:

\begin{pmatrix} 1.5 & 0.9 & 0.6 & \Big| & 264 \\ 0.8 & 1.4 & 0.4 & \Big| & 245 \\ 0.45 & 0.55 & 1.1 &\Big| & 225 \end{pmatrix}

con pivoteo parcial por filas, la matriz aumentada se mantiene igual, pues los valores de la diagonal ya son los mayores posibles según el algoritmo.

Para un método iterativo se despeja una ecuación por cada incógnita.

a = \frac{1}{1.5}(264 - 0.9 b - 0.6 c) b = \frac{1}{1.4}(245 - 0.8 a - 0.4 c) c = \frac{}{1.1}(225 -0.45 a - 0.55 b)

Para los valores iniciales se consideran números mayores que cero, pues existen recursos para los cupos. No se admiten cupos negativos.

X_0 = [50,50,50]

Las iteraciones para el método iterativo de Gauss-Seidel

itera = 0

a = \frac{1}{1.5}(264 - 0.9 (50) - 0.6 (50)) = 126 b = \frac{1}{1.4}(245 - 0.8 (126) - 0.4 (50)) = 88.714 c = \frac{}{1.1}(225 -0.45 (126) - 0.55 (88.714)) =108.642 diferencia = [126-50, 88.714-50, 108.42-50] diferencia = [76, 38.714, 58.642] errado = max|[76, 38.714, 58.642]| =76 X = [126, 88.714, 108.42]

itera = 1

a = \frac{1}{1.5}(264 - 0.9 (88.714) - 0.6 (108.42)) = 79.314 b = \frac{1}{1.4}(245 - 0.8 (79.314) - 0.4 (108.42)) = 98.637 c = \frac{}{1.1}(225 -0.45 (79.314) - 0.55 (98.637)) =122.780 diferencia = [79.314-126, 88, 98.637-88.714, 122.780-108.42] diferencia = [46.685,9.922, 14.137] errado = max| [46.685,9.922, 14.137] | = 46.685

el error disminuye en la iteración

X = [79.314 , 98.637, 122.780]

itera = 2

a = \frac{1}{1.5}(264 - 0.9 (79.314) - 0.6 (122.780)) = 67.705 b = \frac{1}{1.4}(245 - 0.8 (67.705) - 0.4 (122.780)) = 101.230 c = \frac{}{1.1}(225 -0.45 (67.705) - 0.55 (101.230)) =126.232 diferencia = [67.705-79.314, 101.230-98.637, 126.232-122.780] diferencia = [-11.608, 2.594, 3.451] errado = max| [-11.608, 2.594, 3.451] | = 11.608

el error disminuye en la iteración, se considera que el método converge

X = [67.705 , 101.230, 126.232]

con el algoritmo se tiene como resultado:

[126.          88.71428571 108.64285714]
[76.         38.71428571 58.64285714]

[ 79.31428571  98.63673469 122.78033395]
[46.68571429  9.92244898 14.13747681]

[ 67.7058256  101.23086138 126.23218611]
[11.60846011  2.59412669  3.45185216]

[ 64.76860873 101.92302755 127.08769174]
[2.93721688 0.69216617 0.85550564]

[ 64.01110677 102.11145563 127.30336487]
[0.75750196 0.18842808 0.21567312]

[ 63.81178067 102.16373537 127.3587675 ]
[0.1993261  0.05227973 0.05540263]

[ 63.75825178 102.17849398 127.37328637]
[0.05352889 0.01475862 0.01451887]

[ 63.74358906 102.18272443 127.37716953]
[0.01466272 0.00423045 0.00388316]

[ 63.73949753 102.18395297 127.37822907]
[0.00409153 0.00122854 0.00105954]

[ 63.73833659 102.18431364 127.37852366]
[0.00116094 0.00036067 0.0002946 ]

[ 63.73800235 102.18442047 127.37860699]
[3.34240232e-04 1.06824300e-04 8.33224905e-05]

respuesta X: 
[[ 63.73800235]
 [102.18442047]
 [127.37860699]]
verificar A.X=B: 
[[264.00014614]
 [245.00003333]
 [225.        ]]
>>>

se interpreta la respuesta como la parte entera de la solución:

cupos = [ 63, 102 , 127]

Ejercicio: 1Eva_2022PAOII_T1 Esfera flotando en agua

Según el principio de Arquímedes, la fuerza de flotación o empuje es igual al peso de el fluido desplazado por la porción sumergida de un objeto.

F_{empuje} = F_{peso} \rho_{agua} V_{sumergido} \text{ } g = \rho_{esfera}V_{esfera} \text{ } g V_{sumergido} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{esfera} - V_{sobreagua} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{sobreagua} = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) V_{esfera} V_{esfera} = \frac{4}{3}\pi r^3 \frac{\pi h^2}{3}(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) \frac{4}{3}\pi r^3 h^2(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3

El planteamiento para la búsqueda de raíces es f(x) = 0, que para este caso será:

f(h) = h^2(3r-h) - \Big( 1 - \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3 = 0

usando los valores dados para el ejercicio, r=1 y ρ_esfera = 200 Kg/m³ y ρ_agua    = 1000 kg/m³ se tiene que:

f(h) = h^2(3-h) - \Big( 1 - \frac{200}{1000}\Big) 4 f(h) = h^2(3-h) - \frac{16}{5}

Se observa la gráfica de f(h) en el intervalo de h entre[0,2] interpretado como totalmente sumergida y totalmente flotando sobre el agua, confirmando que existe una raíz

Para el caso de aplicar el método del punto fijo se plantea que x=g(x),

h = g(h) h^2(3-h) = \frac{16}{5}

con lo que se puede plantear dos ecuaciones al despejar h

h = \sqrt{ \frac{16}{5(3-h)}} h = 3-\frac{16}{5 h^2}

Iteraciones de la primera ecuación

itera = 0 ; h = h₀ = 0.5 ;

g(h) = \sqrt{ \frac{16}{5(3-0.5)}} = 1.1313 tramo = |1.1313-0.5|=0.6313

itera = 1 ; h = 1.1313 ;

g(h) = \sqrt{ \frac{16}{5(3-1.1313)}} = 1.3086 tramo = |1.3086-1.1313| = 0.1772

itera = 2 ; h = 1.3086 ;

g(h) = \sqrt{ \frac{16}{5(3-1.3086)}} = 1.3754 tramo = |1.3754-1.3086| = 0.0668

Observando los errores o tramos en cada iteración se tiene que se reduce, el método converge.

---

resultados.txt

x,g(x),tramo
0.5 1.131370849898476 0.631370849898476
1.131370849898476 1.308619626317284 0.17724877641880799
1.308619626317284 1.3754802083033437 0.06686058198605971
1.3754802083033437 1.4035002223557855 0.02802001405244181
1.4035002223557855 1.4157629993958152 0.012262777040029649
1.4157629993958152 1.4212317895316 0.005468790135784829
1.4212317895316 1.4236912066694054 0.0024594171378053975
1.4236912066694054 1.424801422465215 0.0011102157958096104
1.424801422465215 1.4253034412081806 0.0005020187429656264
raiz: 1.4253034412081806

Algoritmo en Python

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    print('x,g(x),tramo')
    print(a,b,tramo)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        print(a,b,tramo)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------
# INGRESO
fx = lambda h: h**2*(3-h)-16/5
gx = lambda h: np.sqrt(16/(5*(3-h)))

#fx = lambda h: h**2*(3-h)-16/5
#gx = lambda h: 3-16/(5*(h**2))

x0 = 0.5
tolera = 0.001
iteramax = 50  # itera máximo
a = 0     # intervalo
b = 2
muestras = 51  # gráfico

# PROCEDIMIENTO
respuesta = puntofijo(gx,x0,tolera)

# SALIDA
print('raiz:',respuesta)

hi = np.linspace(a,b,muestras)
fi = fx(hi)
gi = gx(hi)
plt.plot(hi,fi,label='f(h)')
plt.plot(hi,gi,label='g(h)')
plt.plot(hi,hi,label='Identidad')
plt.axhline(0,color='grey')
plt.grid()
plt.xlabel('h')
plt.ylabel('f(h)')
plt.title('esfera sumergida')
plt.legend()
plt.show()