Categoría: Solución 2da Eva

Ejercicio: 2Eva_2024PAOII_T2 Mayoría entre grupos Azules y Rojos

literal a

La población anual del país se describe con x(t), con tasas de natalidad a = 0.018 y mortalidad b = 0.012,

f(t,x,,y) = \frac{\delta x}{\delta t} = 0.018 x - 0.012 x^2

La población de Rojos es minoría y se describe con y(t). Los valore iniciales son: x(0)=2 , y(0)=0.5 y tamaño de paso h=0.5

g(t,x,y) = \frac{\delta y}{\delta t} = 0.026x - 0.017 y^2 +0.19 (0.012) (x-y)

el algoritmo de Runge-Kutta para sistemas de ecuaciones aplicado al ejercicio:

K1x = h f(t,x,y) = h (0.018 x - 0.012 x^2) K1y = h g(t,x,y) = h \Big(0.026x - 0.017 y^2 +0.19 (0.012) (x-y)\Big) K2x = h f(t+h,x+K1x,y+K1y) = h (0.018 (x+K1x) - 0.012 (x+K1x)^2) K2y = h g(t+h,x+K1x,y+K1y) = h \Big(0.026(x+K1x) - 0.017 (y + K1y)^2 +0.19 (0.012) ((x+K1x)-(y + K1y))\Big) x[i+1] = x[i] + \frac{K1x+K2x}{2} y[i+1] = y[i] + \frac{K1y+K2y}{2} t[i+1] = t[i] + h

literal b

Desarrolle tres iteraciones con expresiones completas para x(t), y(t) con tamaño de paso h=0.5.

itera = 0

K1x = 0.5 (0.018 (2) - 0.012 (2)^2) = -0.006 K1y = 0.5 \Big(0.026(2) - 0.017 (0.5)^2 +0.19 (0.012) ((2)-(0.5))\Big) = 0.006085 K2x = 0.5 (0.018 (2-0.006) - 0.012 (2-0.006)^2) = -0.00591 K2y = 0.5 \Big(0.026(2-0.006) - 0.017 (0.5 + 0.006085)^2 +0.19 (0.012) ((2-0.006)-(0.5 + 0.006085))\Big) = 0.006098 x[1] = 2 + \frac{-0.006+-0.00591}{2} = 1.9940 y[1] = 0.5 + \frac{0.006085+0.006098}{2} = 0.5060 t[1] = 0 + 0.5 =0.5

itera = 1

