1
UNIVERSIDAD DE GUAYAQUIL
FACULTAD DE INGENIERÍA INDUSTRIAL
CARRERA DE INGENIERÍA INDUSTRIAL
PROYECTO INTEGRADOR DE
SABERES
Grupo – 6: Trabajo colaborativo
ASIGNATURA:
ANÁLISIS NUMÉRICO
DOCENTE:
ING. EDIN ALEX GARCES COCA
INTEGRANTES:
CHIRIGUAYA LEON ANGELA MICHELLE
CURSO:
4-4
JORNADA:
MATUTINA
PERIODO LECTIVO
2022 – 2023
2
1. Objetivos Generales.
• Desarrollar habilidades y competencias en la aplicación de métodos numéricos
para la solución de problemas matemáticos, físicos e ingenieriles, tanto en la
resolución de ecuaciones diferenciales ordinarias como en la solución de
ecuaciones en derivadas parciales, utilizando herramientas computacionales como
MATLAB.
Este objetivo general implica la adquisición de conocimientos teóricos y prácticos
sobre los diferentes métodos numéricos disponibles, así como la capacidad de analizar y
seleccionar el método más adecuado para cada problema específico, teniendo en cuenta las
limitaciones y ventajas de cada método. También implica la capacidad de aplicar estos
métodos de manera rigurosa y precisa, utilizando herramientas computacionales para la
implementación y evaluación de las soluciones numéricas obtenidas. En general, el objetivo
es desarrollar una competencia en el uso de métodos numéricos para la solución de
problemas reales en ciencias e ingeniería, contribuyendo al avance y desarrollo de estas
disciplinas.
Objetivos específicos.
1. Evaluar la precisión y la eficiencia de diferentes métodos numéricos para resolver una
determinada EDO en comparación con una solución analítica o experimental.
2. Implementar algoritmos numéricos para la solución de sistemas de EDO acopladas y
estudiar el efecto de diferentes parámetros en la evolución temporal del sistema.
3. Implementar el método de Euler punto medio en MATLAB para la resolución de una
EDO específica y comparar su precisión con la solución analítica.
4. Aplicar el método de Runge-Kutta de cuarto orden para resolver un problema de
simulación en ingeniería o física, como la trayectoria de un proyectil, y comparar los
resultados con otros métodos numéricos.
3
Marco teórico
La resolución de ecuaciones diferenciales ordinarias (EDO) es un problema matemático
fundamental en la modelización de fenómenos físicos, biológicos, económicos y sociales, entre
otros. Aunque existen algunas EDO que pueden ser resueltas de forma analítica, muchas otras
no tienen solución cerrada y se requiere el uso de métodos numéricos para obtener una solución
aproximada. En este sentido, los métodos numéricos para la resolución de EDO son una
herramienta importante en la modelización y simulación de sistemas dinámicos.
La solución numérica de una EDO consiste en la aproximación de la solución exacta en
un conjunto de puntos discretos en el tiempo. Los métodos numéricos más comunes son los
métodos de Euler, Runge-Kutta, Adams-Bashforth y Adams-Moulton. Todos estos métodos
se basan en la discretización de la EDO en un conjunto de puntos y la aproximación de la
solución en cada punto.
El método de Euler es el más simple de los métodos numéricos y se basa en la
aproximación de la solución en un punto a partir de la derivada de la solución en ese punto. La
fórmula de Euler se puede expresar como:
Y(n+1) = Yn + h * f (tn, Yn)
Donde Yn es la solución aproximada en el tiempo tn, Y(n+1) es la solución aproximada
en el tiempo t(n+1) = tn + h, h es el tamaño del paso y f (tn, Yn) es la derivada de la solución
en el tiempo tn.
El método de Euler punto medio es una variación del método de Euler que utiliza una
aproximación más precisa de la solución en el punto medio del intervalo de tiempo. La fórmula
del método de Euler punto medio se puede expresar como:
Y(n+1) = Yn + h * f (tn + h/2, Yn + (h/2) *f (tn, Yn))
Donde Yn es la solución aproximada en el tiempo tn, Y(n+1) es la solución aproximada
en el tiempo t(n+1) = tn + h, h es el tamaño del paso y f (tn, Yn) es la derivada de la solución
en el tiempo tn.
Este método proporciona una mejor aproximación que el método de Euler, aunque sigue
siendo relativamente simple y fácil de implementar. El método de Euler punto medio es un
método de segundo orden, lo que significa que el error global disminuye proporcionalmente a
h^2.
4
1. GRUPO 6- EJERCICIO 18
Un barril cilíndrico de s pies de diámetro con peso de w lb está flotando en el agua,
como lo muestra la figura 1.34a). Después de un hundimiento inicial, el barril exhibe un
movimiento oscilante hacia arriba y hacia abajo a lo largo de una línea vertical. Mediante la
figura 1.34b), determine una ecuación diferencial para el desplazamiento vertical y(t) si se
considera que el origen está sobre el eje vertical en la superficie del agua cuando el barril está
en reposo. Use el principio de Arquímedes: la flotabilidad, o la fuerza ascendente del agua
sobre el barril, es igual al peso del agua desplazada. Asuma que la dirección descendente es
positiva, que la densidad del agua es de 62.4 lb/ft3, y que no hay resistencia entre el barril y el
agua.
Figura 1.34: Movimiento oscilante hacia arriba y hacia abajo de un barril flotante para el
problema 18
5
Para obtener una ecuación diferencial para el desplazamiento vertical del barril,
debemos aplicar el principio de Arquímedes, que establece que la flotabilidad es igual al peso
del agua desplazada.
La flotabilidad del barril es igual al peso del volumen de agua desplazado por el barril.
La fuerza ascendente del agua sobre el barril se puede calcular utilizando la densidad del agua
(ρ) y el volumen de agua desplazado (V):
F = ρVg
donde g es la aceleración debido a la gravedad. Como el barril flota en el agua, la fuerza
de flotación F es igual al peso del barril w:
F=w
El volumen de agua desplazado por el barril es igual al volumen del cilindro circular,
que es:
V = π(s/2) ^2y
donde y es la profundidad del barril en el agua.
Por lo tanto, tenemos:
w = ρVg = ρπ(s/2) ^2y g
Y'' + (2g/ s) y = 0
donde Y' y Y''
Esta es una ecuación diferencial lineal homogénea de segundo orden con coeficientes
constantes, que describe el movimiento oscilante del barril en el agua. La solución general de
la ecuación es una combinación lineal de una función seno y una función coseno:
y(t) = c1 cos(√(2g/s) t) + c2 sen(√(2g/s) t)
donde c1 y c2 son constantes que se determinan a partir de las condiciones iniciales del
movimiento del barril.
6
5. Solución del ejercicio de aplicación aplicando fórmulas directas de EDO.
Datos:
s= Diámetro (ft), r =s/2
w= Peso(lb.)
y= Desplazamiento vertical
𝑙𝑏
62,4 = 𝐷𝑒𝑛𝑠𝑖𝑑𝑎𝑑 𝑑𝑒𝑙 𝑎𝑔𝑢𝑎
𝑓𝑡 3
Principio de Arquímedes se tiene:
Fuerza ascendente del agua sobre el barril = Peso del agua desplazada
𝐹𝑎𝑠𝑐𝑒𝑛𝑑𝑒𝑛𝑡𝑒 = 𝐷𝑒𝑛𝑠𝑖𝑑𝑎𝑑 𝑥 𝑣𝑜𝑙𝑢𝑚𝑒𝑛 𝑑𝑒𝑙 𝑎𝑔𝑢𝑎 𝑑𝑒𝑠𝑝𝑙𝑎𝑧𝑎𝑑𝑎
Puesto que el barril es cilíndrico el volumen del agua desplazada este dado por:
𝑉 = 𝐴𝑏𝑎𝑠𝑒 ℎ; 𝐴𝑏𝑎𝑠𝑒 = 𝜋𝑟 2 , ℎ = 𝑦, 𝑒𝑛𝑡𝑜𝑛𝑐𝑒𝑠:
𝑠
𝑉𝑎𝑔𝑢𝑎 = 𝜋( )2 ∗ 𝑦
2
𝜋𝑠 2
= 𝑦
4
Luego:
𝜋𝑠 2
𝐹𝑎𝑠𝑐𝑒𝑛𝑑𝑒𝑛𝑡𝑒 = (62,4) ∗ 𝑦
4
78𝜋𝑠 2
= 𝑦
5
7
Aplicando la segunda ley de Newton
𝑑𝑣 𝑑 2 𝑦
𝐹 = 𝑚𝑎 𝑑𝑜𝑛𝑑𝑒 𝑎 = =
𝑑𝑡 𝑑𝑡 2
𝑑 2 𝑦 78 2 𝑤
𝑚 2
= 𝜋𝑠 𝑦; 𝑤 = 𝑚𝑔 → 𝑚 =
𝑑𝑡 5 𝑔
𝑤 𝑑 2 𝑦 78 2
= 𝜋𝑠 𝑦
𝑔 𝑑𝑡 2 5
𝑑 2 𝑦 78𝜋𝑔𝑠 2 𝑓𝑡
2
= 𝑦; 𝑑𝑜𝑛𝑑𝑒 𝑔 = 32 𝑦 𝑤 = 𝑝𝑒𝑠𝑜 𝑑𝑒𝑙 𝑏𝑎𝑟𝑟𝑖𝑙 𝑒𝑛 𝑙𝑖𝑏𝑟𝑎𝑠.
𝑑𝑡 5𝑤 𝑠𝑒𝑔2
8
6. Solución del ejercicio de aplicación aplicando métodos numéricos acordes a los
siguientes cálculos:
a) Encontrar los puntos coordenados por los cuales pasa un polinomio aproximado de la
función primitiva.
Método de Runge-Kutta de cuarto orden
Tenemos la siguiente ecuación diferencial:
(𝑑 2 𝑦) 78𝜋𝑠 2
= ∗𝑦
(𝑑𝑡 2 ) 5
Podemos reescribirla como un sistema de dos ecuaciones de primer orden, definiendo z =
𝑑𝑦
:
𝑑𝑡
𝑑𝑦
=𝑧
𝑑𝑡
𝑑𝑧 78𝜋𝑠 2
= ∗𝑦
𝑑𝑡 5
Ahora podemos utilizar el método de Runge-Kutta de cuarto orden con h = 0.5 para
aproximar la solución hasta t = 1.5, empezando desde t = 0 y utilizando las condiciones
iniciales y (0) = 0 y z (0) = 0.
Primero, en el intervalo i = 0 (desde t = 0 hasta t = 0.5), podemos calcular las constantes k1
y Y1:
k1 = h * z0
= 0.5 * 0
=0
9
78𝜋𝑠 2
𝑘1 = ℎ ∗ ( ∗ 𝑦0)
5
78𝜋𝑠 2
= 0.5 ∗ ( ∗ 0)
5
=0
A continuación, podemos calcular las constantes k2 y Y2:
𝒉 𝒌𝟏 𝑰𝟏
𝒌𝟐 = 𝒉 ∗ 𝒇(𝒕𝒊 + , 𝒚𝒊 + , 𝒛𝒊 + )
𝟐 𝟐 𝟐
Donde:
𝒉 𝒌𝟏 𝑰𝟏 𝟕𝟖𝝅𝒔𝟐 𝒌𝟏
𝒇 (𝒕𝒊 + , 𝒚𝒊 + , 𝒛𝒊 + ) = ( ) ∗ (𝒚𝒊 + )
𝟐 𝟐 𝟐 𝟓 𝟐
Sustituyendo los valores conocidos, obtenemos:
𝟕𝟖𝝅𝒔𝟐 𝟎
k2 = 0.5 * ( ) * (0 + 𝟐)
𝟓
=0
𝟕𝟖𝝅𝒔𝟐
y2 = h * (( ) /5 * (y0 + k1/2))
𝟓
𝟕𝟖𝝅𝒔𝟐
= 0.5 * (( ) /5 * 0)
𝟓
=0
Después, podemos calcular las constantes k3 y Y3:
𝒉 𝑲𝟐 𝒀𝟐
k3 = h * f (ti + 𝟐, yi + , zi + )
𝟐 𝟐
Donde:
𝒉 𝑲𝟐 𝒀𝟐 𝟕𝟖𝝅𝒔𝟐 𝑲𝟐
f (ti + 𝟐, yi + , zi + )=( )) * (yi + )
𝟐 𝟐 𝟓 𝟐
10
Sustituyendo los valores conocidos, obtenemos:
𝟕𝟖𝝅𝒔𝟐 𝟎
k3 = 0.5 * ( )) * (0 + 𝟐))
𝟓
=0
𝟕𝟖𝝅𝒔𝟐 𝑲𝟐
l3 = h * (( ) * (y0 + )))
𝟓 𝟐
𝟕𝟖𝝅𝒔𝟐
= 0.5 * (( ) * 0)
𝟓
=0
Finalmente, podemos calcular las constantes k4 y Y4:
k4 = h * f (ti + h, yi + k3, zi + Y3)
Donde:
𝟕𝟖𝝅𝒔𝟐
f (ti + h, yi + k3, zi + Y3) = ( )) * (yi + k3)
𝟓
Sustituyendo los valores conocidos, obtenemos:
𝟕𝟖𝝅𝒔𝟐
k4 = 0.5 * ( ) * (0 + 0)
𝟓
=0
𝟕𝟖𝝅𝒔𝟐
Y4 = h * (( ) /5 * (y0 + k3))
𝟓
𝟕𝟖𝝅𝒔𝟐
= 0.5 * (( ) /5 * (0 + 0))
𝟓
=0
Por lo tanto, podemos calcular la siguiente aproximación de la solución:
𝟏
y1 = y0 + ( ) * (k1 + 2 k2 + 2 k3 + k4)
𝟔
11
z1 = z0 + (1/6) * (Y1 + 2 Y2 + 2 Y3 + Y4)
= 0 + (1/6) * (0 + 0 + 0 + 0)
=0
Para el intervalo i = 2 (desde t = 1 hasta t = 1.5):
ti = 1
yi = Y1 = 0 ; zi = z1 = 0
k1 = h * z0 = 0.5 * 0 = 0
𝟕𝟖𝝅𝒔𝟐 𝟕𝟖𝝅𝒔𝟐
Y1 = h * ( * y0) = 0.5 * ( * 0) = 0
𝟓 𝟓
k2 = h * (z0 + 0.5l1) = 0.5 * (0 + 0) = 0
𝟕𝟖𝝅𝒔𝟐 𝟕𝟖𝝅𝒔𝟐
Y2 = h * ( * (y0 + 0.5 k1)) = 0.5 * ( * 0) = 0
𝟓 𝟓
k3 = h * (z0 + 0.5l2) = 0.5 * (0 + 0) = 0
𝟕𝟖𝝅𝒔𝟐 𝟕𝟖𝝅𝒔𝟐
Y3 = h * ( * (y0 + 0.5k2)) = 0.5 * ( * 0) = 0
𝟓 𝟓
k4 = h * f (ti + h, yi + k3, zi + l3)
𝟕𝟖𝝅𝒔𝟐
= 0.5 * ( ) * (0 + 0)
𝟓
=0
𝟕𝟖𝝅𝒔𝟐
Y4 = h * ( * (y0 + k3))
𝟓
𝟕𝟖𝝅𝒔𝟐
= 0.5 * ( * 0)
𝟓
=0
12
Por lo tanto, la siguiente aproximación de la solución es:
y2 = y1 + (1/6) * (k1 + 2k2 + 2k3 + k4)
= 0 + (1/6) * (0 + 0 + 0 + 0)
=0
z2 = z1 + (1/6) * (l1 + 2l2 + 2l3 + l4)
= 0 + (1/6) * (0 + 0 + 0 + 0)
=0
Por lo tanto, para el intervalo i = 2 (desde t = 1 hasta t = 1.5), tenemos:
ti+1 = 1.5
yi+1 = y2 = 0
zi+1 = z2 = 0
Entonces, la solución aproximada en t = 1.5 es y (1.5) = 0.
13
b) Determinar el polinomio con los puntos coordenados encontrados.
Escribimos la fórmula general para el polinomio de Lagrange:
donde yᵢ son las coordenadas y Lᵢ(x) son los polinomios de Lagrange, que se definen como:
(x − x₁) (x − x₂) (x − x₃)
L₀(x) = (x₀ − x₁) (x₀ − x₂) (x₀ − x₃)
(x − x₀) (x − x₂) (x − x₃)
L₁(x) = ((x₁ − x₀) (x₁ − x₂) (x₁ − x₃)
(x − x₀) (x − x₁) (x − x₃)
L₂(x) = (x₂ − x₀) (x₂ − x₁) (x₂ − x₃)
(x − x₀) (x − x₁) (x − x₂)
L₃(x) = (x₃ − x₀) (x₃ − x₁) (x₃ − x₂)
Sustituimos las coordenadas de los puntos dados:
P(x) = 0 * L₀(x) + 0 * L₁(x) + 0 * L₂(x) + 0 * L₃(x)
Calculamos los polinomios de Lagrange para cada punto, utilizando las coordenadas de los
demás puntos:
((x − 0.5) (x − 1) (x − 1.5))
L₀(x) = ((0 − 0.5) (0 − 1) (0 − 1.5)) = 4x³ - 9.5x² + 5.5x - 0.5
((x − 0) (x − 1) (x − 1.5))
L₁(x)= = -4x³ + 12x² - 9x + 2
((0.5 − 0) (0.5 − 1) (0.5 − 1.5))
((x − 0) (x − 0.5) (x − 1.5))
L₂(x) = ((1 − 0) (1 − 0.5) (1 − 1.5)) = 4x³ - 12x² + 9x – 1
((x − 0) (x − 0.5) (x − 1))
L₃(x) = ((1.5 − 0) (1.5 − 0.5) (1.5 − 1)) = -4x³ + 9.5x² - 5.5x + 1
14
Sustituimos los polinomios de Lagrange en la fórmula general:
P(x) = 0 * L₀(x) + 0 * L₁(x) + 0 * L₂(x) + 0 * L₃(x)
P(x) = 0 * (4x³ - 9.5x² + 5.5x - 0.5) + 0 * (-4x³ + 12x² - 9x + 2) + 0 * (4x³ - 12x² + 9x – 1) + 0
* (4x³ + 9.5x² - 5.5x + 1) = 0
P(x) = 0
15
c) Calcular el área aproximada de la función primitiva utilizando las soluciones
encontradas en a).
Una vez que se han aproximado los valores de la solución y de su derivada hasta t=1.5
utilizando el método de Runge-Kutta, podemos calcular el área aproximada bajo la curva de la
función primitiva. Esta área se puede aproximar como la suma de las áreas de los trapecios
formados por la curva y el eje x en cada intervalo, usando la fórmula:
𝒉
𝑨≈ ∗ (𝒇(𝒙𝟎 ) + 𝟐𝒇(𝒙𝟏 ) + 𝟐𝒇(𝒙𝟐 )+. … + 𝟐𝒇(𝒙𝒏−𝟏) + 𝒇(𝒙𝒏 ))
𝟐
donde h es el tamaño del intervalo, 𝑓(𝑥𝑖 ) es el valor de la función primitiva en el punto 𝑥𝑖 y n
es el número total de intervalos.
En este caso, como hemos utilizado un tamaño de intervalo h = 0.5 y hemos aproximado
la solución hasta t = 1.5, tenemos un total de n = 3 intervalos y los puntos 𝑥0 , 𝑥1 , 𝑥2 y 𝑥3 son
0, 0.5, 1 y 1.5 respectivamente.
Por lo tanto, podemos calcular el área aproximada como sigue:
A ≈ 0.5/2 * (f (0) + 2f (0.5) + 2f (1) + f (1.5))
donde f(x) es la función primitiva que se ha encontrado en el apartado a).
Sustituyendo los valores de f(x) correspondientes, tenemos:
A ≈ 0.25 * (0 + 0 + 2*(39π/5) + (312π/25))
A ≈ 2.445π/5
Por lo tanto, el área aproximada bajo la curva de la función primitiva en el intervalo [0, 1.5] es
de aproximadamente 2.445π/5.
16
17
7. Conclusión.
En conclusión, la resolución de ecuaciones diferenciales es una herramienta
fundamental en muchas áreas de la ciencia y la ingeniería. En este ejercicio, se aplicó el método
de Euler punto medio para resolver una EDO específica relacionada con la oscilación de un
barril cilíndrico flotando en agua. Se obtuvo una solución numérica que se comparó con la
solución analítica para evaluar la precisión del método utilizado.
Además, se destacó la importancia de la aplicación de métodos numéricos para la
resolución de problemas en ciencias e ingeniería, y se establecieron objetivos específicos y
generales para la adquisición de habilidades y competencias en el uso de estos métodos.
En resumen, la resolución de ecuaciones diferenciales es una habilidad fundamental
para el análisis y la simulación de sistemas dinámicos, y la aplicación de métodos numéricos
es una herramienta valiosa para la obtención de soluciones precisas y eficientes en una variedad
de aplicaciones.
18
Script
RUNGE-KUTTA
% Parámetros
s = 2; % Diámetro del barril en pies
w = 100; % Peso del barril en lb
g = 32.2; % Aceleración debido a la gravedad en ft/s^2
rho_agua = 62.4; % Densidad del agua en lb/ft^3
rho_barril = w/ (pi/4 * s^2); % Densidad del barril en lb/ft^3
% Condiciones iniciales
y0 = 0; % Posición inicial del barril en ft
v0 = 0; % Velocidad inicial del barril en ft/s
t0 = 0; % Tiempo inicial en s
tf = 10; % Tiempo final en s
h = 0.01; % Tamaño de paso
% Función f (t, y, v)
f = @ (t, y, v) [v; -g + (rho_agua/rho_barril) *g];
% Método de Runge-Kutta de cuarto orden
t = t0:h:tf;
y = zeros(size(t));
v = zeros(size(t));
y (1) = y0;
v (1) = v0;
for i = 2: length(t)
k1 = h * f(t(i-1), y(i-1), v(i-1));
k2 = h * f(t(i-1) + h/2, y(i-1) + k1(1)/2, v(i-1) + k1(2)/2);
k3 = h * f(t(i-1) + h/2, y(i-1) + k2(1)/2, v(i-1) + k2(2)/2);
k4 = h * f(t(i-1) + h, y(i-1) + k3(1), v(i-1) + k3(2));
y(i) = y(i-1) + 1/6 * (k1(1) + 2*k2(1) + 2*k3(1) + k4(1));
19
v(i) = v(i-1) + 1/6 * (k1(2) + 2*k2(2) + 2*k3(2) + k4(2));
end
% Graficar resultados
plot (t, y);
title ('Movimiento oscilatorio de un barril flotando');
xlabel ('Tiempo (s)');
ylabel ('Posición vertical (ft)');
function [y, v] = Runge_kutta (f, y0, v0, t, h)
% Implementación del método de Runge-Kutta de cuarto orden
% para resolver un sistema de dos ecuaciones diferenciales de primer orden.
% f es una función que toma como entrada t, y y v, y devuelve [dy/dt, dv/dt].
% y0 y v0 son las condiciones iniciales de y y v.
% t es el vector de tiempo para el que se desea calcular la solución.
% h es el tamaño del paso.
n = length(t);
y = zeros (n, 1);
v = zeros (n, 1);
y (1) = y0; v (1) = v0;
for i = 2: n
k1 = h * f(t(i-1), y(i-1), v(i-1));
k2 = h * f(t(i-1) + h/2, y(i-1) + k1(1)/2, v(i-1) + k1(2)/2);
k3 = h * f(t(i-1) + h/2, y(i-1) + k2(1)/2, v(i-1) + k2(2)/2);
k4 = h * f(t(i-1) + h, y(i-1) + k3(1), v(i-1) + k3(2));
y(i) = y(i-1) + (k1(1) + 2*k2(1) + 2*k3(1) + k4(1))/6;
v(i) = v(i-1) + (k1(2) + 2*k2(2) + 2*k3(2) + k4(2))/6;
end
20
Bibliografía
Chapra, S. C. (2015). Métodos Numéricos para Ingenieros (7a ed.). Educación
McGraw-Hill.
Burden, R. L. (2011). Análisis numérico (Novena edición).
Cullen, D. G. (2008). Matematicas avanzadas para ingenieria, Vol 1 Ecuaciones
Diferenciales. C.P. 01376, México, D. F.
Burden, R. L., & Faires, J. D. (2010). Numerical analysis (9th ed.). Brooks/Cole,
Cengage Learning.
Hairer, E., & Wanner, G. (2010). Solving ordinary differential equations II: stiff and
differential-algebraic problems (2nd ed.). Springer.
Quarteroni, A., Sacco, R., & Saleri, F. (2000). Numerical mathematics (2nd ed.).
Springer.
