1. No category

Computación I
Unidad 5: Integración de las ecuaciones del movimiento
Unidad 5: Integración de las
ecuaciones del movimiento
• Ecuaciones diferenciales ordinarias (EDO).
• Solución EDO de 1er orden (método de Euler,
método de Runge Kutta).
• Solución EDO de 2o orden (método de Euler,
método de Runge Kutta, Verlet).
Aplicaciones en física:
• Integración práctica de la ecuación de Newton.
Ecuaciones del movimiento con condiciones
iniciales. Modelización de sistemas físicos.
1
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales
Definición
Una ecuación diferencial es una relación matemática entre una función f(x) y sus
derivadas en la forma:
𝑑 𝑛 𝑓 𝑑 𝑛−1 𝑓 𝑑 𝑛−2 𝑓
𝑑𝑓
𝐺
,
,
, … , , 𝑓, 𝑥 = 0
𝑑𝑥 𝑛 𝑑𝑥 𝑛−1 𝑑𝑥 𝑛−2
𝑑𝑥
A la ecuación anterior se la llama ecuación diferencial ordinaria (EDO) de orden n,
siendo n el orden de la derivada mayor. La ecuación es “ordinaria” porque solo
depende de una variable independiente x.
Podemos rescribir la ecuación como una ecuación diferencial ordinaria explícita:
𝑑𝑛 𝑓
𝑑 𝑛−1 𝑓 𝑑 𝑛−2 𝑓
𝑑𝑓
=
𝑔
,
,
…
,
, 𝑓, 𝑥
𝑑𝑥 𝑛
𝑑𝑥 𝑛−1 𝑑𝑥 𝑛−2
𝑑𝑥
Ejemplo:
𝑓 (3) + 𝑥𝑓 ′ + 𝑥 −2 𝑓 = 0
𝑓 (𝑛)
𝑑𝑛 𝑓
=
𝑑𝑥 𝑛
2
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales
Definición
El objetivo de este tema es encontrar numéricamente la función f(x) solución de una
ecuación diferencial ordinaria dada (nos centraremos en las de primer y segundo
orden). En general una ecuación diferencial tiene múltiples soluciones, es decir, hay
múltiples funciones f(x) que la cumplen (ejemplo: f’ = 0). Para tener bien definida la
solución el problema necesita condiciones adicionales, a estas condiciones se las
llama condiciones de contorno.
𝑓 (𝑛) = 𝑔 𝑓 (𝑛−1) , 𝑓 (𝑛−2) , … , 𝑓 ′ , 𝑓, 𝑥
𝑓
𝑛−1
𝑓
𝑛−2
(𝑥0 ) = 𝑓0,𝑛−1
(𝑥0 ) = 𝑓0,𝑛−2
…
𝑓′(𝑥0 ) = 𝑓0,1
𝑓(𝑥0 ) = 𝑓0
Ejemplo: f’ = 0, f(1) = 3. Solución, f(x) = 3.
Multitud de problemas en Física son ecuaciones diferenciales.
Ejemplo: Ley de Newton:
𝑑 2 𝑥(𝑡)
𝑚
=𝐹
𝑑𝑡 2
3
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Interpretación geométrica
Ecuación diferencial ordinaria de primer orden con condición de contorno:
𝑑𝑓(𝑥)
= 𝑔 𝑓, 𝑥
𝑑𝑥
𝑓(𝑥0 ) = 𝑓0
La ecuación solo nos dice cual es valor de la pendiente, g(f,x), de la curva tangente a
f(x):
tangente con pendiente g(f,x)
f(x)
𝑓0
𝑥0
4
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Método de Euler
Ecuación diferencial ordinaria de primer orden con condición de contorno:
𝑑𝑓(𝑥)
= 𝑔 𝑓, 𝑥
𝑑𝑥
𝑓(𝑥0 ) = 𝑓0
Usemos esta información para aproximar el valor de f(x) en una sucesión de puntos xi
separados por un intervalo “pequeño” Δx = xi+1 - xi.
Es decir, usamos la definición de derivada:
pendiente g(f1,x1)
pendiente g(f0,x0)
𝑓2
𝑓1
𝑑𝑓(𝑥)
𝑓 𝑥 + ∆𝑥 − 𝑓(𝑥)
= lim
= 𝑔 𝑓, 𝑥
∆𝑥→0
𝑑𝑥
∆𝑥
suponiendo un valor finito de Δx:
f(x)
𝑓 𝑥 + ∆𝑥 = 𝑓(𝑥) + 𝑔 𝑓, 𝑥 ∗ ∆𝑥
𝑓0
𝑓1 = 𝑓0 + 𝑔 𝑓0 , 𝑥0 ∗ (𝑥1 − 𝑥0 )
𝑥0
𝑥1
𝑥2
𝑓2 = 𝑓1 + 𝑔 𝑓1 , 𝑥1 ∗ (𝑥2 − 𝑥1 )
…
𝑓𝑖+1 = 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ (𝑥𝑖+1 − 𝑥𝑖 )
Método
de
Euler
5
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Método de Euler
Algoritmo para implementación del método de Euler:
𝑓𝑖+1 = 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1 − 𝑥𝑖 , 𝑖 = 0, 𝑁 − 1
initial values x0 & f0
end point xN & number of steps N & Δx
function g(f,x)
i<N
No
Plot fi, xi
Plot analytical
result if possible
Yes
𝑥𝑖+1 = 𝑥𝑖 + ∆𝑥
𝑓𝑖+1 = 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1 − 𝑥𝑖
Error: 𝑓 𝑥𝑖 − 𝑓𝑖 ∝ ∆𝑥 ∗ 𝑥𝑖 − 𝑥0
En algunos casos la solución se vuelve inestable.
6
Computación I
Ecuaciones Diferenciales de 1er orden
Unidad 5: Integración de las ecuaciones del movimiento
Método de Euler
Ejercicio 5.1: Ecuación del movimiento en una dimensión. Calcular mediante el método
de Euler la posición en función del tiempo de una partícula que parte en t0 = 0 s del punto
x0 = 0.5 m y se mueve con una velocidad v(t) = ( 3 – 2 t ) m/s. Comparar la solución
numérica con la exacta.
Ejercicio 5.2: Escribir una función de matlab que resuelva mediante el método de Euler
una ecuación diferencial ordinaria de primer orden, 𝑑𝑓(𝑥)/𝑑𝑥 = 𝑔 𝑓, 𝑥 , con 𝑓(𝑥0 ) =
𝑓0 para N puntos y en el intervalo [𝑥0 , 𝑥𝑁 ]. La función tiene que ser del tipo:
function [f, x] =euler(g,x0,f0,xN,N),
siendo f,x dos vectores con la solución 𝑓𝑖 , 𝑥𝑖 y g la función que define la ecuación
diferencial. Resolver el caso anterior usando esta función.
Ejercicio 5.3: Calcular la velocidad en función del tiempo en un intervalo de 4s de una
partícula que se mueve bajo una fuerza F(t)/m = 1 / ( 3 – 2 t ) m/s2 si i) en t0 = 0 s, v0 = 0.5
m/s; y ii) en t0 = 1.6 s, v0 = 0.5 m/s.
Ejercicio 5.4: Calcular el error cometido por el método de Euler al resolver la ecuación
diferencial 𝑑𝑓(𝑥)Τ𝑑𝑥 = −𝑓 + 2cos(𝑥) con 𝑓 0 = 1 para x de 0 a 10 y i) Δx=0.1, ii)
Δx=0.05. (Solución analítica: 𝑓 𝑥 = sin 𝑥 + cos(𝑥))
7
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Método de Euler modificado
Euler:
pendiente g(fi+1,xi+1)
pendiente g(fi,xi)
𝑓𝑖+2
𝑓𝑖+1
f(x)
𝑓𝑖
𝑥𝑖
𝑥𝑖+1
𝑓𝑖+1 = 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1 − 𝑥𝑖
¿Por qué usar 𝑔 𝑓𝑖 , 𝑥𝑖 ?
Podríamos usar 𝑔 𝑓𝑖+1 , 𝑥𝑖+1 … o su media.
Problema, no conocemos todavía 𝑓𝑖+1 .
Lo solucionamos haciendo una estimación
con el método de Euler:
𝑥𝑖+2
𝑡𝑟𝑖𝑎𝑙
𝑓𝑖+1
= 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1 − 𝑥𝑖
𝑡𝑟𝑖𝑎𝑙
𝑔 𝑓𝑖 , 𝑥𝑖 + 𝑔 𝑓𝑖+1
, 𝑥𝑖+1
𝐺=
2
Método de Euler
modificado
𝑓𝑖+1 = 𝑓𝑖 + 𝐺 ∗ 𝑥𝑖+1 − 𝑥𝑖
8
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Método de Runge-Kutta de orden 2
pendiente g(fi+1,xi+1)
pendiente g(fi,xi)
𝑓𝑖+2
𝑓𝑖+1
f(x)
𝑓𝑖
𝑥𝑖
𝑥𝑖+1
También podríamos usar la pendiente en el
punto medio 𝑔 𝑓𝑖+1/2 , 𝑥𝑖+1/2
𝑓𝑖+1/2 lo estimamos con Euler.
1
Notación: 𝑥𝑖+1/2 = 𝑥𝑖 + 2 ∆𝑥𝑖 ,
𝑓𝑖+1/2 = 𝑓(𝑥𝑖+1/2 )
𝑥𝑖+2
𝑓𝑖+1/2 = 𝑓𝑖 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1/2 − 𝑥𝑖
𝐺=𝑔
𝑡𝑟𝑖𝑎𝑙
𝑓𝑖+1/2
, 𝑥𝑖+1/2
𝑓𝑖+1 = 𝑓𝑖 + 𝐺 ∗ 𝑥𝑖+1 − 𝑥𝑖
Método de RungeKutta de orden 2
(midpoint)
El método de Euler modificado también es un RungeKutta orden dos (trapezoidal). Se llaman de orden dos
porque están desarrollados para que el error sea ∝ ∆𝑥 2
9
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 1er orden
Método de Runge-Kutta de orden 4
Si se impone que el error cometido sea ∝ ∆𝑥 4 se obtiene el método de Runge-Kutta
de orden 4.
𝑘1 = 𝑔 𝑓𝑖 , 𝑥𝑖
𝑘2 = 𝑔 𝑓𝑖 + 𝑘1 ∆𝑥 Τ2 , 𝑥𝑖 + ∆𝑥 Τ2
𝑘3 = 𝑔 𝑓𝑖 + 𝑘2 ∆𝑥 Τ2 , 𝑥𝑖 + ∆𝑥 Τ2
𝑘4 = 𝑔 𝑓𝑖 + 𝑘3 ∆𝑥, 𝑥𝑖 + ∆𝑥
Método de RungeKutta de orden 4
𝐺 = 1Τ6 (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 )
𝑓𝑖+1 = 𝑓𝑖 + 𝐺 ∗ 𝑥𝑖+1 − 𝑥𝑖
10
Computación I
Ecuaciones Diferenciales de 1er orden
Unidad 5: Integración de las ecuaciones del movimiento
Ejercicios
Ejercicio 5.5: Modificar la función euler creada anteriormente para que resuelva
ecuaciones diferenciales con el método de Euler modificado. Repetir el ejercicio 5.4 con
este método. Analizar el comportamiento del error.
Ejercicio 5.6: Modificar la función euler creada anteriormente para que resuelva
ecuaciones diferenciales con el método de Runge-Kutta orden 2 (midpoint). Repetir el
ejercicio 5.4 con este método. Analizar el comportamiento del error.
Ejercicio 5.7: Modificar la función euler creada anteriormente para que resuelva
ecuaciones diferenciales con el método de Runge-Kutta orden 4. Repetir el ejercicio 5.4
con este método. Analizar el comportamiento del error. Comparar todos los métodos
estudiados.
Ejercicio 5.8: Movimiento de un cuerpo en caída libre. Calcular mediante el método de
Euler la velocidad en función del tiempo de un cuerpo de masa m = 1 g que parte del
reposo y además de la fuerza de la gravedad sufre una fuerza de rozamiento proporcional
a la velocidad, -kv, con k = 10-4 Ns/m. Repetir el cálculo usando el método de Runge-Kutta
de orden 4 y el mismo intervalo y número de pasos. Comparar ambas soluciones con el
resultado exacto.
11
Computación I
Ecuaciones Diferenciales de 1er orden
Unidad 5: Integración de las ecuaciones del movimiento
Ejercicios
Ejercicio 5.9: Suponer que el cuerpo del ejercicio anterior se lanza desde el origen de
coordenadas con una velocidad inicial de módulo v0 = 1 m/s formando un ángulo de 45o
con la horizontal. Representar la velocidad y posición en función del tiempo, también la
trayectoria y(x). Comparar la solución exacta con la numérica.
Ejercicio 5.10: Método de Verlet (métodos multistep). Los métodos multistep usan varios
“pasos” de la sucesión de puntos fi, xi para calcular el siguiente “paso” fi+1. El más sencillo
es el método de Verlet: 𝑓𝑖+1 = 𝑓𝑖−1 + 𝑔 𝑓𝑖 , 𝑥𝑖 ∗ 𝑥𝑖+1 − 𝑥𝑖−1 . a) (Avanzado) Deducir el
método a partir del desarrollo de Taylor de la función f(x) (Ayuda: utilizar el desarrollo
hasta orden dos para aproximar 𝑓𝑖+1 y 𝑓𝑖−1 y combinar ambas ecuaciones para cancelar el
término proporcional a ∆𝑥 2 . b) Usar el método de Verlet para resolver el ejercicio 5.4
(Ayuda: para iniciar el algoritmo hace falta conocer no solo la condición inicial, 𝑓(𝑥0 ) =
𝑓0 , si no también el siguiente paso 𝑓1 , estimarlo usando Euler).
12
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 2º orden y
sistemas de ecuaciones diferenciales.
𝑑 2 𝑓(𝑥)
′
=
𝑔
𝑓
, 𝑓, 𝑥
𝑑𝑥 2
EDO 2º orden
𝑓(𝑥0 ) = 𝑓0
𝑓′(𝑥0 ) = 𝑓′0
Vamos a resolver este problema convirtiéndolo en un sistema de dos ecuaciones
diferenciales de primero orden. Para entender mejor el procedimiento vamos a resolver el
problema de un oscilador armónico:
𝑥ሶ = 𝑣
Con las condiciones de contorno
𝑘
𝑥ሷ + 𝑥 = 0
𝑘
𝑥 𝑡0 = 𝑥0 , 𝑣 𝑡0 = 𝑣0
𝑚
𝑣ሶ = − 𝑥
𝑚
Lo que queremos calcular numéricamente es el conjunto (vector) de funciones 𝑥(𝑡), 𝑣(𝑡).
Tenemos que:
𝑑 𝑥(𝑡)
𝑣
𝑔1 (𝑥, 𝑣, 𝑡)
=
= − 𝑘Τ𝑚 𝑥
𝑔2 (𝑥, 𝑣, 𝑡)
𝑑𝑡 𝑣(𝑡)
Para resolverlo aplicaremos al mismo tiempo a x y v los métodos de resolución de EDO’s
de primer orden explicados anteriormente.
13
Unidad 5: Integración de las ecuaciones del movimiento
Computación I
Ecuaciones Diferenciales de 2º orden y
sistemas de ecuaciones diferenciales.
𝑥ሶ = 𝑣
𝑘
𝑥ሷ + 𝑥 = 0
𝑚
Solución por Euler:
Con las condiciones de contorno
𝑥 𝑡0 = 𝑥0 , 𝑣 𝑡0 = 𝑣0
𝑘
𝑣ሶ = − 𝑥
𝑚
𝑥𝑖+1 = 𝑥𝑖 + 𝑣𝑖 ∆𝑡
𝑣𝑖+1
En general:
Oscilador armónico
𝑘
= 𝑣𝑖 − 𝑥𝑖 ∆𝑡
𝑚
𝑑 2 𝑓(𝑥)
= 𝑔 𝑓 ′ , 𝑓, 𝑥
2
𝑑𝑥
∈
𝑥𝑖+1 = 𝑥𝑖 + 𝑔1 (𝑣𝑖 , 𝑥𝑖 , 𝑡)∆𝑡
𝑣𝑖+1 = 𝑣𝑖 + 𝑔2 (𝑣𝑖 , 𝑥𝑖 , 𝑡)∆𝑡
𝑓(𝑥0 ) = 𝑓0
𝑑𝑓(1) (𝑥)
= 𝑔1 𝑓(1) , 𝑓(2) , 𝑥 = 𝑓(2)
𝑑𝑥
𝑑𝑓(2) (𝑥)
= 𝑔2 𝑓(1) , 𝑓(2) , 𝑥 = 𝑔 𝑓 ′ , 𝑓, 𝑥
𝑑𝑥
𝑓′(𝑥0 ) = 𝑓′0
con:
𝑓 1 𝑥 = 𝑓(𝑥)
𝑑𝑓(𝑥)
𝑓2 𝑥 =
𝑑𝑥
14
Computación I
Unidad 5: Integración de las ecuaciones del movimiento
Ecuaciones Diferenciales de 2º orden
Oscilador armónico
Ejercicio 5.11: Calcular numéricamente mediante Euler la solución de un oscilador
armónico con k = 1 kg/s2, m = 2 kg, t0 = 0 s x0 = -1 m, v0 = 2 m/s. Representar las
soluciones numéricas y analíticas en el intervalo de tiempo [0s,80s] (usar N=1000).
Calcular y representar la energía total.
Ejercicio 5.12: Repetir el ejercicio anterior usando el método de Runge-Kutta de orden 2
(midpoint). Comparar las dos soluciones numéricas y la analítica.
Ejercicio 5.13: Oscilador armónico amortiguado. Añadir al problema anterior una fuerza
de fricción proporcional a la velocidad, -bv, con b = 0.1 Ns/m y resolver el problema en el
mismo intervalo temporal.
Ejercicio 5.14: Resolver numéricamente en el intervalo [0,10] la ecuación:
𝑑2𝑦
+ 𝑦2 − 𝑥 = 0
2
𝑑𝑥
con las condiciones de contorno y(0) = 0, y’(0) = 1
15
Computación I
Ecuaciones Diferenciales de 2º orden
Unidad 5: Integración de las ecuaciones del movimiento
Oscilador armónico
Ejercicio 5.15: Método de Verlet para fuerzas conservativas. La ecuación de movimiento
de Newton 𝑑 2 𝑥(𝑡)Τ𝑑𝑡 2 = 𝑎 𝑡 = 𝐹(𝑥, 𝑡)/𝑚 para fuerzas no dependientes de la
velocidad se puede resolver usando el siguiente algoritmo: 𝑥𝑖+1 = 2𝑥𝑖 − 𝑥𝑖−1 + 𝑎 𝑡𝑖 ∗
(∆𝑡)2 . Además de la condición inicial, 𝑥0 , 𝑣0 , necesitamos conocer la posición 𝑥1 . La
estimamos usando: 𝑥1 = 𝑥0 + 𝑣0 ∆𝑡 + 1Τ2 𝑎(𝑡0 )(∆𝑡)2 . Demostrar que este algoritmo
conserva la energía mecánica del sistema aplicándolo a la resolución del ejercicio 11.
Ejercicio 5.16: Contrariamente a lo que ocurre en la mecánica newtoniana, en relatividad
especial la aceleración de un objeto depende del sistema de referencia inercial en el que
se mide. Supongamos que un astronauta experimenta una aceleración constante igual a g
en su sistema en reposo instantáneo y que parte del reposo desde la Tierra, se pude
demostrar que 𝑎𝑥 𝑡 = 𝑑𝑣Τ𝑑𝑡 = 𝑔(1 − 𝑣 2 Τ𝑐 2 )3Τ2 .
Calcular su velocidad (medida desde la Tierra) en un instante de tiempo t (terrestre).
Calcular también qué distancia ha recorrido al cabo de un tiempo terrestre t. ¿Es posible
que supere la velocidad de la luz al cabo de un cierto tiempo?
16
Computación I
Ecuaciones Diferenciales de 2º orden
Unidad 5: Integración de las ecuaciones del movimiento
Cargas
Ejercicio 5.17: El movimiento en 1Dimension de una partícula cargada bajo la acción de
una fuerza de Coulomb no tiene solución exacta. Una carga puntual Q = 1μC se encuentra
en el origen de coordenadas. Otra carga puntual q = 1nC de masa m = 1g se lanza contra
Q en t = 0, desde x0 = -1m con velocidad v0 = 0.1m/s. Calcular numéricamente su posición
y velocidad como función del tiempo durante 20s.
Nota: Aplicando el principio de conservación de la energía, calcular analíticamente la
distancia mínima y la velocidad máxima que q puede alcanzar. Analizar si la solución
numérica cumple estos criterios.
17
Computación I
Unidad 5: Integración de las ecuaciones del movimiento
Ecuaciones Diferenciales de 2º orden
Péndulo
Ejercicio 5.18: Tenemos un péndulo de longitud L que parte del reposo desde un ángulo
inicial 𝜃0 y con una masa M (que supondremos puntual). El péndulo tiene su soporte
suspendido a cierta altura H < L de un fluido, de tal forma que al oscilar parte del
recorrido lo realiza sumergido en este. Mientras se encuentra en el fluido el objeto
experimenta una fuerza de rozamiento FR = − μ0 v , siendo la velocidad v = L dθ/dt y μ0 es
el coeficiente de rozamiento. Cuando no está en el fluido, la masa M no experimenta
ningún tipo de rozamiento.
La ec. del movimiento es 𝑀 𝑑 2 𝜃Τ𝑑𝑡 2 = − 𝑔𝑀Τ𝐿 × 𝑠𝑒𝑛𝜃 − 𝜇 𝑑𝜃Τ𝑑𝑡,
con 𝜇 = ቊ
𝜇0 dentro del fluido, 𝜃 ≤ 𝛼
0 fuera del fluido, 𝜃 > 𝛼
Realizar un script de nombre que para tres
oscilaciones completas:
A. Realice la integración de la ecuación del
movimiento.
B. Dibuje en un gráfico la posición angular 𝜃 en
función del tiempo y en otro diferente la
energía cinética, la potencial gravitatoria y la
total en función del tiempo.
Datos: L=1.0 m, H=0.6 m, M=0.1 Kg,
𝜇0 =0.025 Kg/s ,𝜃0 =1.5 rad.
18
Computación I
Unidad 5: Integración de las ecuaciones del movimiento
19