K1x = 0.5 (0.018 (1.9940) - 0.012 (1.9940)^2) =-0.005911 K1y = 0.5 \Big(0.0261.9940 - 0.017 (0.5060)^2 +0.19 (0.012) (1.9940-y)\Big) = 0.006098 K2x = 0.5 (0.018 (1.9940-0.005911) - 0.012 (1.9940-0.005911)^2 = -0.005823 K2y = 0.5 \Big(0.026(1.9940-0.005911) - 0.017 (0.5060 + 0.006098)^2 +0.19 (0.012) ((1.9940-0.005911)-(0.5060 + 0.006098))\Big) = 0.006111 x[2] = 1.9940 + \frac{-0.005911-0.005823}{2} = 1.988178 y[2] = 0.5060 + \frac{0.006098+0.006111}{2} = 0.512196 t[2] = 0.5 + 0.5 = 1

itera = 2

K1x = 0.5 (0.018 (1.988178) - 0.012 (1.988178)^2) =-0.005824 K1y = 0.5 \Big(0.026(0.512196) - 0.017 (0.512196)^2 +0.19 (0.012) (1.988178-0.512196)\Big) = 0.006111 K2x = 0.5 (0.018 (1.988178-0.005824) - 0.012 (1.988178-0.005824)^2) =-0.005737 K2y = 0.5 \Big(0.026(0.512196+0.006111) - 0.017 (0.512196 + 0.006111 )^2 +0.19 (0.012) ((1.988178-0.005911)-(0.512196 + 0.006111)\Big) =0.006124 x[3] = (1.988178) + \frac{-0.005824-0.005737}{2} = 1.982398 y[3] = 0.512196 + \frac{0.006111+0.006124}{2} = 0.518314 t[3] = 1 + 0.5 = 1.5

literal c

Los resultados para x(t) de las primeras iteraciones indican que la población total del país disminuye en el intervalo observado. Se puede comprobar al usar el algoritmo para mas iteraciones como se muestra en la gráfica.

literal d

Los resultados para y(t) muestran que la población clasificada como Rojo aumenta en el intervalo observado. Sin embargo la mayoría sigue siendo Azul para las tres iteraciones realizadas.

literal e

La población y(t) alcanza la mitad de la población cuanto t cambia de 30.5 a 31. Tiempo en el que de realizar elecciones, ganarían los Rojos.

Usando el algoritmo, se añade una columna con la condición que MayoriaY = yi>xi/2, que se convierte a 1 o verdadero cuando se cumple la condición. Con lo que se encuentra el cambio de mayoría a Rojo.

Runge-Kutta Segundo Orden
i  [ ti,  xi,  yi ]
   [ K1x,  K1y,  K2x,  K2y , MayoriaY]
0 [0.  2.  0.5]
   [0. 0. 0. 0. 0.]
1 [0.5      1.994045 0.506092]
  [-0.006     0.006085 -0.00591   0.006098  0.      ]
2 [1.       1.988178 0.512196]
  [-0.005911  0.006098 -0.005823  0.006111  0.      ]
3 [1.5      1.982398 0.518314]
  [-0.005824  0.006111 -0.005737  0.006124  0.      ]
4 [2.       1.976702 0.524443]
…
   [-0.002761  0.005927 -0.002727  0.005907  0.      ]
60 [30.        1.755801  0.869617]
   [-0.002728  0.005907 -0.002695  0.005887  0.      ]
61 [30.5       1.753123  0.875494]
   [-0.002695  0.005887 -0.002662  0.005867  0.      ]
62 [31.        1.750476  0.88135 ]
   [-0.002663  0.005867 -0.002631  0.005846  1.      ]
63 [31.5       1.747861  0.887185]
   [-0.002631  0.005846 -0.002599  0.005824  1.      ]
64 [32.        1.745277  0.892998]
   [-0.002599  0.005824 -0.002568  0.005802  1.      ]
65 [32.5       1.742724  0.898789]
   [-0.002568  0.005802 -0.002538  0.00578   1.      ]
66 [33.        1.740201  0.904558]
   [-0.002538  0.00578  -0.002508  0.005757  1.      ]
67 [33.5       1.737708  0.910303]
   [-0.002508  0.005757 -0.002478  0.005734  1.      ]

 

Ejercicio: 2Eva_2024PAOII_T1 Área de incendio forestal en Cerro Azul

literal a

Calcular los tamaños de paso dxi en cada frontera y plantear la integración con fórmulas compuestas.

Usando los datos de las coordenadas de obtiene cada dxi = xi[i+1]-xi[i]. De forma semejante se encuentra cada dxj, Seleccionando los métodos según se disponga de tamaño de paso iguales y consecutivos como se muestra en la tabla ampliada.

Frontera superior

		Trapecio		Trapecio			Trapecio		Trapecio			Simpson 1/3
	Trapecio		Trapecio		Simpson 1/3			Trapecio		Simpson 1/3
dxi	40	100	-30	66	20	20	79	125	54	50	50	20	20	—
xi	410	450	550	520	586	606	626	705	830	884	934	984	1004	1024
yi	131	194	266	337	402	483	531	535	504	466	408	368	324	288

Frontera Inferior

				Trapecio
	Simpson 3/8
dxj	190	190	190	44	—
xj	410	600	790	980	1024
yj	131	124	143	231	288

Desarrollando con instrucciones sobre el arreglo en Python con la instrucción np.diff(xi).

>>> xi = [410, 450, 550, 520, 586, 606, 626, 705, 830,
          884, 934, 984, 1004, 1024]
>>> xi = np.array(xi)
>>> dxi = np.diff(xi)
>>> dxi
array([ 40, 100, -30, 66, 20, 20, 79, 125, 54,
        50, 50, 20, 20])
>>> xj = [410, 600, 790, 980, 1024]
>>> dxj = np.array(xj)
>>> dxj = np.diff(xj)
>>> dxj
array([190, 190, 190, 44])
>>>

literal b

Desarrollar las expresiones del área para las coordenadas de la frontera superior, según el literal a. Cuando se tienen dos tamaños de paso iguales se usa Simpson de 1/3.

I_{superior} = 40\Big(\frac{131+194}{2}\Big) +100\Big(\frac{194+266}{2}\Big) + -30\Big(\frac{266+337}{2}\Big) +66\Big(\frac{337+402}{2}\Big)+ + \frac{20}{3}\Big(402+4(483)+531\Big) + +79\Big(\frac{531+535}{2}\Big) +125\Big(\frac{535+504}{2}\Big) + 54\Big(\frac{504+466}{2}\Big) + + \frac{50}{3}\Big(466+4(408)+368\Big) + + \frac{20}{3}\Big(368+4(324)+288\Big) I_{superior} = 254753,33

literal c

Realice los cálculos para la frontera inferior y encuentre el área afectada. Con tres tamaños de paso iguales se usa Simpson de 3/8.

I_{Inferior} = \frac{3}{8}(190)\Big(131+3(124)+3(143)+231\Big) +44\Big(\frac{231+288}{2}\Big) = 94281,75

Área total afectada:

A_{afectada} = I_{superior} - I_{Inferior} = 254753,33 -94281,75 A_{afectada} = 160.471,58

literal d

Estime la cota de error en los cálculos.

Considere usar las unidades en Km en lugar de metros para los tamaños de paso.

Error_{truncaSup} = O( 0.040^3) + O(0. 100^3)+ O( (-0.030)^3) + O( 0.066^3)+ + O( 0.020^5) + O(0. 079^3)+ O( 0.125^3) + O( 0.054^3)+ + O( 0.050^5) + O( 0.020^5) Error_{truncaInf} = O( 0.190^5) + O(0. 044^3)

Ejercicio: 2Eva_2024PAOI_T1 EDO Salto Bungee tiempo extiende cuerda

Como referencia par el ejercicio, solo se realizará la primera parte, acorde a:

Encuentre el tiempo t_c y la velocidad de la persona cuando se alcanza la longitud de la cuerda extendida y sin estirar (30 m), es decir y<L, aún se entra cayendo signo(v)=1. (solo primera ecuación)

\frac{d^2y}{dt^2} = g - signo(v) \frac{C_d}{m}\Big( \frac{dy}{dt}\Big)^2 y \leq L

Suponga que las condiciones iniciales son:

y(0) = 0

\frac{dy(0)}{dt} = 0

literal a

Realice el planteamiento del ejercicio usando Runge-Kutta de 2do Orden

\frac{d^2y}{dt^2} = 9.8 - signo(v) \frac{0.25}{68.1}\Big( \frac{dy}{dt}\Big)^2

Usando los valores de las constantes g, C_d, m, haciendo el cambio de variable dy/dt=v. Adicionalmente, en la caída inicial, signo(v)=1 y se mantiene constante hasta el 2do tramo con la cuerda estirada.

v' = 9.8 - \frac{0.25}{68.1} v^2

Se plantea el ejercicio como Runge-Kutta para 2da derivada

v = y' = f(t,y,v) v' = g(t,y,v) = 9.8 - \frac{0.25}{68.1} v^2

siendo h=0.5, las iteraciones se realizan como:

K_{1y} = 0.5 v K_{1v} = 0.5 \Big( 9.8 - \frac{0.25}{68.1} v^2\Big) K_{2y} = 0.5 \Big( v_i + K_{1v}\Big) K_{2v} = 0.5 \Big(9.8 - \frac{0.25}{68.1}(v_i + K_{1v})^2 \Big) y_{i+1}=y_i+\frac{K_{1y}+K_{2y}}{2} v_{i+1}=v_i+\frac{K_{1v}+K_{2v}}{2} t_{i+1} = t_i +0.5

literal b

Desarrolle tres iteraciones para y(t) con tamaño de paso h=0.5

itera = 0

K_{1y} = 0.5 (0) = 0 K_{1v} = 0.5 \Big( 9.8 - \frac{0.25}{68.1} (0)^2\Big) = 4.9 K_{2y} = 0.5 ( 0 + 4.9) = 2.45 K_{2v} = 0.5 \Big(9.8 - \frac{0.25}{68.1}(0 + 2.45)^2 \Big) =4.8559 y_{i+1}=0+\frac{0+2.45}{2}=1.225 v_{i+1}=0+\frac{4.9+4.8889}{2}=4.8779 t_{i+1} =0 +0.5 = 0.5

itera=1

…

Continuar con las iteraciones como tarea.

literal c

Usando el algoritmo, aproxime la solución para y en el intervalo entre [0,t_c], adjunte sus resultados.txt

las iteraciones con el algoritmo son:

 [ ti, yi, vi,K1y,K1v,K2y,K2v]
[[0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
  0.000000e+00 0.000000e+00]
 [5.000000e-01 1.225000e+00 4.877964e+00 0.000000e+00 4.900000e+00
  2.450000e+00 4.855929e+00]
 [1.000000e+00 4.878063e+00 9.669162e+00 2.438982e+00 4.856324e+00
  4.867144e+00 4.726071e+00]
 [1.500000e+00 1.089474e+01 1.429311e+01 4.834581e+00 4.728391e+00
  7.198776e+00 4.519513e+00]
 [2.000000e+00 1.917255e+01 1.868062e+01 7.146557e+00 4.525013e+00
  9.409063e+00 4.249997e+00]
 [2.500000e+00 2.957773e+01 2.277738e+01 9.340309e+00 4.259461e+00
  1.147004e+01 3.934054e+00]
 [3.000000e+00 4.195334e+01 2.654573e+01 1.138869e+01 3.947708e+00
  1.336254e+01 3.589005e+00]

con lo que se observa que para alcanzar los 30m de la cuerda sin estirar se alcanzan entre los [2.5, 3.0] s.

literal d

Indique el valor de t_c, muestre cómo mejorar la precisión y realice sus observaciones sobre los resultados.

Para mejorar la precisión del algoritmo se reduce el valor del tamaño de paso h,  de esta forma se puede obtener una mejor lectura del tiempo de la primera fase del ejercicio.

Por ejemplo haciendo el tamaño de paso h=0.5/4, el tiempo entre [2.5, 2.625]. que tiene un error segundos menor al valor encontrado inicialmente y se mejora la precisión.

 

Ejercicio: 2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B

Literal a

La ecuación diferencial a resolver es:

\frac{d^2y}{dx^2} = \frac{w_0}{T_0} \Big[ 1+ \sin \Big(\frac{\pi x}{2l_B} \Big) \Big]

donde w₀ = 1 000 lbs/ft, T₀. = 0.588345×10⁶.
dy(0)/dx = 0 y l_B=200 de la gráfica presentada.

\frac{d^2y}{dx^2} = \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

Para usar Runge Kutta para segunda derivada:

z= y' = f(x,y,z) z' = (y')' = 0z + \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big] g(x,y,z) = \frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

los valores iniciales para el ejercicio acorde al enunciado son: x₀ = 0, y₀=0, z₀ = 0, con h=0.5

literal b

para itera 0

K1y = h*z = 0.5*0 = 0 K1z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0)}{2(200)} \Big) \Big] = 0.0008498 K2y = h*(z+K1z) = (0.5) (0+0.00084984) = 0.0004249 K2z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0+0.5)}{2(200)} \Big) \Big] =0.0008531 y = 0+\frac{0+0.0004249}{2} = 0.0002124 z = 0+\frac{0.0008498+0.0008531}{2} = 0.0008515 x = 0 + 0.5 = 0.5

