Análisis Numérico para Ingeniería

Clase Nro. 4

“Ver lo que está delante

de nuestros ojos requiere
un esfuerzo constante.”

Eric Arthur Blair

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 1

Temas a tratar en esta clase :

Sistemas de EDOs de 1er. Orden.

Transformación de EDOs de orden n.

Resolución de Sistemas de 1er. Orden.

Estrategias de ajuste del paso h.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 2

 Sistemas de EDOs de 1er. orden

Infinidad de problemas de diversas áreas del conocimiento

como la física, la química, la biología, la economía, la
ingeniería, etc. se expresan a menudo en términos de
ecuaciones diferenciales ordinarias.
Estas ecuaciones pueden estar relacionadas entre sí, por lo
que conforman sistemas de ecuaciones diferenciales.
Las ecuaciones diferenciales de orden n, pueden
transformarse en sistemas de ecuaciones de primer orden,
para su resolución.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 4

 Ejemplo de Sistemas de EDOs
El modelo SIR, desarrollado por el químico escocés William
Ogilvy Kernack y el matemático A.G. McKendrick , modela
la evolución de una enfermedad infecto-contagiosa.

El modelo divide MODELO SIR

a la población en SUCEPTIBLE
individuos dS
suceptibles de =−β⋅S (t )⋅I (t)
dt
ser infectados, INFECTADO dI
individuos =β⋅S(t)⋅I (t )− γ⋅I (t )
infectados e dt
individuos dR
= γ⋅I (t )
recuperados de RECUPERADO dt
la enfermedad.
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 5
 Ejemplo de Sistemas de EDOs
En el año 1978 en la revista científica British Medical Journal
se publicaron los datos de un brote de gripe en un internado
del norte de Inglaterra que se extendio del 22 de enero al 4
de febrero, infectando a 512 de las 763 personas que allí
estudiaban.

Debido a que en este caso concreto se cuenta con datos

reales y la cantidad de individuos permanece constante, es
un escenario ideal para corroborar la exactitud del modelo
SIR para predecir la evolución de la infección.

Precisamente para este caso en particular se han

determinado los siguientes valores: γ = 0.4477 y β = 0.0022.
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 6
 Ejemplo de Sistemas de EDOs
En este caso los valores iniciales son:
Individuos Suceptibles: S0 = 762,
Individuos Infectados: I0 = 1,
Individuos Recuperados: R0 = 0

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 7

Modelo SEIR
El modelo SEIR, es una adaptación del modelo SIR explicado
anteriormente

El modelo divide
SUCEPTIBLE MODELO SEIR
a la población en
individuos dS
=−β⋅S (t )⋅I (t)/ N
suceptibles de EXPUESTO dt
ser infectados, dE
individuos =β⋅S (t)⋅I (t )/ N − σ⋅E(t)
dt
expuestos, dI
individuos INFECTADO dt
= σ⋅E(t )− γ⋅I (t )
infectados e
individuos dR
= γ⋅I (t )
recuperados de dt
la enfermedad. RECUPERADO
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 8
 Modelo SEIR para el COVID-19

Para los que

deseen saber
más sobre el
tema y quieran
experimentar
con estos
modelos, les
dejo el enlace
a este Blog.

BLOG DEL INSTITUTO DE MATEMÁTICAS DE LA UNIVERSIDAD DE SEVILLA

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 9
 EDOs de orden n
Una EDO de orden n puede transformarse en un sistema de n
ecuaciones de primer orden.
n 2 n−1
d y dy d y d y
=f  x , y , , 2, , n−1 
dx n
dx dx dx
y∣x= x = y 0
0

dy
∣
dx x= x
= y '0
0

d(2) y CONDICIONES
2
dx x = x
………………
∣
= y ' '0
0
INICIALES

d(n−1) y
dx n−1 x= x
= y 0
∣
(n−1)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 10

 Algoritmo de transformación
Las EDOs de orden n se transforman en sistemas de n ecuaciones
de primer orden para poder resolverlas con los métodos ya vistos.

dn y dy d 2 y d n−1 y
=f  x , y , , 2, , n−1 =f  x , y , p , q , , w 
dx n
dx dx dx

