Ecuaciones Diferenciales Ordinarias
Motivación
● Problema del paracaidista
d vd t= g−
cm
v t
● v: variable dependiente (incógnita)
● t: variable independiente (dato)
● EDOs: ● De primer orden● De segundo orden
d vd t= g−
cm
v t
md 2 x
d t 2 cd xd tk x=0
EDOs de segundo orden
● Se transforman en un sistema de EDOs de primer orden con una sustitución:
y=d xd t
⇒d yd x=d 2 x
d t 2
● reemplazando
md yd tc yk x=0 ⇒
d yd t=−
c yk xm
● El sistema,d xd t= y
d yd t=−
c yk xm
Solución de EDOs sin computadora
● En algunos casos se obtiene por integración indefinida:
d vd t= g−
cm
v t ⇒ v=∫ g− cm
v t d t ⇒ v t =g mc
1−e−
cm
t ● Técnica habitual: linealización
● Ejemplo: péndulo● EDO original (no lineal)
d 2
d t 2 gl
sin=0
● Si θ es pequeño (EDO lineal)d 2
d t 2 gl=0
Ejemplos de EDOs en Ingeniería● Segunda ley de Newton del movimiento:
d vd t=
Fm
● Ley del calor de Fourier:
q=−kd Td x
● Ley de difusión de Fick:
J=−Dd Cd x
● Ley de Faraday:
V=Ld id t
Solución de problemas
Antecedentes matemáticos
● Solución de una EDO: función de la v.i. Y de las condiciones iniciales.
● Supongamos
y=−0.5 x44 x3−10 x28.5 x1
● Derivando,
d yd x=−2 x3
12 x2−20 x8.5
● Graficando,
Antecedentes matemáticos
Antecedentes matemáticos
● Ahora suponemos que debemos resolver la EDO:
d yd x=−2 x3
12 x2−20 x8.5
● Integrando,
y=∫ −2 x312 x2−20 x8.5 dx
y= ...=−0.5 x44 x3
−10 x28.5 xC
Solución general
● Cuál de todas las curvas es? Falta más información...
Antecedentes matemáticos
● ... por ejemplo saber que cuando x = 0, y = 1 (condición inicial). Entonces:
y=−0.5 x44 x3−10 x28.5 x1Solución particular
● En general, una EDO de orden n requiere n condiciones:● Si todas se fijan en el mismo punto (p. ej. x = 0) es
un problema de valores iniciales● Si se fijan en distintos puntos, es un problema de
contorno
Orientación
Métodos de Runge-Kutta
● Solución de EDOs de la formad yd x= f x , y
● El problema del paracidista se resolvió como
Nuevo valor = valor anterior + pendiente x paso
y i1= y i h
=d yd x
.... donde??
Método de Euler
● La pendiente se estima como:
= f x i , y i
● Por lo tanto,
y i1= y i f x i , y i h
Ejemplo 25.1 pag. 720
● Resolver, para x entre 0 y 4 con h = 0.5, d yd x=−2 x3
12 x2−20 x8.5 , y 0=1
● solución exacta: y=−0.5 x44 x3−10 x28.5 x1
● Solución. En el primer paso,y 0.5= y 0 f 0,1×0.5=1.08.5×0.5=5.25
● La solución verdadera esy 0.5=−0.5×0.54
4×0.53−10×0.52
8.5×0.51=3.21875
● El error E t=3.21875−5.25=−2.03125
t=−63.1%
Ejemplo 25.1 pag. 720
● Solución completa: 25_1.ods
Análisis del error en el método de Euler
● Errores de● Truncamiento
– Local– Global
● Redondeo
● Serie de Taylor
y i1= y i y ' i hy ' ' i2!
h2
y iIII
3!h3 ...
y in
n!hnRn
● Es deciry i1= y i f x i , y i h
Euler
f ' x i , y i
2!h2 ...
error
⇒ Ea=O h2
Error local
Ejemplo 25.2
● Calcular el error del método de Euler en el primer paso del ejemplo 25.1
● Solución E t=f ' x i , y i
2!h2
f ' ' x i , y i
3!h3
f III x i , y i
4!h4
● dondef ' x i , y i =−6 x2
24 x−20 ⇒ f ' 0,1=−20
f ' ' x i , y i=−12 x24 ⇒ f ' ' 0,1=24
f III x i , yi =−12
● luego
E t=−202
0.52
246
0.53−
1224
0.54=−2.50.5−0.03125=−2.03125
Análisis del error en el método de Euler
● Error local O(h²) --> error global O(h)● El error se reduce reduciendo h● El método es exacto para y lineal (f' = 0)
Ejemplo 25.3
● Repita el ejemplo 25.1 con h = 0.25
● Solución: 25_3.ods
Algoritmo del método de Euler
● Código en Octave: euler.m
Ejemplo 25.4
● Resolver el problema del paracaidistad vd t= g−
cm
v t lineal
d vd t= g−
cm [ va v
vmax b
] No lineal
● Donde: a = 8.3 , b = 2.2 y vmax
= 46 son
constantes empíricas
● Datos: m = 68.1 kg ; c = 12.5 kg/s; v = 0 en t = 0
● Código en Octave: p25_4.m
Método de Heun
● Se predice por Euler:
y i10 = y i f x i , y ih
predictor
● Para estimar
y ' i1= f x i1 , y i10
● Que sirve para calcular
y i1= y iy ' i y ' i1
2h
corrector
y i1= y if x i , y i f x i1 , y i1
0
2h
Método de Heun
● Se puede plantear el esquema predictor – corrector en forma iterativa, hasta que
∣ t∣=∣ y i1j − y i1
j−1
y i1j ∣100% s
Ejemplo 25.5● Resolver, con h = 1
y '=4 e0.8 x−0.5 y , 0≤ x≤4 , y 0=2● Solución analítica:
y=4
1.3 e0.8 x
−e−0.5 x 2 e−0.5 x
● Primer paso: y ' 0=4 e0−0.5×2=3
predictor y10=23×1=5
y ' 1= f x1, y10=4 e0.8×1−0.5×5=6.402164
correctory1=2
36.4021642
×1=6.701082
● Solución completa: 25_5.ods
Error del método de Heun● Supongamos que y' = f(x,y) = f(x) entonces
y i1= y if x i f x i1
2h
regla del trapecio
d yd x= f x ⇒ ∫
y i
yi1
d y=∫x i
x i1
f x d x ⇒ y i1− y i=∫x i
x i1
f x d x
Por la regla del trapecio,
∫x i
x i1
f x d x=f x i f x i1
2h−
f ' ' 12
h3
y i1= y if x i f x i1
2hO h3 Error local
--> error global O(h²)
Error del método de Heun
d yd x=−2 x3
12 x2−20 x8.5Solución de
Método del punto medio
● Se calcula
y i1 / 2= y i f x i , y ih2
y ' i1 / 2= f x i1 / 2 , y i1 / 2
y i1= y i f x i1 / 2 , y i1 / 2 h
● Se basa en la fórmula de integración del punto medio:
∫x i
xi1
f x d x=h f x i1 / 2
Error local: O(h³)Error global: O(h²)
Algoritmos de los métodos de Heun y del punto medio
● Códigos en Octave:● heun.m● puntomedio.m
● Se resuelve el problema
y '=−2 x− y , y 0=1 , 0≤ x≤1
Métodos de Runge - Kutta
● Forma general y i1= y i x i , y i , hfuncion incremento
h
=a1k1a2k2 ...an kn
k1= f x i , y i
k2= f x i p1h , y iq11 k1h
k3= f x i p2h , y iq21k1hq22 k2h...
kn= f x i pn−1h , yi∑j=1
n−1
qn−1, j k j h
Métodos de Runge – Kutta de segundo orden
● Forma general
y i1= y ia1k1a2k2h
k1= f x i , y i
k2= f x i p1h , y iq11 k1h
● De la serie de Taylor, y i1= y i f x i , y i hf ' x i , y i
2!h2
f ' x i , y i=∂ f x i , y i
∂ x∂ f x i , y i
∂ yd yd x
Regla de la cadena
● reemplazando, y i1= y i f x i , y i h ∂ f∂ x∂ f∂ y
f h2
2!
● De la serie de Taylor, g x r , y s= g x , y r∂ g∂ x s∂ g∂ y ...
Métodos de Runge – Kutta de segundo orden
● Se tiene,
f x i p1h , y iq11 k1h= f x i , y i p1h∂ f∂ xq11k1h
∂ f∂ yO h2
● reemplazando,
y i1= y ia1h f x i , y i a2 p1h2 ∂ f∂ xa2q11h
2 f x i , y i∂ f∂ yO h3
● reordenando,
y i1= y i a1a2 f x i , y i h [ a2 p1∂ f∂ xa2q11 f x i , y i
∂ f∂ y ] h2
O h3
● comparando, a1a2 =1 , a2 p1=12
, a2q11=12
● Es decir a1=1−a2 , p1=q11=1
2a2
Métodos de Runge – Kutta de segundo orden
● Si a2=12⇒ a1=
12⇒ p1=q11=1
se tiene y i1= y i 12k1
12k 2 h
con k1= f x i , y i , k 2= f x ih , yik1h Método de Heun
● Si a2=1 ⇒ a1=0 ⇒ p1=q11=12
se tiene y i1= y ik 2h
con k1= f x i , y i , k 2= f x i12h , y i
12k1h Método del punto
medio
Método de Ralston
● Si a2 = 2/3, se minimiza el error de
truncamiento.
● Si a2=23⇒ a1=
13⇒ p1=q11=
34
se tiene y i1= y i 13k1
23k 2 h
k1= f x i , y i , k 2= f x i34h , y i
34k1h Método de Ralston
Ejemplo 25.6 pag. 744
● Comparación de varios esquemas de RK de segundo orden d y
d x=−2 x3
12 x2−20 x8.5 , y 0=1
● Solución en 25_6.ods
Métodos de Runge-Kutta de tercer orden
● Una versión común es
y i1= y i16 k14 k2k3 h
● donde
k2= f x i12h , y i
12k1h
k3= f x ih , y i−k1h2k2h
k1= f x i , y i
● Si f = f(x) se transforma en la Regla 1/3 de Simpson --> error local O(h4) y error global O(h³)
Métodos de Runge-Kutta de cuarto orden
● Una versión común es
y i1= y i16 k12k22 k3k4 h
● donde
k2= f x i12h , y i
12k1h
k1= f x i , y i
k3= f x i12h , y i
12k2h
k4= f x ih , y ik3h
Ejemplo 25.7 pag. 747
● Resolver, con h = 0.5d yd x=−2 x3
12 x2−20 x8.5 , y 0=1
● Se calculan (...):
k1=8.5 , k2=4.21875 , k3=4.21875 , k4=1.25
● reemplazando:
y 0.5=1168.52×4.218752×4.218751.25 0.5=3.21875
Solución exacta
Ejemplo 25.7 pag. 747● Resolver, con h = 0.5
y '=4 e0.8 x−0.5 y , 0≤ x≤0.5 , y 0=2
● Solución.k1= f 0,2=4 e0.8×0
−0.5×2=3 ⇒ y 0.25=23×0.25=2.75
k 2= f 0.25,2.75=4 e0.8×0.25−0.5×2.75=3.510611
y 0.25=23.510611×0.25=2.877653k 3= f 0.25,2.877653=4 e0.8×0.25
−0.5×2.877653=3.446785
y 0.5=23.446785×0.5=3.723392k 4= f 0.5,3.723392=4 e0.8×0.5
−0.5×3.723392=4.105603
=16[ 32×3.5106112×3.4467854.105603 ]=3.503399
y 0.5=23.503399×0.5=3.751669 Exacta : 3.751521
Métodos de Runge-Kutta de orden superior● Método RK de quinto orden de Butcher:
y i1= y i1
90 7 k132k 312k 432 k57 k6 h
● con k1= f x i , y i
k 2= f x i14h , y i
14k1h
k 3= f x i14h , y i
18k1h
18k 2h
k 4= f x i12h , y i−
12k 2hk 3h
k 5= f x i34h , y i
316
k 1h9
16k 4h
k 6= f x ih , y i−37k1h
27k 2 h
127k 3h−
127
k 4h87k 5h
Algoritmo de los métodos de Runge-Kutta
● Código en Octave para el método de RK de 4to orden
● rk4.m
Sistemas de ecuaciones
● Requiere n condiciones iniciales en x0
d y1
d x= f 1 x , y1, y2, ... , yn
d y2
d x= f 2 x , y1, y2, ... , yn
...d yn
d x= f n x , y1, y2, ... , yn
Ejemplo 25.9, pag. 752
● Resolver el sistema de ecuacionesd y1
d x=−0.5 y1 ;
d y2
d x=4−0.3 y2−0.1 y1
● Por el método de Euler, con h = 0.5, sabiendo que en x = 0, y
1 = 4 y y
2 = 6.
● Solución en planilla de cálculo: 25_9.ods
Ejemplo 25.10 pag. 752
● Resolver por el método de Runge-Kutta de 4to orden el problema anterior.
● Solución en planilla de cálculo: 25_10.ods
Métodos adaptativos
● Los métodos de paso constante pueden ser ineficientes
● Se puede “adaptar” el paso estimando el error de truncamiento local en cada paso.
● También se pueden aplicar al cálculo de integrales
Método adaptativo de RK o de mitad de paso
● Cada paso se calcula dos veces:
● De un solo paso h, obteniendo y1.
● En dos pasos h/2, obteniendo y2.
● El error Δ se representa por: Δ = y2 - y
1
● Se realiza la corrección y2 = y
2 + Δ/15
● Exactitud de quinto orden
Problemas 25.1 a 25.26, pag. 764
Problemas de valores en la frontera (problemas de contorno) y de
valores propios● Problemas de EDOs:
● Problemas de valores iniciales
● Problemas de contorno– Valores propios,
autovalores o eigenvalores
Ejemplo● Distribución de
temperaturas en una barra uniforme no aislada
d 2T
d x2h' T a−T =0
● Condiciones de contorno
T 0=T 1
T L=T 2
● Si L = 10 m, Ta = 20 °C, T
1 = 40 °C, T
2 = 200 °C y
h' = 0.01 m-2, la solución analítica es:
T=73.4523 e0.1 x−53.4523 e−0.1 x20
El método del disparo
● Consiste en convertir el problema de contorno en un problema de valores iniciales.
● Se realizan dos estimaciones y luego una interpolación.
Ejemplo 27.1
● Resolver con
se transforma en un sistema de EDOs de primer orden:
d 2T
d x2h' T a−T =0
T 0=40T 10=200
d Td x= z ,
d zd x=h' T−T 0
Suponemos un valor inicial para z, z(0) = 10
Aplicando RK4 con h = 2, (27_1.ods) se obtiene T(10) = 168.3797
Suponemos z(0) = 20, y con lo que se obtiene (27_1.ods) T(10) = 285.8980
Ejemplo 27.1
● Como la EDO es lineal, se interpola linealmente el valor de z(0) para obtener T(10) = 200:
z 0=200−168.3797
285.8980−168.379720
200−285.8980168.3797−285.8980
10
z 0=12.6907
● Con este valor se resuelve el sistema de EDOs de primer orden para determinar T a lo largo de la barra
Problemas no lineales
● Se reformula como un problema de raíces, considerando que
T 10= f z0
● Se busca la raíz deg z0= f z0−200
Ejemplo 27.2
● Suponer la siguiente EDO no lineal para la barra calentada, con h'' = 5 x 10-8:
d 2T
d x2h' ' T a−T
4=0
● Se reduce a un sistema de EDOs de 1er orden:
d Td x= z ,
d zd x=h' T−T 0
4
● Solución en: 27_2.ods
Ejemplo 27.2
Método de diferencias finitas*
● Se sustituyen las derivadas por diferencias finitas divididas
● La EDO se transforma en un sistema de ecuaciones algebraicas● Lineales, si la EDO es lineal● No lineales, si la EDO es no lineal
● Para resolverlo se aplican los métodos de la parte 3
Ejemplo 27.3
● Resolver por diferencias finitas, el problema de la transmisión de calor en la barra
● Se sustituyen las derivadas:
d 2T
d x2h' T a−T =0 ⇒
T i1−2T iT i−1
x2−h' T i−T a=0
Ecuación sustituta o molécula de cálculo
● reordenando:
−T i−1 2h' x2 T i−T i1=h' x2T a
Ejemplo 27.3
● Se aplica la molécula de cálculo a cada nodo:
● Se genera un sistema de n-1 ecuaciones:● Cada una con 3 incógnitas,● Excepto la primera y la última con 2 incógnitas
Ejemplo 27.3
● Eligiendo Δx = 2, se tiene:● Nodo 1:● Nodo 2:● Nodo 3:● Nodo 4:
● En forma matricial
2.04T 1−T 2=0.01 x 2 x 4040=40.8
−T 12.04T 2−T 3=0.8
−T 22.04T 3−T 4=0.8
−T 32.04T 4=0.8200=200.8
[2.04 −1 0 0−1 2.04 −1 00 −1 2.04 −10 0 −1 2.04 ] {
T 1
T 2
T 3
T 4} = {
40.80.80.8
200.8}
Ejemplo 27.3
● Sistema tridiagonal → algoritmo de Thomas● Solución● Algoritmo en Octave: diferencias.m
{T }T= [ 65.9698 93.7785 124.5382 159.4795 }
Condiciones de borde en la derivada*
● Resolver● Para● Con
● Ahora hay n+1 incógnitas (incluido yn+1
):
y ' '−2 y ' y= x
0.0≤ x≤1.0
h=0.1 , y 0.0=−1, y ' 1.0=1
Condiciones de borde en la derivada*
● La condición de borde en 1.0 se iguala a una diferencia finita dividida centrada:
y ' 1.0= y ' n=yn1− yn−1
2h=1
● La (n+1)-ésima ecuación surge de la anterior :
yn1− yn−1=2h
Problemas de valores propios
● El sistema tiene solución única si det(A) ≠ 0
● El sistema tiene solución única trivial si det(A) ≠ 0, y soluciones no triviales (∞) si det(A) = 0
● Los problemas de valores propios son de la forma:
[ A ] {X }={B}
[ A ] {X }=0
a11− x1 a12 x2 ...a1n xn=0a21 x1 a22− x2 ... a2n xn=0
........... ..........an1 x1 an2 x2 ... ann− xn=0
Problemas de valores propios
● O más brevemente,● Los valores λ que hacen
se denominan valores propios del sistema, y su solución vector propio
[ [ A]− [ I ] ] {X }=0
det [ [ A]− [ I ] ]=0
Ejemplo: oscilación masa-resorte
● Sabiendo que
m1
d 2 x1
d t2=−k1 xk x2− x1
m2
d 2 x2
d t2=−k x2− x1−k x2
x i=Ai sin t
x ' ' i=−Ai2 sin t
● Se llega a
[2 km1
− 2 −km1
−km2
2km2
−2 ] { A1
A2}= { 00 }
Ejemplo: pandeo
● Deformación de la columna:
● Donde M = -P.y● Reemplazando,
d 2 y
d x2=
ME I
d 2 y
d x2 p2 y=0
● con p2=
PE I
● Para y(0) = y(L) = 0
Ejemplo: pandeo
● Solución analítica: ● Para x = 0,● Para x = L,● Reemplazando(...),
y=Asin pxB cos px
y=Asin 0B cos0=0 ⇒ B=0
y=Asin pL=0 ⇒ pL=n , n=1, 2, 3,...
P=n2
2 E I
L2
● Si n = 1 (primer modo):
P=
2 E I
L2 Fórmula de Euler
Ejemplo: pandeo
● Ej. 27.5● Datos:
● E = 10 x 109 Pa● I = 1.25 x 10-5 m4
● L = 3 m
1 1,0472 137,0782 2,0944 548,3113 3,1416 1233,7014 4,1888 2193,2455 5,2360 3426,9466 6,2832 4934,8027 7,3304 6716,8148 8,3776 8772,982
n p (m^-2) P (kN)
Método del polinomio
● Se reemplaza la EDO por una ecuación en diferencias finitas, se iguala a 0 el determinante de la matriz de coeficientes y se resuelve el polinomio resultante.
● Para el ejemplo anterior, se tiene
● reordenando,
y i1−2 y i y i−1
h2 p2 y i=0
y i−1−2−h2 p2 y i y i1=0
Método del polinomio
● Con 5 tramos (4 nodos interiores) se tiene:
[2−h2 p2 −1 0 0
−1 2−h2 p2 −1 0
0 −1 2−h2 p2 −1
0 0 −1 2−h2 p2] {
y1
y2
y3
y4} =0
Ejemplo 27.6
● Determinar los valores propios del ej. 27.5, con a) 1, b) 2, c) 3, y d) 4 nodos interiores.
● a) h = 3/2−2−2.25 p2 y i=0
det A=2−2.25 p2=0
p=±0.9428 ⇒ t≈10%
Ejemplo 27.6
● b) h = 3/3 = 1
[ 2− p2−1
−1 2− p2 ] { y1
y2}=0
det A= 2− p2 2−1=0
p=±1 ⇒ t≈4.5%
p=±1.73205 ⇒ t≈17%
Ejemplo 27.6
● c) h = 3/4
[2−0.5625 p2
−1 0−1 2−0.5625 p2
−10 −1 2−0.5625 p2 ] {
y1
y2
y3}=0
det A= 2−0.5625 p2 3−2 2−0.5625 p2 =0
p=±1.0205 ⇒ t≈2.5%
p=±1.8856 ⇒ t≈10%
2−0.5625 p2=0
2−0.5625 p2= 2
p=±2.4637 ⇒ t≈22%
Ejemplo 27.6
El método de potencias
● El problema de valores propios se escribe como:
[ A] {X }= {X }
● Puesto en forma iterativa:
[ A] {X k−1}= k {X k }
● El método permite obtener el mayor valor propio
Ejemplo 27.7
● Determine el valor propio mayor del punto c) del ejemplo 27.6
● Se escribe como:
3.5556 x1 −1.7778 x2 −1.7778 x1 3.5556 x2 −1.7778 x3
−1.7778 x2 3.5556 x3
===
x1
x2
x3
y i1−2 y i y i−1
h2 p2 y i=0
Ejemplo 27.7
● Se asume {x0}= [ 1 1 1 ]T
● Se reemplaza
3.5556 1 −1.77781 −1.7778 1 3.5556 1 −1.77781 −1.77781 3.5556 1
===
1.77780
1.7778
● Se normaliza el vector de la derecha y se obtiene la primera estimación de λ:
{1.7778
01.7778 }=1.7778 {
101 } ⇒ 1=1.7778
Ejemplo 27.7
● Ahora se reemplaza [1 0 1]T en el lado izquierdo:
3.5556 1 −1.77780 −1.7778 1 3.5556 0 −1.77781 −1.77780 3.55561
===
3.5556−3.55563.5556
● Normalizando,
{3.5556−3.55563.5556 }=3.5556 {
1−11 } ⇒ 2=3.5556
● El error relativo es, ∣ a∣=∣ 3.5556−1.77783.5556 ∣100%=50 %
Ejemplo 27.7● Tercera iteración
[3.5556 −1.7778 0−1.7778 3.5556 −1.7778
0 −1.7778 3.5556 ] {1−11 }= {
5.334−7.1125.334 }=−7.112 {
−0.751
−0.75 }● Cuarta iteración
[3.5556 −1.7778 0−1.7778 3.5556 −1.7778
0 −1.7778 3.5556 ] {−0.75
1−0.75 }= {
−4.4456.223−4.445 }=6.223 {
−0.7141
−0.714 }● Quinta iteración
[3.5556 −1.7778 0−1.7778 3.5556 −1.7778
0 −1.7778 3.5556 ] {−0.714
1−0.714 }= {
−4.3176.095−4.317 }=6.095 {
−0.7081
−0.708 }
Determinación del valor propio menor
● Se aplica el método de potencias a A-1
● El método converge al valor mayor de 1/λ, es decir al menor valor de λ
● Algoritmo en Octave: potencias.m
Problemas 27.1 a 27.29, pag. 822