…

Desarrollar dos iteraciones adicionales como tarea.

Para las primeras iteraciones de un total de 400+1, los valores con Python y en resultados.txt :

estimado[xi,yi,zi,K1y,K2y,K1z,K2z]
[0.0000 0.0000e+00 0.0000e+00 
0.0000e+00 0.0000e+00 
0.0000e+00 0.0000e+00]

[0.5000 2.124603761398499073e-04 8.515101601627471980e-04 
0.000000000000000000e+00 4.249207522796998146e-04 
8.498415045593996292e-04 8.531788157660947667e-04]

[1.0000 8.515101601627470896e-04 1.706357605799455881e-03 
4.257550800813735990e-04 8.523444879644209281e-04 
8.531788157660947667e-04 8.565160755073224913e-04]

[1.5000 1.918817981939305653e-03 2.564542259712321911e-03 
8.531788028997279406e-04 1.281436840653389295e-03 
8.565160755073224913e-04 8.598532323184093513e-04]

[2.0000 3.416052419875068892e-03 3.426063993239661116e-03 
1.282271129856160955e-03 1.712197746015365740e-03 
8.598532323184093513e-04 8.631902347362692754e-04]
...

literal c

resultado en archivo.txt al ejecutar el algoritmo.

literal d

Algoritmo con Python