dy
= y '= p
dx
2
d y Sustituimos cada
2
= y ' '= p '= q derivada por una
dx

nueva variable.
d n−1 y
n−1
=z '= w
dx

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 11

 Ejemplo (I)
Ecuación diferencial de orden 3
y ' ' 'x⋅y ' '−x⋅y '−2⋅y=x
x 0=0 ; y x = y ' ' x =0 ; y ' x =1
0 0 0

Sistema de 3 ecuaciones diferenciales de 1er. orden

y '= p Sustituimos cada

derivada por una
p '= q nueva variable.
q '= x− x⋅q x⋅p2⋅y

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 12

 Resolución completa de un problema
Sistema Masa – Resorte - Amortiguador

Sistema de 2 ecuaciones diferenciales de orden 2

M 1⋅y ' ' 1 +B 1⋅y ' 1 +k 1⋅y 1−B 1⋅y ' 2−k 2⋅y 2=F 1 (t)
−B 1⋅y ' 1−k 1⋅y 1 + M 2⋅y ' ' 2 +B 1⋅y ' 2 +(k 1 +k 2 )⋅y 2=0
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 16
 Transformación del Ejemplo (II)
Sistema de 2 ecuaciones diferenciales de orden 2
M 1⋅y ' ' 1B 1⋅y ' 1k 1⋅y 1−B 1⋅y ' 2−k 2⋅y 2=F 1 t 
−B 1⋅y ' 1−k 1⋅y 1M 2⋅y ' ' 2B 1⋅y ' 2k 1k 2 ⋅y 2=0

Sistema de 4 ecuaciones diferenciales de 1er. orden

y 1 '= v 1
F1  t − B1⋅v 1  B1⋅v 2 − k 1⋅y 1  k 2⋅y 2
v 1 '=
M1
y 2 '= v 2
B1⋅v 1  k 1⋅y 1− B1⋅v 2 − k1  k 2 ⋅y 2
v 2 '=
M2

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 17

 Implementación con vectores

La utilización de vectores simplifica notablemente la

implementación de los métodos para resolver sistemas
de ecuaciones diferenciales.

v= t y1 v1 y2 v2 Vector de variables
v(0) v(1) v(2) v(3) v(4)

vp = 1 y'1 v'1 y'2 v'2 Vector de derivadas

vp(0) vp(1) vp(2) vp(3) vp(4)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 18

 Implementación con vectores (II)
y 1 '= v 1
F1  t − B1⋅v 1  B1⋅v 2 − k 1⋅y 1  k 2⋅y 2
v 1 '=
M1
y 2 '= v 2
B1⋅v 1  k 1⋅y 1− B1⋅v 2 − k 1  k 2 ⋅y 2
v 2 '=
M2

v= t y1 v1 y2 v2 Vector de variables
v(0) v(1) v(2) v(3) v(4)

vp = 1 y'1 v'1 y'2 v'2 Vector de derivadas

vp(0) vp(1) vp(2) vp(3) vp(4)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 19

 Implementación con vectores (III)
vp(1)= v (2)
F1 ( v ( 0))− B1⋅v ( 2)+ B1⋅v ( 4)− k 1⋅v ( 1)+ k 2⋅v ( 3)
vp(2)=
M1
vp(3)= v (4)
B1⋅v ( 2)+ k 1⋅v (1)− B1⋅v ( 4)−( k 1 + k 2 )⋅v ( 3)
vp(4)=
M2

v= t y1 v1 y2 v2
v(0) v(1) v(2) v(3) v(4)
Vector de variables

vp = 1 y'1 v'1 y'2 v'2

vp(0) vp(1) vp(2) vp(3) vp(4) Vector de derivadas

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 20

 Definición de la Función Derivada
MODULE sistema_MRA
! Cantidad de Ecuaciones del Sistema
INTEGER, PARAMETER :: cant_ec = 4

CONTAINS

FUNCTION v_prima(v)
! Definición del Sistema de Ecuaciones Diferenciales
REAL(8), PARAMETER :: m1=1.0, m2=0.05, k1=5.0, k2=5.0, b1=1.0
REAL(8), DIMENSION(0:cant_ec) :: v, v_prima