# 2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3+4),dtype=float)

    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K2y,K1z,K2z]
    return(estimado)

# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
T0 = 0.588345e6
LB = 200
f = lambda x,y,z: z
g = lambda x,y,z: (1000/T0)*(1+np.sin(np.pi*x/(2*LB)))
x0 = 0
y0 = 0
z0 = 0
h  = 0.5
muestras = 401

# PROCEDIMIENTO
puntosRK2 = rungekutta2_fg(f,g,x0,y0,z0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi,zi,K1y,K2y,K1z,K2z]')
print(puntosRK2)
np.savetxt("tablaRk2.txt",puntosRK2)


# Gráfica
import matplotlib.pyplot as plt

plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         color='m',
         label ='y Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

 

Ejercicio: 2Eva_2023PAOII_T1 Volumen por solido de revolución

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

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

literal a y c

Para el volumen con f(x) con al menos 3 tramos y un método de Simpson, directamente se puede usar 3/8. Por lo que se Se reemplaza en la fórmula de volumen del sólido de revolución f(x) con:

f(x) = \sqrt{\sin (x/2)}

obteniendo:

V_{fx} = \int_{a}^{b} \pi \Big(\sqrt{\sin (x/2)} \Big)^2 dx = \int_{a}^{b} \pi \sin (x/2) dx

La expresión dentro del integral se denomina como f_v:

f_v (x)= \pi \sin (x/2)

en el intervalo [0.1, 1.8],  con al menos 3 tramos, se requieren 4 muestras con tamaño de paso h_f: y truncando a 4 decimales los resultados calculados con Python.

h_f =\frac{b-a}{tramos} = \frac{1.8-0.1}{3} = 0.5666

los puntos de muestra quedan np.linspace(0.1,1.8,3+1):

x_is= [0.1, 0.6666, 1.2333, 1.8 ]

El integral se calcula con los puntos de muestra,

V_{fx} = \frac{3}{8} (0.5666) \Big( f_v(0.1) +3 f_v(0.6666) + + 3 f_v(1.2333)+ f_v(1.8)\Big)

recordando que se usa en radianes,

V_{fx} = \frac{3}{8} (0.5666) \Bigg( \pi \sin \Big(\frac{0.1}{2}\Big) +3 \pi \sin \Big(\frac{0.6666}{2}\Big) + + 3 \pi \sin\Big(\frac{1.2333}{2}\Big)+ \pi \sin \Big(\frac{1.8}{2}\Big)\Bigg) = \frac{3}{8} (0.5666) \Big( 0.1570+3 (1.0279) + + 3 (1.8168)+ 2.4608\Big)

literal d. el volumen generado por f(x) tiene como resultado:

V_{fx} = 2.3698

la cota de error para fx es el orden de O(h⁵) = O(0.5666⁵) = O(0.05843), queda como tarea completar la cota de error total.

literal b y c

Para el volumen con g(x) con al menos 2 tramos y Cuadratura de Gauss de dos puntos, se reemplaza en la fórmula de volumen de sólido de revolución:

g(x) = e^{x/3} - 1 V_{gx} = \int_{a}^{b} \pi (e^{x/3} - 1)^2 dx

La expresión dentro del integral se denomina como g_v:

g_v = \pi (e^{x/3} - 1)^2

en el intervalo [0.1, 1.8],  con al menos 2 tramos, se requieren 3 muestras con tamaño de paso h_g:

h_g =\frac{b-a}{tramos} = \frac{1.8-0.1}{2} = 0.85

x_ic = [0.1, 0.95, 1.8 ]

tramo 1: [0.1, 0.95] , a = 0.1 , b= 0.95, truncando a 4 decimales

x_a = \frac{0.95+0.1}{2} - \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.2796 x_b = \frac{0.95+0.1}{2} + \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.7703 g_v(0.2796) = \pi (e^{0.2796/3} - 1)^2 = 0.02998 g_v(0.7703) = \pi (e^{0.7703/3} - 1)^2 = 0.2692 V_{c1} = \frac{0.95-0.1}{2}(g_v(0.2796) + g_v(0.7703)) V_{c1} = \frac{0.95-0.1}{2}(0.02998 + 0.2692) V_{c1} = 0.1271

tramo 2: [0.95, 1.8] , a = 0.95 , b= 1.8

x_a = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.1296 x_b = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.6203 g_v(1.1296) = \pi (e^{1.1296/3} - 1)^2 = 0.6567 g_v(1.6203) = \pi (e^{1.6203/3} - 1)^2 = 1.6115 V_{c2} = \frac{1.8-0.95}{2}(g_v(1.1296) + g_v(1.6203)) V_{c2} = \frac{1.8-0.95}{2}(0.6567 + 1.6115) V_{c2} = 0.9640

literal d. volumen generado por g(x)

V_{gx} = V_{c1} + V_{c2} = 0.1271 + 0.9640 = 1.0912

completar la cota de error para cuadratura de Gauss de dos puntos.

literal e. El volumen de revolución se genera como la resta del volumen de f(x) y volumen g(x)

V = V_{fx} - V_{gx} = 2.3698 - 1.0912 = 1.2785

---

Algoritmo con Python

Los resultados usando el algoritmo con las operaciones usadas en el planteamiento son:

para f(x):
xis= [0.1        0.66666667 1.23333333 1.8       ]
fiv= [0.15701419 1.02791246 1.81684275 2.46089406]
Volumenfx:  2.369836948864926

para g(x):
Por tramos: [0.1  0.95 1.8 ]
xab= [0.2796261355944091, 0.770373864405591,
      1.129626135594409, 1.620373864405591]
gab= [0.02998177327492598, 0.26928904479784566,
      0.6567986343358181, 1.6115494735531555]
Vc1= 0.12719009768092793  ; Vc2= 0.964047945852814
Volumengx:  1.0912380435337419

Volumen solido revolucion: 1.2785989053311841

Considerando realizar los cálculos para cada sección:

# 2Eva_2023PAOII_T1 Volumen por solido de revolución
import numpy as np

# INGRESO
fx = lambda x: np.sqrt(np.sin(x/2))
gx = lambda x: np.exp(x/3)-1
a = 0.1
b = 1.8
tramosSimpson = 3
tramosCGauss = 2

# PROCEDIMIENTO
# Volumen para f(x) con Simpson
fv = lambda x: np.pi*np.sin(x/2)
hs = (b-a)/tramosSimpson
xis = np.linspace(a,b,tramosSimpson +1)
fiv = fv(xis)
Vs = (3/8)*hs*(fiv[0]+3*fiv[1]+3*fiv[2]+ fiv[3])

# Volumen para g(x) con Cuadratura de Gauss
gv = lambda x: np.pi*(np.exp(x/3)-1)**2
hc = (b-a)/tramosSimpson
xic = np.linspace(a,b,tramosCGauss +1)
# tramo 1
ac = xic[0]
bc = xic[1]
xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
Vc1 = (bc-ac)/2*(gv(xa)+gv(xb))
xab = [xa,xb]
gab = [gv(xa),gv(xb)]
# tramo 2
ac = xic[1]
bc = xic[2]
xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
Vc2 = (bc-ac)/2*(gv(xa)+gv(xb))
Vc = Vc1+Vc2
xab.append(xa)
xab.append(xb)
gab.append(gv(xa))
gab.append(gv(xb))

# Volumen solido revolucion
Volumen = Vs - Vc

# SALIDA
print("para f(x):")
print("xis=", xis)
print("fiv=", fiv)
print("Volumenfx: ",Vs)
print()
print("para g(x):")
print("Por tramos:",xic)
print("xab=", xab)
print("gab=", gab)
print("Vc1=",Vc1," ; Vc2=",Vc2) 
print("Volumengx: ",Vc)
print()
print("Volumen solido revolucion:",Volumen)

para la gráfica presentada en el enunciado (no requerida) , se complementa con las instrucciones:

# para grafica -------------------
import matplotlib.pyplot as plt
muestras = 21 # grafica
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
xig = np.linspace(a,b,tramosCGauss+1)
fis = fx(xis)
gig = gx(xig)

# grafica
plt.plot(xi,fi, label="f(x)")
plt.plot(xi,gi, label="g(x)")
plt.plot([0.0,2.0],[0,0], marker=".", color="blue")
plt.fill_between(xi,fi,gi,color="lightgreen")
plt.axhline(0)
plt.axvline(a, linestyle="dashed")
plt.axvline(b, linestyle="dashed")
plt.xlabel('x')
plt.ylabel('f(x), g(x)')
plt.legend()
plt.plot(xis,fis,'.b')
plt.plot(xig,gig,'.r')
plt.grid()
plt.show()

---

Gráfica de sólido de revolución en 3D

Instrucciones en Python

# 2Eva_2023PAOII_T1 Volumen por solido de revolución
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO
f = lambda x: np.sqrt(np.sin(x/2))
g = lambda x: np.exp(x/3)-1

# eje x
xa = 0.1
xb = 1.8
xmuestras = 31
# angulo w de rotación
w_a = 0
w_b = 2*np.pi
w_muestras = 31

# PROCEDIMIENTO
# muestreo en x y angulo w
xi = np.linspace(xa, xb, xmuestras)
wi = np.linspace(w_a, w_b, w_muestras)
X, W = np.meshgrid(xi, wi)

# evalua f(x) en 3D
Yf = f(xi)*np.cos(W)
Zf = f(xi)*np.sin(W)

# evalua g(x) en 3D
Yg = g(xi)*np.cos(W)
Zg = g(xi)*np.sin(W)

# SALIDA

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

grafica.plot_surface(X, Yf, Zf,
                     color='blue', label='f(x)',
                     alpha=0.3, rstride=6, cstride=12)