v_prima(0) = 1.0
v_prima(1) = v(2)
v_prima(2) = (-b1*v(2)+b1*v(4)-k1*v(1)+k2*v(3))/m1
v_prima(3) = v(4)
v_prima(4) = (-b1*v(4)+b1*v(2)-k1*v(3)+k1*v(1)-k2*v(3))/m2

END FUNCTION v_prima

END MODULE
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 21
 Euler Simple (Implementación)
SUBROUTINE EulerSimple(vi, cant_ec, max_iter, h, formato)
!Metodo de Euler Simple
INTEGER, INTENT(IN) :: cant_ec, max_iter
REAL(8), INTENT(IN), DIMENSION(0:cant_ec) :: vi
REAL(8), INTENT(IN) :: h
CHAR*15, INTENT(IN) :: formato
INTEGER iter
REAL(8), DIMENSION(0:cant_ec) :: v

WRITE (*, formato) 0, vi

v = vi
DO iter = 1, max_iter
v = v + h*v_prima(v)
WRITE (*, formato) iter, v
END DO

END SUBROUTINE EulerSimple

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 22

 Solución Euler Simple (h=0.1)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 23

 Euler Modificado (Implementación)
SUBROUTINE EulerModificado(vi, cant_ec, max_iter, h, formato)
!Metodo de Euler Modificado
INTEGER, INTENT(IN) :: cant_ec, max_iter
REAL(8), INTENT(IN), DIMENSION(0:cant_ec) :: vi
REAL(8), INTENT(IN) :: h
CHAR*15, INTENT(IN) :: formato
INTEGER iter
REAL(8), DIMENSION(0:cant_ec) :: v
WRITE (*, formato) 0, vi
v = vi
DO iter = 1, max_iter
vp = v_prima(v)
v = v + h*(vp + v_prima(v + h*vp))/2.0
WRITE (*, '(I5, formato) iter, v
END DO

END SUBROUTINE EulerModificado

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 24

 Solución Euler Modificado (h=0.1)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 25

 Runge-Kutta (4to. Orden)
SUBROUTINE RungeKutta(vi, cant_ec, max_iter, h, formato)
!Metodo de Runge-Kutta (de 4to. orden)
REAL(8), INTENT(IN), DIMENSION(0:cant_ec) :: vi
REAL(8), INTENT(IN) :: h
INTEGER, INTENT(IN) :: cant_ec, max_iter
CHAR*15, INTENT(IN) :: formato
INTEGER iter
REAL(8), DIMENSION(0:cant_ec) :: v, k1, k2, k3, k4
WRITE (*, formato) 0, vi
v = vi
DO iter = 1, max_iter
k1 = h*v_prima(v)
k2 = h*v_prima(v+k1/2.0)
k3 = h*v_prima(v+k2/2.0)
k4 = h*v_prima(v+k3)
v = v + (k1 + 2.0*k2 + 2.0*k3 +k4)/6.0
WRITE (*, formato) iter, v
END DO
END SUBROUTINE RungeKutta
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 26
 Solución Runge-Kutta 4to. Orden (h=0.1)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 27

 Runge-Kutta-Fehlberg
SUBROUTINE RungeKuttaFehlberg(vi, cant_ec, max_iter, h, formato)
!Metodo de Runge-Kutta-Fehlberg (de 6to. orden)
REAL(8), INTENT(IN), DIMENSION(0:cant_ec) :: vi
REAL(8), INTENT(IN) :: h
INTEGER, INTENT(IN) :: cant_ec, max_iter
CHAR*15, INTENT(IN) :: formato
INTEGER iter
REAL(8), DIMENSION(0:cant_ec) :: v, k1, k2, k3, k4, k5, k6, e

WRITE (*, formato) 0, vi, 0

v = vi
DO iter = 1, max_iter
k1 = h*v_prima(v)
k2 = h*v_prima(v + k1/4.0)
k3 = h*v_prima(v + (3.0*k1 + 9.0*k2)/32.0)
k4 = h*v_prima(v + (1932.0*k1 - 7200.0*k2 + 7296.0*k3)/2197.0)
k5 = h*v_prima(v + 439.0*k1/216.0 - 8.0*k2 + 3680.0*k3/513.0 - 845.0*k4/4104.0)
k6 = h*v_prima(v - 8.0*k1/27.0 + 2.0*k2 - 3544.0*k3/2565.0 + 1859.0*k4/4104.0 - 11.0*k5/40.0)
v = v + (25.0*k1/216.0 + 1408.0*k3/2565.0 + 2197.0*k4/4104.0 - k5/5.0)
e = k1/360.0 - 128.0*k3/4275.0 - 2197.0*k4/75240.0 + k5/50.0 + 2.0*k6/55.0
WRITE (*, formato) iter, v, e
END DO
END SUBROUTINE RungeKuttaFehlberg

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 28

 Runge-Kutta-Fehlberg (h=0.1)

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 29

 Estrategias de ajuste del paso h

La exactitud de los métodos numéricos de

resolución de ecuaciones diferenciales
depende fuertemente del paso h elegido.

Existen diferentes estrategias para ajustar

automáticamente el paso h y de esta forma
mantener el error de la solución por debajo de
una cierta tolerancia establecida.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 45

Estimación del Error
En aquellos métodos que no tienen una fórmula específica para la
estimación del error, es posible estimarlo por medio de la diferencia
del valor calculados con un paso de avance h y del obtenido al
avanzar dos pasos h/2.

Valor calculado en un
solo paso h.
yn+1
ERROR
ESTIMADO
y*n+1 Valor calculado en dos
pasos h/2.

yn

h/2 h/2

xn h X n+1

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 46

 Estrategia de ajuste 1
Se calcula el valor yn+1 utilizando un paso de avance h.

Se calcula el valor y*n+1 utilizando dos pasos de avance h/2.

Si el valor absoluto de la diferencia entre ambos valores es

mayor a la tolerancia establecida, entonces se establece h/2
como nuevo valor de h y se vuelve al paso 1.

Si el valor absoluto de la diferencia entre ambos valores es

menor a la mitad de la tolerancia establecida, entonces se
establece 2*h como nuevo valor de h y se vuelve al paso 1.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 47

 Diagrama de Flujo de la estrategia 1

CALCULAR YN+1
CON UN PASO h

CALCULAR Y*N+1
HACER h=h/2 CON DOS PASOS h/2

SI
|YN+1 – Y*N+1| > tol HACER h=2*h

NO

SI
|YN+1 – Y*N+1| < (tol / 2)

NO

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 48

 Variación del error con la Estrategia 1
El paso se reduce de ser necesario, pero nunca se aumenta.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 50

 Variación del paso h con la Estrategia 1
El paso se reduce de ser necesario, pero nunca se aumenta.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 51

 Variación del error con la Estrategia 1
El paso se reduce o se aumenta para lograr la cota de error deseada.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 52

 Variación del paso h con la Estrategia 1
El paso se reduce o se aumenta para lograr la cota de error deseada.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 53

 Estrategia de ajuste 2

A partir de una estimación del error obtenida por medio

de un método como por ejemplo Runge-Kutta-Fehlberg,
podemos ajustar el valor de h, en cada paso.


tol problema
hnuevo =h actual⋅ ∣  RKF ∣
El valor de α se ajusta en base a la relación entre el valor
estimado del error y la tolerancia establecida.
Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 54
 Diagrama de la estrategia 2

SI
tol >= errorRKF

NO
HACER α=0.2
HACER α=0.22


tol
h=h∗∣ ∣
error RKF

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 55

Variación del error con la Estrategia 2
El paso se ajusta de acuerdo al valor de error calculado.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 56

 Variación del paso h con la Estrategia 2
El paso se ajusta de acuerdo al valor de error calculado.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 57

Comparación de ambas estrategias

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 58

 Resumen de la Estrategia de ajuste 1

Ventajas
No se necesita de un método que calcule
explícitamente el error en cada paso.
Los valores de h obtenidos son siempre
múltiplos del h original.

Desventajas
La implementación es un poco más complicada
que la de la estrategia (II).

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 59

 Estrategia de ajuste 2

Ventajas
Su implementación es muy sencilla.

Desventajas
El valor de h que se obtiene generalmente no
es un múltiplo del h original.

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 60

Preguntas . . .

Mg.Ing. Francisco A. Lizarralde Facultad de Ingeniería - UNMDP - 2021 61