grafica.plot_surface(X, Yg, Zg,
                     color='orange', label='g(x)',
                     alpha=0.3, rstride=6, cstride=12)

grafica.set_title('Solidos de revolución')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
# grafica.legend()
eleva = 45
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: 2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

Dada la EDP elíptica,

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = \Big( x^2 + y^2 \Big) e^{xy} 0 \lt x \lt 1 0 \lt y \lt 0.5

Se convierte a la versión discreta usando diferencias divididas centradas, según se puede mostrar con la gráfica de malla.

literal b

literal a

---

\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = \Big( x^2 + y^2 \Big) e^{xy}

---

literal c

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy²

\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = \Big( x^2 + y^2 \Big) e^{xy}\frac{\Delta y^2}{\Delta x^2}

---

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1

---

se simplifica el coeficiente en λ =1

u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

---

Se agrupan los términosiguales

u[i-1,j]-4u[i,j]+ u[i+1,j] + u[i,j-1] +u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

---

Se desarrollan las iteraciones para tres rombos y se genera el sistema de ecuacioens a resolver.
para i=1,j=1

u[0,1]-4u[1,1]+ u[2,1] + u[1,0] +u[1,2] = \Big( x[1]^2 + y[1]^2 \Big) e^{x[1]y[1]} 1-4u[1,1]+ u[2,1] + 1 + \frac{1}{8} = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} -4u[1,1]+ u[2,1] = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} - \frac{1}{8}

para i=2, j=1

u[1,1]-4u[2,1]+ u[3,1] + u[2,0] +u[2,2] = \Big( x[2]^2 + y[1]^2 \Big) e^{x[2]y[1]} u[1,1]-4u[2,1]+ u[3,1] + 1 + 0.25 = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} u[1,1]-4u[2,1]+ u[3,1] = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} -1.25

para i=3, j=1

u[2,1]-4u[3,1]+ u[4,1] + u[3,0] +u[3,2] = \Big( x[3]^2 + y[1]^2 \Big) e^{x[3]y[1]} u[2,1]-4u[3,1]+ 0.25 + 0 + \frac{3}{8} = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} u[2,1]-4u[3,1] = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} - 0.25 - \frac{3}{8}

con lo que se puede crear un sistema de ecuaciones y resolver el sistema para cada punto desconocido

\begin{pmatrix} -4 & 1 & 0 & \Big| & 0.008061807364732415 \\ 1 & -4 & 1 & \Big| & -0.8958911084166168 \\0 & 1 & -4 &\Big| & 0.1288939058881129 \end{pmatrix}

se obtiene los resultados para:

u[1,1] = 0.05953113
u[2,1] = 0.24618634
u[3,1] = 0.02932311

>>> import numpy as np
>>> (0.25**2+0.25**2)*np.exp(0.25*0.25) - 1/8
0.008061807364732415
>>> (0.5**2+0.25**2)*np.exp(0.5*0.25) - 1.25
-0.8958911084166168
>>> (0.75**2+0.25**2)*np.exp(0.75*0.25) - 0.25 -3/8
0.1288939058881129
>>> A=[[-4,1,0],[1,-4,1],[0.0,1.0,-4.0]]
>>> B = [0.008061807364732415, -0.8958911084166168, 0.1288939058881129]
>>> np.linalg.solve(A,B)
array([0.05953113, 0.24618634, 0.02932311])

Ejercicio: 2Eva_2023PAOI_T2 Péndulo vertical amortiguado

literal a

\frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)

Se simplifica su forma a:

\frac{d\theta}{dt}= z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' = -\mu z -\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

se usan los valores dados: g = 9.81 m/s², L = 2 m

f_t(t,\theta,z) = z g_t(t,\theta,z) = - 0.5 z -\frac{9.81}{2}\sin (\theta)

y los valores iniciales para la tabla: θ(0) = π/4 rad, θ’ (0) = 0 rad/s, se complementan los valores en la medida que se aplica el desarrollo.

ti	θ(ti)	θ'(ti)=z
0	π/4	0
0.2	0.7161	-0.6583
0.4	0.5267	-0.1156
0.6	0.2579	-0.1410
…

literal b

Iteración 1:  ti = 0 ; yi = π/4 ; zi = 0

K1y = h * ft(ti,yi,zi) 
    = 0.2*(0) = 0
K1z = h * gt(ti,yi,zi) 
    = 0.2*(-0.5(0) -(9.81/2)sin (π/4) = -0.6930
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(0-0.6930)= -0.1386
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.5(0-0.6930) -(9.81/2)sin(π/4-0) 
    = -0.6237

yi = yi + (K1y+K2y)/2 
   = π/4+ (0+(-0.1386))/2 = 0.7161
zi = zi + (K1z+K2z)/2 
   = 0+(-0.6930-0.6237)/2 = -0.6583
ti = ti + h = 0 + 0.2 = 0.2

estimado[i] = [0.2,0.7161,-0.6583]

Iteración 2:  ti = 0.2 ; yi = 0.7161 ; zi = -0.6583

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-0.6583) = -0.1317
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5 ( -0.6583) -(9.81/2)sin (0.7161) 
    = -0.5775
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.6583 -0.5775)= -0.2472
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 (-0.6583 -0.5775) -(9.81/2)sin(0.7161-0.1317) 
    = -0.4171

yi = yi + (K1y+K2y)/2 
   = 0.7161 + (-0.1317-0.2472)/2 = 0.5267
zi = zi + (K1z+K2z)/2 
   = -0.6583+(-0.5775-0.4171)/2 = -0.1156
ti = ti + h = 0.2 + 0.2 = 0.4

estimado[i] = [0.4,0.5267,-0.1156]

Iteración 3:  ti = 0.4 ; yi = 0.5267 ; zi = -1.156

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-1.156) = -0.2311
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5(-1.156) -(9.81/2)sin (0.5267) 
    = -0.3771
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-1.156 -0.3771)= -0.3065
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 ( -1.156 -0.3771) -(9.81/2)sin(0.5267-0.2311) 
    = -0.1322

yi = yi + (K1y+K2y)/2 
   = 0.5267 + (-0.2311-0.3065)/2 = 0.2579
zi = zi + (K1z+K2z)/2 
   = -1.156+(-0.3771-0.1322)/2 = -1.410
ti = ti + h = 0.4 + 0.2 = 0.6

estimado[i] = [0.6,0.2579,-1.410]

literal c

resultados del algoritmo:

[ t, 		 y, 	 dyi/dti=z,  K1y,	 K1z,	    K2y,      K2z]
[[ 0.000e+00  7.854e-01  0.000e+00  0.000e+00  0.000e+00  0.000e+00   0.000e+00]
 [ 2.000e-01  7.161e-01 -6.583e-01  0.000e+00 -6.930e-01 -1.386e-01  -6.237e-01]
 [ 4.000e-01  5.267e-01 -1.156e+00 -1.317e-01 -5.775e-01 -2.472e-01  -4.171e-01]
 [ 6.000e-01  2.579e-01 -1.410e+00 -2.311e-01 -3.771e-01 -3.065e-01  -1.322e-01]
 [ 8.000e-01 -3.508e-02 -1.377e+00 -2.820e-01 -1.089e-01 -3.038e-01    1.756e-01]
...

con h=0.2 se tienen 1/0.2 = 5 tramos por segundo, por lo que para 10 segundo serán 50 tramos. La cantidad de muestras = tramos + 1(valor inicial) = 51

con lo que se puede usar el algoritmo en EDO Runge-Kutta d2y/dx2

literal d

Se observa que la respuesta es oscilante y amortiguada en magnitud como se esperaba según el planteamiento. Con el tiempo se estabilizará en cero.

Instrucciones en Python

# 2Eva_2023PAOI_T2 Péndulo vertical amortiguado
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,7),dtype=float)
    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
    return(estimado)

# INGRESO theta = y
g = 9.8
L  = 2
ft = lambda t,y,z: z
gt = lambda t,y,z: -0.5*z +(-g/L)*np.sin(y)

t0 = 0
y0 = np.pi/4
z0 = 0
h = 0.2
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)

# SALIDA
np.set_printoptions(precision=3)
print(' [ t, \t\t y, \t dyi/dti=z, K1y,\t K1z,\t K2y,\t K2z]')
print(tabla)

# Grafica
ti = np.copy(tabla[:,0])
yi = np.copy(tabla[:,1])
zi = np.copy(tabla[:,2])
plt.subplot(121)
plt.plot(ti,yi)
plt.grid()
plt.xlabel('ti')
plt.title('yi')
plt.subplot(122)
plt.plot(ti,zi, color='green')
plt.xlabel('ti')
plt.title('dyi/dti')
plt.grid()
plt.show()

Ejercicio: 2Eva_2023PAOI_T1 Material para medalla de academia

f(x) = 2-8\Big( \frac{1}{2} - x \Big)^2

0 \le x \lt 1 g(x) = -\Big( 1-x\Big)\ln \Big( 1- x \Big)

Para f(x) se usará Simpson de 1/3 que requiere al menos dos tramos para aplicar:

a. Realice el planteamiento de las ecuaciones para el ejercicio.

I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

b. Describa el criterio usado para determinar el número de tramos usado en cada caso.

h = \frac{b-a}{2} = \frac{1-0}{2} = 0.5

c. Desarrolle las expresiones completas del ejercicio para cada función.

I_{fx}\cong \frac{0.5}{3}[f(0)+4f(0.5) + f(1)] f(0) = 2-8\Big( \frac{1}{2} - (0) \Big)^2 = 0 f(0.5) = 2-8\Big( \frac{1}{2} - (0.5) \Big)^2 = 2 f(1) = 2-8\Big( \frac{1}{2} - (1) \Big)^2 = 0 I_{fx} = \frac{1}{6}[0+4(2) + 0] = \frac{8}{6} = \frac{4}{3} = 1.3333

cota de error O(h5) = O(0.5⁵)= O(0.03125)

Para g(x) se usará Simpson de 3/8 que requiere al menos tres tramos para aplicar:

I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] h = \frac{b-a}{3} = \frac{1-0}{3} = 0.3333 I_{gx}\cong \frac{3(0.3333)}{8}[f(0)+3f(0.3333) +3 f(0.6666)+f(1)] g(0) = -\Big( 1-0\Big)\ln \Big( 1- 0 \Big) = 0 g(0.3333) = -\Big( 1-0.3333\Big)\ln \Big( 1- 0.3333 \Big) = 0.2703 g(0.6666) = -\Big( 1-0.6666\Big)\ln \Big( 1- 0.6666 \Big) = 0.3662 g(0.9999) = -\Big( 1-0.9999\Big)\ln \Big( 1- 0.9999 \Big) = 0

para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.

I_{gx}\cong \frac{3(0.3333)}{8}[0+3(0.2703) + 3(0.3662)+0] = 0.2387

d. Indique el resultado obtenido para el área requerida y la cota de error
Area = I_{fx} – I_{gx} = 1.3333 – 0.2387 = 1.0945

cota de error = O(0.03125) + O(0.00411) = 0.03536

e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032

Si el factor de mayor error es de Simpson 1/3, se considera como primera aproximación que:

cota de error O(h⁵) = O(0.00032), h = (0.00032)^(1/5) = 0.2
es decir el número de tramos es de al menos (b-a)/tramos = 0.2 , tramos = 5.
El número de tramos debe ser par en Simpson de 1/3, por lo que se toma el entero mayor tramos=6 y el tamaño de paso recomendado es al menos 1/6. EL error al aplicar 3 veces la formula es 3(O((1/6)⁵)) = 0.0003858.

Lo que podría indicar que es necesario al menos dos tramos adicionales con h=1/8 y error O(0,00012) que cumple con el requerimiento.

Se puede aplicar el mismo criterio para Simpson 3/8 y se combinan los errores para verificar que cumplen con el requerimiento.

Algoritmo con Python

Resultados

Ifx:  1.3333332933359998
Igx:  0.238779092876627
Area:  1.094554200459373

Instrucciones en Python usando las funciones

# 2Eva_2023PAOI_T1 Material para medalla de academia
import numpy as np
import matplotlib.pyplot as plt

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

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)

# INGRESO
fx = lambda x: 2-8*(0.5-x)**2
gx = lambda x: -(1-x)*np.log(1-x)
a = 0
b = 1-1e-4
muestras1 = 2+1
muestras2 = 3+1

# PROCEDIMIENTO
xi1 = np.linspace(a,b,muestras1)
xi2 = np.linspace(a,b,muestras2)
fi = fx(xi1)
gi = gx(xi2)

Ifx = integrasimpson13_fi(xi1,fi)
Igx = integrasimpson38_fi(xi2,gi)
Area = Ifx - Igx

# SALIDA
print('Ifx: ', Ifx)
print('Igx: ', Igx)
print('Area: ', Area)

plt.plot(xi1,fi,'ob',label='f(x)')
plt.plot(xi2,gi,'or', label='g(x)')
plt.grid()
plt.legend()
plt.xlabel('xi')

# curvas suave con mas muestras (no en evaluación)
xi = np.linspace(a,b,51)
fxi = fx(xi)
gxi = gx(xi)
plt.fill_between(xi,fxi,gxi,color='navajowhite')
plt.plot(xi,fxi,color='blue',linestyle='dotted')
plt.plot(xi,gxi,color='red',linestyle='dotted')

plt.show()

Ejercicio: 2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

\frac{\partial^2 u}{\partial x^2} = b \frac{\partial u}{\partial t}

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = b\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

agrupando variables

\frac{\Delta t}{b} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = \frac{\Delta t}{b}b\frac{u_{i,j+1}-u_{i,j}}{\Delta t} λ = \frac{\Delta t}{b(\Delta x)^2} λ = \frac{0.002}{2(0.2)^2} =0.025

como λ<0.5 el método converge.

\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] = u[i,j+1]-u[i,j]

por el método explícito:

u[i,j+1] =\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] + u[i,j] u[i,j+1] =\lambda u[i+1,j]+(1-2\lambda)u[i,j]+\lambda u[i-1,j]

iteración i=1, j=0

u[1,1] =\lambda u[0,0]+(1-2\lambda)u[1,0]+\lambda u[2,0] u[1,1] =0.025(1) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.2\Big)+0.025 \cos \Big( \frac{3π}{2}0.4\Big)

iteración i=2, j=0

u[2,1] =\lambda u[1,0]+(1-2\lambda)u[2,0]+\lambda u[3,0] u[2,1] =0.025\cos \Big( \frac{3π}{2}0.2\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.4\Big) +0.025 \cos \Big( \frac{3π}{2}0.6\Big)

iteración i=3, j=0

u[3,1] =\lambda u[2,0]+(1-2\lambda)u[3,0]+\lambda u[4,0] u[3,1] =0.025\cos \Big( \frac{3π}{2}0.4\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.6\Big) +0.025 \cos \Big( \frac{3π}{2}0.8\Big)

iteración i=4, j=0

u[4,1] =\lambda u[3,0]+(1-2\lambda)u[4,0]+\lambda u[5,0] u[4,1] =0.025\cos \Big( \frac{3π}{2}0.6\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.8\Big) +0.025 (0)

continuar con las iteraciones en el algoritmo

Resultados con el algoritmo

Tabla de resultados
[[ 1.    1.    1.    1.    1.    1.    1.    1.    1.    1.  ]
 [ 0.59  0.58  0.56  0.55  0.54  0.53  0.53  0.52  0.51  0.5 ]
 [-0.31 -0.3  -0.3  -0.29 -0.28 -0.28 -0.27 -0.27 -0.26 -0.26]
 [-0.95 -0.93 -0.91 -0.89 -0.88 -0.86 -0.84 -0.82 -0.81 -0.79]
 [-0.81 -0.79 -0.78 -0.76 -0.74 -0.73 -0.71 -0.7  -0.68 -0.67]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 1
Tb = 0
#T0 = 25
fx = lambda x: np.cos(3*np.pi/2*x)
# longitud en x
a = 0
b = 1
# Constante K
K = 2
# Tamaño de paso
dx = 0.2
dt = dx/100
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
    
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

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