Apunte EDO Version 2018 PDF
Apunte EDO Version 2018 PDF
Apunte EDO Version 2018 PDF
Prohibida su reproducción
UNIVERSIDAD DE CHILE
FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS
DEPARTAMENTO DE INGENIERÍA MATEMÁTICA
ECUACIONES DIFERENCIALES
ORDINARIAS
Versión 2018
AXEL OSSES
Centro de Modelamiento Matemático
Departamento de Ingenierı́a Matemática
[email protected]
Derecho de autor. Prohibida su reproducción
Se concede permiso para imprimir o almacenar una única copia de este documento.
Salvo por las excepciones más abajo señaladas, este permiso no autoriza fotocopiar o re-
producir copias para otro uso que no sea el personal, o distribuir o dar acceso a copias
electrónicas de este documento sin permiso previo por escrito del Director del Departa-
mento de Ingenierı́a Matemática (DIM) de la Facultad de Ciencias Fı́sicas y Matemáticas
(FCFM) de la Universidad de Chile.
Las excepciones al permiso por escrito del párrafo anterior son: (1) Las copias electróni-
cas disponibles bajo el dominio uchile.cl, (2) Las copias distribuidas por el cuerpo docente
de la FCFM en el ejercicio de las funciones que le son propias.
Este documento fue financiado a través de los recursos asignados por el DIM para la
realización de actividades docentes que le son propias.
Derecho de autor. Prohibida su reproducción
El texto original fue elaborado por el autor en el periodo 2002-2010, fue revisado
por Juan Peypouquet el 2008 quien corrigió y agregó varias secciones y ejercicios
a lo largo del texto y luego en 2010 por el propio autor. Este texto y el material
docente de ejercicios que lo acompaña recibió aportes de los siguientes alumnos de
la Facultad de Ciencias Fı́sicas y Matemáticas de la Universidad de Chile en el
periodo 2002-2017: Francisco Ortega, Oscar Peredo, Andre De Laire, Jorge Lemus,
Nicolás Carreño, Felipe Serrano, Esteban Iglesias, Paulina Arena y Rolando Rogers.
Otras correcciones han sido enviadas por académicos y alumnos que han encontrado
útil este apunte y lo han usado en sus cursos en esta u otras universidades. A todos
ellos el autor les agradece enormemente.
Derecho de autor. Prohibida su reproducción
Derecho de autor. Prohibida su reproducción
Índice general
Capı́tulo 1
1
Derecho de autor. Prohibida su reproducción
2 1. NOCIONES BÁSICAS Y MÉTODOS ELEMENTALES DE RESOLUCIÓN
t=0 t>0
CB
CA
t
Figura 2. Evolución de las concentraciones durante la osmosis
que convergen asintóticamente al promedio M de las concentracio-
nes de ambos medios.
Ley de Osmosis
“El aumento de concentración es proporcional en cada ins-
tante a la diferencia de concentración entre el promedio
asintótico de concentraciones y la concentración actual. La
constante de proporcionalidad cuantifica la permeabilidad
de la membrana.”
La ley de osmosis representada por la ecuación diferencial (2) provee una inter-
pretación más intuitiva y profunda del proceso de osmosis. Más adelante veremos
que (2) tiene como solución (1). La otra ventaja esta ley es que, aparte de su sim-
plicidad, sirve para modelar por analogı́a otros problemas similares, como lo es la
ley de enfriamiento o algunos modelos de población, que veremos más adelante.
Ejercicio Propuesto 1.1. Si en el gráfico del experimento de la Figura 2, el
aumento de la concentración en A hubiese resultado con una forma sigmoide (esto es
estrictamente creciente hacia la ası́ntota pero con un punto de inflexión) ¿Qué mo-
delo diferencial propondrı́a usted para este fenómeno? Indicación: averiguar sobre
el modelo de crecimiento logı́stico.
1.2. Modelamiento de una ley fı́sica: ley de gravitación universal.
Muchos fenómenos de la naturaleza pueden modelarse a través de relaciones entre
cantidades y sus variaciones3, esta idea revolucionaria para la ciencia fue intro-
ducida en el siglo XVII entre otros por Fermat, Newton y Leibniz. En términos
3o fluxiones en la terminologı́a original de Newton, que da la idea de continuo temporal.
Derecho de autor. Prohibida su reproducción
4 1. NOCIONES BÁSICAS Y MÉTODOS ELEMENTALES DE RESOLUCIÓN
Esta ley del siglo XVII lleva a dos ecuaciones con dos derivadas de la forma:
x y
x′′ = − , y ′′ = − 2
(x2 + y 2 )3/2 (x + y 2 )3/2
donde (x, y) es la posición del planeta en el plano de origen en el Sol. Estas ecuacio-
nes tienen como solución elipses en el caso de un solo planeta. En el caso de varios
planetas, en realidad también ejercen fuerza sobre un planeta los demás planetas y
no solamente el Sol, lo que lleva a órbitas muchı́simo más complicadas y al estudio
posterior de órbitas caóticas en el siglo XX por Poincaré entre otros.
El movimiento de los planetas hace parte de los esos fenómenos en los que
se conoce la ley y por lo tanto las ecuaciones diferenciales que representan dicho
fenómeno, pero no se logra comprender completamente la solución a dichas ecua-
ciones. La ley de movimiento en este caso resulta más simple que el movimiento en
si.4
Este ejemplo nos muestra a demás que las ecuaciones diferenciales se pueden
presentar en grupos de dos o más formando sistemas de ecuaciones diferenciales.5
Ejercicio Propuesto 1.2. Averigue quién recopiló los datos a partir de los
que Kepler postuló la idea de órbitas elı́pticas.
1.3. Modelamiento de algunos problemas geométricos. Las ecuacio-
nes diferenciales pueden ser también útiles para plantear y resolver problemas
geométricos. Puede servir por ejemplo para agrupar una familia de curvas. Vea-
mos el caso de las familias ortogonales.
Consideremos la familia F de todas las parábolas con vértice en el origen defi-
nidas por la ecuación y = ax2 , con a ∈ R. Derivando con respecto a x obtenemos
y ′ = 2ax. Si usamos la expresión anterior para eliminar el parámetro a obtenemos
La ecuación diferencial que representa la familia F :
y
y′ = 2 .
x
4Es el caso también de las ecuaciones de Navier-Stokes que modelan el movimiento de un
fluido tridimensional, y no se sabe si dan lugar o no a una única solución; o el de la ecuación de
Schrödinger, que modela la compleja dualidad onda-partı́cula de manera increı́blemente simple.
5Ver Capı́tulo 5.
Derecho de autor. Prohibida su reproducción
2. DEFINICIONES BÁSICAS Y NOCIÓN DE SOLUCIÓN 5
6Si se deriva con respecto a varias variables, se habla de Ecuación Diferencial Parcial (EDP).
Notemos que existen también las inecuaciones diferenciales, donde la igualdad en la Definición 1.1
es reemplazada por una desigualdad.
7Si a (x) = 0 para valores discretos de x, se puede normalizar por intervalos.
n
Derecho de autor. Prohibida su reproducción
3. RESOLUCIÓN DE EDO ELEMENTALES 7
n orden de la EDO
Q(x) = 0 homogénea
Q(x) 6= 0 no homogénea
∀i, ai = cte coeficientes constantes
∃i, ai = ai (x) coeficientes variables
an = 1 normalizada
n
X
Cuadro 1. Clasificación de la EDO lineal ai (x)y (i) = Q(x).
i=0
10Ver Capı́tulo 2 donde se estudia este problema conocido como problema de Cauchy.
11Para ecuaciones de orden n necesitaremos conocer los valores de la función y sus derivadas
hasta el orden n − 1 pues las integraciones sucesivas arrojarán más constantes.
12El Teorema de la Función Implı́cita, que se estudia en el curso de Cálculo en Varias Varia-
bles, da las condiciones para que esto sea posible.
Derecho de autor. Prohibida su reproducción
10 1. NOCIONES BÁSICAS Y MÉTODOS ELEMENTALES DE RESOLUCIÓN
y ′ = xy.
Aquı́ f (x) = x y g(y) = y. Observemos primero que la función nula y(x) ≡ 0 define
una solución de la EDO. Si y 6= 0 hacemos
Z Z
dy
= x dx + C.
y
y ′ = cos2 (y).
En este caso f (x) = 1 y g(y) = cos2 (y). Las soluciones constantes son de la forma
y(x) ≡ π2 + kπ con k ∈ Z. Para encontrar las demás soluciones,
Z Z
dy
= dx + C
cos2 (y)
Z
sec2 (y)dy = x+C
tan(y) = x + C,
y = arctan(x + C).
y = kπ + arctan(x + C),
con C ∈ R y k ∈ Z.
Derecho de autor. Prohibida su reproducción
3. RESOLUCIÓN DE EDO ELEMENTALES 11
3π
2
π
2
− π2
− 3π
2
en la arena para minimizar su tiempo. De modo que su trayectoria debiera ser una
lı́nea quebrada de A a B con vértice en C situado en la interfaz agua-arena. Si a y b
son las distancias verticales de A y B a C y c y x son las distancias horizontales de
A a B y C respectivamente, entonces en tiempo total de la carrera está dado por:
√ p
x2 + a2 (c − x)2 + b2
T (x) = +
v1 v2
Haciendo T ′ (x) = 0 se obtiene la coordenada x (lo que determina el punto C) de
tiempo mı́nimo
x c−x
√ = p
v1 x2 + a2 v2 (c − x)2 + b2
de donde se obtiene a su vez, introduciendo los ángulos de incidencia respecto de
la normal a la interfaz α1 y α2 , la conocida Ley de Snell:
sin α1 sin α2
= .
v1 v2
Si se agrega ahora un tramo extra de cemento a la carrera de para ir de A
a B, primero nadando un tramo de agua desde A, luego atravesando la arena, y
finalmente corriendo en el cemento al punto final B, donde el deportista tiene una
velocidad v3 corriendo en el cemento con v1 < v2 < v3 , entonces de manera similar
es fácil verificar que la trayectoria óptima es una lı́nea con dos quiebres en C1 y C2
en ambas interfaces y se cumple que:
sin α1 sin α2 sin α3
= = = cte.
v1 v2 v3
De manera análoga, la argolla de masa m que cae bajo la gravedad g tiene una
velocidad v que va en aumento de acuerdo a la distancia vertical recorrida y. En
efecto, igualando energı́as cinética y potencial se obtiene:
1 p
mv 2 = mgy ⇒ v = 2gy.
2
Por esta razón, Jean de Bernoulli conjeturó la Ley de Snell generalizada siguiente:
sin α
= cte
v
donde α es el ángulo que instantáneamente forma la tangente a la curva en la posi-
ción de la argolla y la vertical. Usando la propiedad geométrica de que la tangente
de π/2 − α es la derivada de y se obtiene que:
1 1
sin α = cos(π/2 − α) = p 2
=p
1 + tan (π/2 − α) 1 + (y ′ )2
por lo que, usando la Ley de Snell generalizada, se obtiene
1 p
v= p = 2gy.
cte 1 + (y )′ 2
Y 0
−0.2
−0.4
−0.6
−0.8
−1
0 0.5 1 1.5
X
cado por el estudio de la curva cicloide, lo que quedó plasmado en sus obras com-
pletas de 1742 donde aparece su retrato sosteniendo un dibujo de la curva.
Gota de lluvia
“Una gota de lluvia que cae por efecto de su peso va au-
mentando su volumen a medida que captura gotas más pe-
queñas en su superficie inferior. Esto es, a una tasa que
es proporcional a su velocidad de caı́da y a su superficie
inferior.
Supondremos que i) la gota es esférica, ii) la gota alcanza una aceleración constante,
iii) esta aceleración lı́mite es menor que la de gravedad. Si el radio de la esfera es r(t)
entonces su volumen es proporcional a r3 y su superficie media a r2 . Si la densidad
es constante, entonces la masa es m(t) es proporcional a r3 de donde despejando r(t)
resulta proporcional a m1/3 . Con esto, suponiendo que y es la distancia recorrida
verticalmente hacia abajo por la gota, la EDO queda:
m′ (t) = Km2/3 y ′ , K > 0 constante.
Además, la segunda ley de Newton (atención que la masa es variable) es
(my ′ )′ = mg.
Al reemplazar en la primera EDO obtenemos
m′ y ′ + my ′′ = mg
2/3
Km (y ′ )2 + my ′′ = mg
−1/3 ′ 2 ′′
Km (y ) + y = g
K(y ′ )2
m1/3 = .
g − y ′′
Derivando esta expresión vemos que
′
1 −2/3 ′ K(y ′ )2
m m = ,
3 g − y ′′
de donde
′
′ (y ′ )2
y = 3
g − y ′′
2(g − y ′′ )y ′ y ′′ + y ′′′ (y ′ )2
y′ = 3
(g − y ′′ )2
y ′ (g − y ′′ )2 = 6(g − y ′′ )y ′ y ′′ + 3y ′′′ (y ′ )2
3y ′′′ y ′ = (g − y ′′ )(g − y ′′ − 6y ′′ ) = (g − y ′′ )(g − 7y ′′ ).
Suponiendo que y ′′ < g y que la aceleración es constante (y ′′′ = 0) se obtiene que
g
y ′′ = .
7
Derecho de autor. Prohibida su reproducción
3. RESOLUCIÓN DE EDO ELEMENTALES 15
Ahora integrando entre 0 y t una vez (y suponiendo que la gota parte del reposo)
se obtiene la velocidad
gt
y′ = .
7
Reemplazando este valor en la EDO original de la masa se obtiene
gK
m′ = t m2/3 ,
7
que es una EDO a variables separables, que podemos resolver:
Z
gK t2
m−2/3 dm = +C
7 2
gK 2
3m1/3 = t + C.
14
1/3
Si la masa inicial era m0 , entonces C = 3m0 . Luego
3
gK 2 1/3
m(t) = t + m0 .
42
a0 (x) Q(x)
Si a1 (x) 6= 0 podemos normalizar a0 (x) = a1 (x) , Q(x) = a1 (x) y reescribir la
ecuación como
y ′ + a0 (x)y = Q(x).
R
Si multiplicamos ambos lados de la ecuación por el factor integrante µ(x) = exp a0 (x)dx
obtenemos
(µ(x)y(x))′ = µ(x)Q(x)
Z
µ(x)y(x) = µ(x)Q(x)dx + C
Z
C 1
(13) y(x) = + µ(x)Q(x)dx.
µ(x) µ(x)
0 0
CA + CB
Introduciendo la concentración promedio = M , se tiene
2
′
CA + σCA = σM
Z Z Z
′
CA exp σdt + σCA exp σdt = σM exp σdt
′ σt
CA e + σCA eσt = σM eσt
′
CA eσt = σM eσt
Z
CA eσt = σM eσt dt + C
Z
CA = Ce−σt + σM e−σt eσt dt
Z
1 σt
y como eσt dt = e , se tiene que
σ
CA (t) = Ce−σt + M
con C ∈ R. El resultado es una familia de curvas (indexadas por el parámetro C).
Si evaluamos en el tiempo inicial t = 0, se puede encontrar el valor de la constante
C 0 − CB0
C, es decir, CA (0) = C + M y CA (0) = CA 0
. Luego C = CA 0
−M = A . Por
2
lo tanto, la solución es
0 0
CA − CB C 0 + CB 0
CA (t) = e−σt + A
2 2
T = Ce−kt + TA
de donde, evaluando en t = 0, se obtiene C = T0 − TA . Con esto,
T (t) = (T0 − TA )e−kt + TA .
La temperatura del mar tiende exponencialmente a la temperatura ambiente. Más
rápidamente a mayores valores de k. Ver Figura ??.
10
k=3
4
k =2
3
k =1
2
0
0 1 2 3 4 5 6 7 8 9 10
21
20
19
18
17
16
0 2 4 6 8 10 12
ω k
Si consideramos sen φ = √ y cos φ = √ podemos escribir
2
k +ω 2 k + ω2
2
Ak
T (t) = Ce−kt + TA0 + √ (sen(ωt) cos(φ) − cos(ωt) sen(φ)).
k + ω2 |
2 {z }
sen(ωt−φ)
−kt 0 Ak
| {z } + TA + √k 2 + ω 2 sen(ωt − φ). Las variaciones de la
Finalmente, T (t) = Ce
yh | {z }
yp
temperatura del mar se encuentran asintóticamente retrasadas o con desfase posi-
tivo con respecto a las del ambiente (ver Figura 5). Además, notar que la amplitud
Derecho de autor. Prohibida su reproducción
4. ECUACIONES QUE SE REDUCEN A CASOS ELEMENTALES 21
dx dy
= −b cos θ, = b sen θ − a.
dt dt
dy dx dy
Aplicando regla de la cadena tenemos que = . Recordando que cos θ =
dx dt dt
x −y
p y sen θ = p vemos que
2
y +x2 y 2 + x2
−a + b √ −y p
dy b sen θ − a y 2 +x2 −a x2 + y 2 − by
= = = .
dx −b cos θ −bx
−b √ 2x 2
y +x
y
Esta última ecuación es homogénea. Hacemos el cambio de variable z = x ⇒
xz + z = y ′ . Recordando que x e y son positivas tenemos que
′
ap
xz ′ + z = 1 + z2 + z
b
z′ a
√ =
1+z 2 bx
Z
dz a
√ = ln x + C.
1 + z2 b
Reescribimos la constante C como C = ab ln k con k > 0 y obtenemos
p a
ln(z + 1 + z 2 ) = ln(kx) b
p a
z + 1 + z 2 = (kx) b
2a a
1 + z2 = (kx) b − 2z(kx) b + z 2 .
Despejando z y deshaciendo el cambio de variables obtenemos
1 a a
z = (kx) b − (kx)− b
2
y 1 a a
= (kx) b − (kx)− b
x 2
x a a
y = (kx) b − (kx)− b .
2
1
De la condición de borde y(c) = 0 vemos que k = .
c
Ley logı́stica
“El aumento de una población es proporcional al produc-
to entre la población misma y su diferencia respecto a un
valor máximo que es función de los recursos disponibles
limitados.”
z ′ = −M σz + σ,
de donde
Z t Z t Z s
1 1
= exp − M σ(s)ds + exp M σ(s)ds σ(s)ds .
P 0 P0 0 0
Reordenando se obtiene
P0 M
(16) P = .
P0 + (M − P0 ) exp(−M σ t)
Notemos que P (0) = P0 y que P → M si t → ∞.
z′
y1′ − y reemplazando en (17),
z2
2
z′ 1 1
y1′ − 2 = p(x) y1 + + q(x) y1 + + r(x)
z z z
z′ y1 p(x) q(x)
y1′ − 2 = p(x)y12 + 2p(x) + 2 + q(x)y1 + + r(x)
z z z z
z′ y1 p(x) q(x)
y1′ − 2 = [p(x)y12 + q(x)y1 + r(x)] + 2p(x) + 2 +
z z z z
z ′ = −2p(x)y1 z − p(x) − q(x)z,
de donde
z ′ + (2p(x)y1 + q(x))z = −p(x),
que resulta ser una EDO lineal de primer orden no homogénea en la variable z.
1 1
en x = 3 ln(3C), de modo que |y(x)| → ∞ si x → 3 ln(3C).
Ley de Hooke
“Para pequeños desplazamientos en torno a la posición de
equilibrio, la fuerza de restitución del resorte es proporcio-
nal al desplazamiento”.
r
′′ ′′ 2 k
Esto es my = −ky. Es decir, y + w y = 0, con w = . Si z = y ′ entonces
m
dz ′
z′ = y y la ecuación se reescribe como
dy
dz
z + w2 y = 0
dy
dz
z = −w2 y.
dy
z2 y2
= −w2 +C
2 2
(y ′ )2 = −w2 y 2 + 2C
p
′
y = 2C − w2 y 2 .
ρLy ′′ = ρgy
g
y ′′ − y = 0.
L
Derecho de autor. Prohibida su reproducción
28 1. NOCIONES BÁSICAS Y MÉTODOS ELEMENTALES DE RESOLUCIÓN
m
y
r
g
Definiendo σ = se tiene la EDO y ′′ − σ 2 y = 0. El mismo cambio de variable
L
nos da
dz
z − σ2 y = 0
dy
dz
z = σ2 y
dy
z2 y2
= σ2 + C
2 p 2
z = σ 2 y 2 + 2C
p 2C
y′ = σ y 2 + a2 con a= .
σ2
Esto es Z Z
dy
p = σ dt = σt + φ.
y 2 + a2
Haciendo el cambio de variable y = a senh θ en el lado izquierdo,
Z
a cosh θ
dθ = σt + φ
a cosh θ
θ = σt + φ.
Por lo tanto, y = a senh(σt + φ), con φ ∈ R constante.
que es una ecuación en variables separables. Para hallar la solución general hacemos
Z Z
dw
= dx,
w2 + 1
de donde arctan(w) = x + C. Esto es
y = −x + tan(x + C), con C ∈ R.
Capı́tulo 2
1. Definiciones preliminares
Para establecer los teoremas generales de existencia y unicidad de EDO, ası́ co-
mo métodos prácticos de resolución por computador, requerimos primero precisar
la noción de solución de una EDO.
Definición 2.1. Dado un intervalo real I (no reducido a un punto), denota-
remos por C n (I) el conjunto de las funciones y : I → R cuyas derivadas hasta el
orden n existen y son continuas en I y las llamaremos funciones de clase C n .
Es fácil ver que C n (I) es un espacio vectorial con las operaciones usuales de
suma y ponderación por escalar.
Consideremos ahora una función F : I × Rn+1 → R. La noción1 de solución
que estudiaremos es la siguiente:
Definición 2.2. Diremos que una función y : I → R es solución de la ecuación
diferencial
F (x, y, y ′ , y ′′ , . . . , y (n) ) = 0
en el intervalo I si y ∈ C n (I) y
F (x, y(x), y ′ (x), y ′′ (x), . . . , y (n) (x)) = 0
para todo x ∈ I. Por otra parte, la solución general de una ecuación diferencial
es una expresión que agrupa, de manera compacta y explı́cita, la familia de to-
das las soluciones. Una solución cualquiera se denomina también a veces solución
particular.
Ejemplo 2.1. Consideremos la ecuación diferencial y ′ (x) = 0 en un intervalo
I. Recordemos que la derivada de una función en un intervalo es cero si, y sólo si,
la función es constante. Por lo tanto, la solución general de la ecuación es
y(x) = C con C ∈ R.
En muchas ocasiones las ecuaciones diferenciales revelan información cualitativa
acerca de solución, como vemos en el siguiente ejemplo:
Ejemplo 2.2. Consideremos la ecuación diferencial y ′ = 1 − |y| en R. En el
capı́tulo precedente estudiamos las técnicas que nos permiten resolver esta ecua-
ción. Sin embargo, sin necesidad de resolver, solamente al examinar la ecuación
diferencial podemos tener una idea de cómo son las soluciones. En efecto, vemos
1Más adelante, en el capı́tulo de Transformada de Laplace, veremos un caso más general de
soluciones continuas y con derivada continua por pedazos.
31
Derecho de autor. Prohibida su reproducción
32 2. EL PROBLEMA DE CAUCHY: EXISTENCIA, UNICIDAD Y MÉTODOS NUMÉRICOS
2. El problema de Cauchy
En el ejemplo anterior, si bien la ecuación diferencial tiene más de una solución,
el gráfico sugiere que por cada punto (x0 , y0 ) de R2 pasa una y sólo una de estas
curvas. Determinar una (o la) solución de una EDO que pasa por cierto punto dado
es lo que se conoce como resolver el problema de Cauchy.
Definición 2.3. El problema de Cauchy asociado a una ecuación de primer
orden consiste en encontrar y ∈ C 1 (I) tal que
′
y (x) = f (x, y(x)) para todo x ∈ I,
(P C)
y(x0 ) = y0 ,
donde I es un intervalo, x0 ∈ I e y0 ∈ R. La igualdad y(x0 ) = y0 se llama con-
dición inicial o de borde según la variable x se interprete como tiempo o espacio
respectivamente.
Las soluciones del problema de Cauchy se ajustan a las condiciones iniciales o
de borde. Es importante no confundirlas con la familia de soluciones de la ecuación
diferencial sin condiciones iniciales o de borde (solución general).
Dependiendo de f , x0 e y0 , el problema de Cauchy puede o no tener solución
y, en caso de tener, ésta puede o no ser única. Los Teoremas 2.1 y 2.3 más adelante
garantizarán la existencia y unicidad de solución del problema de Cauchy bajo
ciertas condiciones.
Veamos bajo qué hipótesis es razonable esperar que el problema de Cauchy
asociado a una EDO tenga una única solución.
En primer lugar, para que y ′ (x) sea continua es necesario que y(x) y f (x, y)
(como función de dos variables) lo sean. Sin más información sobre la función y
es poco probable que la composición f (x, y(x)) sea continua si f no lo es3. Es de
esperar que necesitemos imponer que la función f sea continua para garantizar
que el problema tenga solución. Sin embargo, esto no es suficiente para obtener la
unicidad, como muestra el siguiente ejemplo:
Ejemplo 2.3. Consideremos el problema de Cauchy:
p
y ′ (x) = y(x) para todo x ∈ R,
y(0) = 0.
La función constante y(x) ≡ 0 es solución, al igual que la función
0 si x < 0,
y(x) = x2
4 si x ≥ 0.
√
Observemos que la derivada de la función g(y) = y diverge cerca del valor
y = 0. Más precisamente,
√ √
y− 0
lı́my→0 y−0 = ∞.
E = C 0 (I) con k · k∞
E = y ∈ C 0 (I)|ψ1 (x) ≤ y ≤ ψ2 (x) (ψ1 , ψ2 continuas en I)
4Cuya versión más general será el Teorema 5.2 que se demostrará más adelante en este texto.
Derecho de autor. Prohibida su reproducción
3. LOS TEOREMAS DE EXISTENCIA Y UNICIDAD 35
Esto haciendo uso del Teorema fundamental del calculo pues f ∈ C 0 (f continua
con respecto a x y como es Lipschitz con respecto a y, luego es continua).
Z x0
(2) y(0) = y0 + f (x, y)dx = y0
x0
Teniendo esto definimos el siguiente operador integral:
φ : C 0 (I) → C 0 (I)
Z x
y 7→ φ(y) = y0 + f (x, y)dx
x0
Podemos ver que φ esta bien definida pues φ(y) ∈ C 0 (I) (Por composición de
funciones).
Derecho de autor. Prohibida su reproducción
36 2. EL PROBLEMA DE CAUCHY: EXISTENCIA, UNICIDAD Y MÉTODOS NUMÉRICOS
5Cuya versión más general será el Teorema 5.4 que se demostrará más adelante también.
Derecho de autor. Prohibida su reproducción
38 2. EL PROBLEMA DE CAUCHY: EXISTENCIA, UNICIDAD Y MÉTODOS NUMÉRICOS
√
Ejemplo 2.8. Retomando el Ejemplo 2.3, notemos que la función f (x, y) = y
no es localmente Lipschitz con respecto a y en un entorno de y0 = 0, de modo que el
Teorema 2.3 no permite esperar que exista una única solución. Sin embargo, se pue-
de probar fácilmente que la función sı́ es localmente Lipschitz si tomamos y0 > 0.
Se deja al lector verificar esto y encontrar el intervalo de existencia y unicidad.
×
×
×
× ×
x0 x1 x2 x3 ··· x0 + L
Solución exacta y su aproximación en el intervalo [x0 , x0 + L].
4.1. Orden del algoritmo. Error local y global. En las siguientes seccio-
nes veremos varios algoritmos de discretización usando la identidad integral (19).
Como ya mencionamos, al discretizar nos interesará que los errores de discretización
globales
en = yn − y(xn )
sean pequeños, esto es, que la solución discretizada se mantenga cerca de la solu-
ción exacta. Sin embargo, el error se va acumulando al pasar de un intervalo Ii al
siguiente y eso quiere decir que los errores en I0 , I1 , . . . , In van a incidir en el error
en al final del intervalo en xn+1 . Por eso llamamos a en error global.
Para saber cuál es el error de discretización local en cada intervalo In , hay que
considerar el problema de Cauchy local
z′ = f (x, z), t ∈ In
z(xn ) = yn
es decir, considerar la solución exacta si ésta coincidiera con la numérica en x = xn .
En general z(xn+1 ) difiere de yn+1 lo que lleva a definir el error de discretización
local:
dn = yn+1 − z(xn+1 ),
que también difiere del error global en .
Definición 2.6. El máximo p tal que
|dn | ≤ Chp+1
para una constante C independiente de h se define como el orden del método.
Que el orden sea p y no p + 1 se explica pues el error de discretización global
eN luego de N = h/L iteraciones será comparable con la suma de los N errores
locales dn , n = 0, . . . , N − 1:
N
X −1
C p
|eN | ≤ dn ≤ CN hp+1 = h
n=0
L
por lo que se pierde un order al sumar y el error acumulado al final del cálculo
medido en el último intervalo se comporta solamente como hp .
No analizaremos el orden de cada uno de los métodos que veremos en detalle,
aunque sı́ lo indicaremos. Es fácil convencerse de que el orden del método depende
principalmente del error de la cuadratura que se escoge para aproximar el término
integral de la derecha en (19). Entonces, cuadraturas por rectángulos, trapecios y
el método de Simpson generarán métodos de cada vez mayor orden.
4.2. Métodos de primer orden. La manera más obvia de aproximar la
integral en (19) es usar un rectángulo. Esto nos da los métodos de Euler de primer
orden, de acuerdo con el extremo del intervalo que usemos para definir la altura del
rectángulo.
y1 = y0 + hf (x0 , y0 )).
yn+1 − yn
= f (xn , yn ) =: dn .
h
dn
xn xn+1
R xn+1
Método de Euler progresivo: xn
f (s, y(s))ds ∼ dn (xn+1 − xn ).
lo que equivale a
yn+1 − yn
= f (xn+1 , yn+1 ).
h
Se dice que este es un método implı́cito pues es necesario resolver la ecuación (21),
que generalmente es no lineal, para determinar yn+1 . Este proceso puede ser difı́cil,
lo cual representa una desventaja de este método. Sin embargo, el método de Euler
retrógrado (y los métodos implı́citos en general) goza de mejores propiedades de
estabilidad, concepto que presentaremos más adelante.
Derecho de autor. Prohibida su reproducción
42 2. EL PROBLEMA DE CAUCHY: EXISTENCIA, UNICIDAD Y MÉTODOS NUMÉRICOS
dn+1
xn xn+1
R xn+1
Método de Euler retrógrado: xn
f (s, y(s))ds ∼ dn+1 (xn+1 − xn ).
cn+1
xn x
bn+1 xn+1
R xn+1
Método de Euler modificado: xn
f (s, y(s))ds ∼ cn+1 (xn+1 − xn ).
Método de Heun. En esta ocasión utilizamos la regla de los trapecios para apro-
ximar la integral:
Z xn+1
hh i
f (s, y(s))ds ∼ f (xn , yn ) + f (xn+1 , y(xn+1 )) .
xn 2
Derecho de autor. Prohibida su reproducción
4. MÉTODOS NUMÉRICOS 43
dn+1
dn
xn xn+1
R xn+1 1
Método de Heun: xn
f (s, y(s))ds ∼ 2
(dn+1 + dn )(xn+1 − xn ).
0 1 1 0 1 0 1 0 1 0
0.05 1.104 1.100 0.004 1.108 0.004 1.104 0.000 1.104 0.000
0.10 1.216 1.208 0.008 1.224 0.009 1.216 0.000 1.216 0.000
0.15 1.336 1.323 0.013 1.350 0.014 1.336 0.000 1.336 0.000
0.20 1.466 1.447 0.018 1.485 0.020 1.465 0.000 1.465 0.000
0.25 1.605 1.581 0.024 1.631 0.026 1.605 0.000 1.605 0.000
0.30 1.755 1.724 0.031 1.788 0.033 1.754 0.001 1.754 0.000
0.35 1.916 1.878 0.038 1.957 0.041 1.915 0.001 1.915 0.001
0.40 2.089 2.043 0.046 2.138 0.050 2.088 0.001 2.088 0.001
0.45 2.274 2.219 0.055 2.333 0.059 2.273 0.001 2.273 0.001
0.50 2.473 2.409 0.064 2.543 0.070 2.472 0.001 2.472 0.001
0.55 2.687 2.612 0.075 2.768 0.082 2.685 0.001 2.685 0.001
0.60 2.915 2.829 0.086 3.010 0.094 2.914 0.002 2.914 0.001
0.65 3.161 3.061 0.099 3.269 0.108 3.159 0.002 3.159 0.001
0.70 3.423 3.310 0.113 3.547 0.124 3.421 0.002 3.422 0.002
0.75 3.705 3.577 0.128 3.845 0.140 3.702 0.002 3.703 0.002
0.80 4.006 3.861 0.145 4.165 0.159 4.003 0.003 4.004 0.002
0.85 4.328 4.166 0.163 4.507 0.178 4.325 0.003 4.326 0.002
0.90 4.673 4.491 0.182 4.873 0.200 4.670 0.003 4.671 0.003
0.95 5.042 4.838 0.204 5.266 0.224 5.038 0.004 5.039 0.003
1 5.437 5.210 0.227 5.686 0.250 5.432 0.004 5.433 0.003
Calculamos los valores de yn que nos entrega el algoritmo y los comparamos con
los valores de la solución exacta del problema de Cauchy. Diremos que un método
de aproximación es
1. Incondicionalmente estable si lı́mn→∞ |yn − y(xn )| = 0 para cualquier valor
de h > 0.
2. Condicionalmente estable si lı́mn→∞ |yn − y(xn )| = 0 para algunos valores
de h > 0.
3. Inestable si |yn − y(xn )| no converge a cero para algún valor de h > 0.
Para entender mejor este concepto analizaremos la estabilidad de los métodos
de primer orden descritos arriba al momento de aproximar ciertas ecuaciones.
con a > 0 e y0 > 0. Se trata de una ecuación lineal de primer orden cuya solución
es y(x) = y0 e−ax . Observemos que lı́mx→∞ y(x) = 0 para todos los valores de a e
y0 . Por lo tanto, lı́mn→∞ y(xn ) = 0.
Derecho de autor. Prohibida su reproducción
4. MÉTODOS NUMÉRICOS 45
y0
× ×
×
× × × × × × × × ×
1 1
h< a h= a
×
1 2
a <h< a ×
× ×
2
×
h≥ a
×
× ×
×
×
×
1 1
yn+1 = yn = y0 .
1 + ah (1 + ah)n+1
Como 1 + ah > 1 para todo h > 0 tenemos que lı́mn→∞ yn = 0 y por lo tanto
lı́mn→∞ |yn − y(xn )| = 0 sin importar el valor de h. Esto nos dice que el método es
incondicionalmente estable.
y0 ×
×
×
×
× × ×
con a > 0 e y0 > 0. Se trata de una ecuación lineal de primer orden cuya solución
es y(x) = y0 eax . Esta vez no es cierto que lı́mx→∞ y(x) = 0 para ningún valor de a.
donde xn+ 21 = xn + h/2 se inspira el siguiente algoritmo que hace parte de los
métodos de Runge-Kutta de orden 4 explı́citos. Este algoritmo es uno de los más
utilizados para resolver ecuaciones diferenciales ordinarias:
Derecho de autor. Prohibida su reproducción
5. EJEMPLO NUMÉRICO: LA LEY DE OMORI EN SISMOLOGÍA 47
y0 = y(0)
g1 = f (xn , yn )
h
g2 = f xn+ 12 , yn + g1
2
h
g3 = f (xn+ 12 , yn + g2 )
2
g4 = f (xn+1 , yn + h g3 )
h
yn+1 = yn + (g1 + 2g2 + 2g3 + g4 )
6
No daremos aquı́ mayores detalles de la deducción del método, para lo cual puede
consultar un libro especializado en métodos numéricos. Sin embrago, para acercarse
a su comprensión haga el siguiente ejercicio.
Ejercicio Propuesto 2.1. Verifique que si f no depende de x se tiene que
g2 = g3 y se recupera en este caso la idea de la cuadratura de Simpson.
1200
1100
1000
900
800
700
Euler
600 Euler Modificado
Heun
Sol Analitica
500
0 20 40 60 80 100 120 140 160 180 200
14 for i=1:N
15 ti=t0+(i−1)*h;
16 y Euler(i+1) = y Euler(i) + h*f(ti);
17 end
18 % 2.− Metodo de Euler modificado
19 y Euler Mod(1)=yM(t0);
20 for i=1:N
21 tim=t0+(i−1)*h+h/2;
22 y Euler Mod(i+1) = y Euler Mod(i) + h*f(tim);
23 end
24 % 3.− Metodo de Heun
25 y Heun(1)=yM(t0);
26 for i=1:N
27 ti=t0+(i−1)*h;
28 tf=t0+i*h;
29 y Heun(i+1) = y Heun(i) + h/2*(f(ti)+f(tf));
30 end
31 figure(2)
32 hold on
33 plot(t,y Euler,'x');
34 plot(t,y Euler Mod,'+');
35 plot(t,y Heun,'o');
36 plot(t,yM(t),'k');
37 legend('Euler','Euler Modificado','Heun','Sol ...
Analitica','Location','SouthEast')
38 title(sprintf('Comparacion de soluciones h= %3.2f',h))
39 hold off
Los datos de las réplicas quedan almacenados como y{1}.data, y{2}.data, y{3}.data,
etc. Las magnitudes mı́nimas respectivas son M(i)=y{i}.magnitud min. Una uni-
dad de tiempo corresponde a 6 horas, y{i}.data(1) es 1 y corresponde al sismo
inicial y y{i}.data(5) es el número de réplicas de magnitud mayor o igual que
M(i) luego de trascurridas (5 − 1) ∗ 6 = 24 horas del terremoto. Si desea ver las
curvas tiene que hacer simplemente:
plot(y{i}.data) o bien
plot(y{i}.fechas,yi.data); datetick(’x’,’dd/mm’)
Derecho de autor. Prohibida su reproducción
5. EJEMPLO NUMÉRICO: LA LEY DE OMORI EN SISMOLOGÍA 51
para verlo con fechas. Si graficamos los datos reales de las réplicas y{1}, y{2},
y{3}, etc. para el sismo del 27 de Febrero del 2010 en función de una nueva variable
τ = ln(ept − 1) − pt para distintos valores de p. Se obtiene aproximadamente un
haz de rectas que se intersectan en un punto, esto es y − y0 = kM (τ − τ0 ) donde
y0 ,τ0 se pueden estimar.
Num acumulado de replicas en funcion del tiempo Grafico de los datos en funcion de tau para p=0.02
1200 1200
Magn>=4.0 (y1)
Magn>=4.2 (y2)
1000 1000
Magn>=4.4 (y3)
800 800
numero de replicas
Magn>=4.6 (y4)
600 600
Magn>=4.8 (y5)
400 400
Magn>=5.0 (y6)
200 200
Magn>=5.2 (y7)
Magn>=5.4 (y8)
Magn>=5.6
Magn>=5.8 (y9)
(y10)
Magn>=6.0 (y11)
0 0
21/02 28/02 07/03 14/03 21/03 28/03 04/04 11/04 18/04 −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0
tiempo
y{i} M kM y0
1 4.0 290 396
2 4.2 258 367
3 4.4 211 322
4 4.6 150 254
5 4.8 93 190
6 5.0 51 139
7 5.2 19 89
8 5.4 9 45
9 5.6 6 25
10 5.8 4 15
1000
800
600
400
200
−200
0 20 40 60 80 100 120 140 160 180 200
5.4. Modificación del modelo. Observe del gráfico anterior que el dı́a 11
de marzo se produce una discontinuidad en la derivada de las curvas y M . Podemos
mejorar el modelo aproximando el haz de rectas por una función lineal por pedazos,
esto es kM es una constante diferente antes y después del 11 de marzo. Al hacer esto,
se obtienen ajustes como el de la figura, donde en azul se grafica la solución numérica
usando el método de Euler modificado y en rojo los datos reales: Notemos que
Ajuste mejorado
1200
1000
800
600
400
200
−200
0 20 40 60 80 100 120 140 160 180 200
aunque todavı́a es posible obtener una solución analı́tica en este caso, si kM = kM (t)
fuera una función complicada del tiempo, de modo tal que kM (t) eptp−1 no tenga
una primitiva analı́tica, entonces no podrı́a obtenerse una solución anal´?tica y el
único recurso serı́an los métodos numéricos. Esto suele suceder en la práctica en
la ingenierı́a y las ciencias, en que los modelos simples con que se introducen y se
estudian los fenómenos tienen solución analı́tica, pero para aplicarlos en la práctica
Derecho de autor. Prohibida su reproducción
es necesario introducir modelos más complejos que suelen no tener dichas soluciones
y es necesario recurrir a las aproximaciones numéricas.
5.5. Capacidad predictiva. ¿Serı́a posible estimar con nuestro mejor mo-
delo el número de réplicas con magnitud mayor o igual a 5 que se producirı́an del
16 de abril al 16 de mayo del 2010? El modelo puede dar una predicción. Tomando
M = 5 y calculando yM (16/05/2010) − yM (16/04/2010) con el modelo de sección
anterior se obtienen 2 réplicas de magnitud mayor o igual a 5.
44 t=1:length(y{i}.data);
45 tau2=log(exp(p2*t)−1)−p2*t;
46 kM2(i)=(y{i}.data(tf2)−y{i}.data(ti2))/(tau2(tf2)−tau2(ti2));
47 y02(i)=y{i}.data(ti2);
48 end
49 tau02=tau(ti2);
50 % Grafico
51 figure(6)
52 hold on
53 for i=1:10
54 t=1:length(y{i}.data);
55 tau=log(exp(p*t)−1)−p*t;
56 tau2=log(exp(p2*t)−1)−p2*t;
57 y analitica=y0(i)+kM(i)*(tau−tau0);
58 y analitica2=y02(i)+kM2(i)*(tau−tau02);
59 y analitica3=y analitica. *(t<ti2)+y analitica2.*(t>=ti2);
60 plot(t,y{i}.data,'r.',t,y analitica3,'−');
61 end
62 title('Ajuste mejorado')
63 hold off
8Esta sección puede leerse también junto con el Capı́tulo 5, aunque no se requiere.
Derecho de autor. Prohibida su reproducción
6. EJEMPLO NUMÉRICO: LOTKA-VOLTERRA Y LA PESCA EN EL MAR ADRIÁTICO 55
y(t)
III II
IV I
x(t)
(x0 , y0 )
= (x(0), y(0))
g1= f (xn , yn )
h
g2 = f (xn , yn ) + g1
2
h
g3 = f (xn , yn ) + g2
2
g4 = f ((xn , yn ) + h g3 )
h
(xn+1 , yn+1 ) = (xn , yn ) + (g1 + 2g2 + 2g3 + g4 )
6
De hecho, es posible verificar que un método de tipo Euler de orden 1 o de
Runge-Kutta de orden 2 no provee suficiente precisión para que las órbitas sean
cerradas, y es por esto que es necesario utilizar un método de orden 4.
Los resultados numéricos obtenidos con el método de Runge-Kutta de orden
4 descrito más arriba se muestran en la Figura 3, donde se hizo la simulación del
paso de una situación con pesca (lı́neas oscuras) a una situación sin pesca (lı́neas
claras)9. Observe que el punto de equilibrio que está en el centro de las órbitas se
desplaza arriba y hacia la izquierda como predice la teorı́a.
10B. van der Pol and J. van der Mark, The heartbeat considered as a relaxation oscillation,
and an electrical model of the heart. The London, Edinburgh, and Dublin Philosophical Magazine
and Journal of Science Ser.7, 6, 763-775, 1928
11There?s No Quiet Without Noise, by David Colman, New York Times, July 8, 2011.
Derecho de autor. Prohibida su reproducción
1 % Funcion que implementa el sistema de EDO de Van der Pol Caso No ...
Forzado.
2 % dX = F(t,X,mu) entrega el vector de dimension 2 correspondiente ...
a la
3 % derivada de X dada por la EDO de Van der Pol con el parametro ...
mu, caso
4 % no forzado, ie, f = 0.
5
6 function dX = F(t,X,mu)
7 x=X(1);
8 y=X(2);
9 dX=[y;mu*(1−xˆ2)*y−x];
10 end
y0 = y(0)
g1 = f (xn , yn )
h
g2 = f (xn+ 21 , yn + g1 )
2
h
g3 = f (xn+ 21 , yn + g2 )
2
g4 = f (xn+1 , yn + hg3 )
h
yn+1 = yn + (g1 + 2g2 + 2g3 + g4 )
6
En este caso:
x ← t
y ← X = (x, y)de la EDO
f ← dX = F (t, X, µ)
I(2) − I(1)
h ←
N
Figura 4. Solución Van der Pol para µ = 0,1, (x0 , y0 ) = (3, 1).
Derecho de autor. Prohibida su reproducción
62 2. EL PROBLEMA DE CAUCHY: EXISTENCIA, UNICIDAD Y MÉTODOS NUMÉRICOS
1 % Funcion que implementa el sistema de EDO de Van der Pol Caso Forzado.
2 % dX = F2(t,X,mu,A,omega) entrega el vector de dimension 3 que
3 % corresponde a la derivada de X del sistema de EDO de Vander Pol ...
con el
4 % parametro mu y una funcion f(t) = A*sin(omega*t).
5
6 function dX = F2(t,X,mu,A,omega)
7 x=X(1);
8 y=X(2);
9 z=X(3);
10 dX=[y;mu*(1−xˆ2)*y−x+A*cos(z);omega];
11 end
22 subplot(3,2,p)
23 plot(t, x); % grafica (t,x)
24 xlabel('t'); ylabel('solucion')
25 title(['van der Pol, mu = ',num2str(mu)]);
26 subplot(3,2,p+1)
27 plot(x, y); % grafica (x,y)
28 title(['van der Pol, mu = ',num2str(mu)]);
29 xlabel('x'); ylabel('y = dx')
30
31 p = p+2; % se pasa a la siguiente fila del subplot
32 soundsc(x) % con esto se 'escucha' la solucion.
33
34 end
35
36 end
el que se escucha un sonido más irregular, como un ruido. Para µ = 6, cuya solución
se escucha como un sonido regular, el ciclo lı́mite se observa con gran nitidez, al
contrario del caso µ = 8,53, cuya solución se escucha como un ruido irregular, donde
se tiene que la gráfica de (x,y) está delimitada, pero es difı́cil determinar a partir
de la gráfica si existe un ciclo lı́mite o no. De lo anterior se deduce que el orden
(solución periódica) corresponde al caso µ = 6 y el caos, al caso µ = 8,53.
Ejercicio Propuesto 2.2. El siguiente es un modelo epidemiológico intro-
ducido por Kermack & McKendrick, ambos británicos, en 192712. Modela la com-
posición de una población constante N = x + y + z ante una enfermedad, donde
x(t) representa el número de personas suceptibles de contagiarse, y(t) el número de
infectados y z(t) el número de inmunes a la enfermedad para cada t ≥ 0:
x′ = −β xy
y′ = β xy − γ y
z′ = γy
x(0) = x0 > 0, y(0) = y0 > 0, z(0) = 0, x0 + y0 = N,
donde β > 0 es un parámetro de infección y γ > 0 la tasa de inmunización. Todo esto
suponiendo que las soluciones x, y, z existen, son únicas, no negativas, continuas y
con derivada continua para t ≥ 0. Pruebe los siguientes hechos a partir del modelo:
La epidemia crece inicialmente en el número de infectados si el número
inicial de suceptibles supera el umbral ρ = βγ , esto es x0 > ρ ⇒ y ′ (0) > 0.
Suponiendo que existe el lı́mite
lı́m (x, y, z) = (x∞ , y∞ , z∞ )
t→∞
se puede demostrar que y∞ = 0, x∞ = N − z∞ y z∞ es raı́z de la función:
−z
F (z) = N − z − x0 exp .
ρ
Esto puede obtenerse de dividir la primera y última de las EDO.
Se puede usar el Teorema del Valor intermedio para probar que la raı́z
anterior existe y es única.
Discretice estas ecuaciones para valores de los parámetros que elija, tal como se
hizo en el ejemplo anterior (notar que la variable dependiente es t no x como
antes, además hay tres componentes y no dos como antes). Tomando un número
de individuos suceptibles inferior y superior al umbral teórico ρ y simule.
Capı́tulo 3
67
Derecho de autor. Prohibida su reproducción
68 3. EDO LINEALES DE ORDEN SUPERIOR
n
X
P (x, D) = ak (x)Dk
k=0
= an (x)Dn + an−1 (x)Dn−1 + . . . + a1 (x) D + a0 (x),
P (x, D)f (x) = an (x)f (n) (x) + · · · + a1 (x)f ′ (x) + a0 (x)f (x).
Si los coeficientes ak son todos constantes se dice que el operador diferencial lineal
es a coeficientes constantes
El operador diferencial nos servira para definir las EDOs lineales de orden n de
mejor forma.
P (D)y = Q(x),
Pn
donde P (x) = k=0 ak Dk , ak ∈ R, an = 1
Para esto nos resultara útil mostrar algunas propiedades de este operador.
Pn
Propiedades 3.1. (Factorización) Sea P (D) = k=0 ak Dk operador diferen-
cial de orden n, luego puede ser escrito como:
l
Y
P (D)y = (D − λj )mj y
j=1
P (D)y = D2 y + a1 Dy + a0 y = (D − λ1 )(D − λ2 )y = Q,
donde λ1 y λ2 son las raı́ces del polinomio caracterı́stico
p(λ) = λ2 + a1 λ + a0 .
Introduciendo la variable auxiliar z = (D − λ2 )y se obtiene el siguiente sistema de
dos EDO lineales de primer orden que sabemos resolver:
(D − λ1 )z = Q
(D − λ2 )y = z.
De la primera ecuación encontramos z, que sustituimos en la segunda ecuación
para obtener y. Más precisamente, los resultados de la Sección 3 nos dan
Z
z = C1 eλ1 x + eλ1 x e−λ1 x Q(x)dx
Z
λ2 x λ2 x
y = C2 e +e e−λ2 x z(x)dx,
de donde
Z
y(x) = C2 eλ2 x + C1 eλ2 x e(λ1 −λ2 )x dx
| {z }
Solución Homogénea yh (x)
Z Z
+ eλ2 x e(λ1 −λ2 )x e−λ1 t Q(t)dt dx .
| {z }
Solución Particular yp (x)
3.1. Solución homogénea. Se tienen tres casos para los valores caracterı́sti-
cos:
k +
F(t) q
E(t) -
C
m
b y R
2. Si ∆ = 0 las raı́ces
√ son reales e iguales. Decimos que el sistema está crı́ticamente
amortiguado (b = 2 km. En este caso la solución es
yh = C1 e−Bt + C2 te−Bt .
Derecho de autor. Prohibida su reproducción
3. ESTUDIO COMPLETO DE LA EDO DE ORDEN DOS 73
3. Finalmente, si ∆ < 0 las √raı́ces son complejos conjugados. Decimos que el sistema
está subamortiguado (b < 2 km) y la solución es
√ √
yh = C1 e−Bt sen −∆ t + C2 e−Bt cos −∆ t .
√
De manera equivalente podemos escribir yh = C3 e−Bt sen −∆ t + C4 , pues
sen(a + b) = sen a cos b + sen b cos a.
20
15
10
5
0 10 20 30 40 50 60 70 80 90 100
30
20
10
0
0 10 20 30 40 50 60 70 80 90 100
10
−10
0 1 2 3 4 5 6 7 8 9 10
3.2. Condiciones de borde. Hemos visto que una ecuación lineal de se-
gundo orden tiene una única solución si especificamos el valor de la función y su
derivada en un punto. Otro problema interesante es determinar si existe alguna
solución que parta de un punto y llegue a otro predefinido. Este tipo de problema
aparece con frecuencia en fı́sica, ingenierı́a y economı́a. Consideremos, por ejemplo,
un problema del tipo ′′
y + λy = 0
y(0) = 0
y(1) = 0.
Deseamos saber para cuáles valores de λ este problema tiene soluciones y cuántas.
De acuerdo con el signo de λ podemos distinguir 3 casos:
Derecho de autor. Prohibida su reproducción
74 3. EDO LINEALES DE ORDEN SUPERIOR
Y
p
S
Y
1
Y
2
Gracias al Lema 3.1, podemos expandir la última fila. Todos los sumandos,
salvo el último, tienen una fila repetida, por lo que la derivada del Wronskiano
queda de la forma
(27) W ′ (x) = −an−1 (x)W (x).
Al integrar la ecuación (27) obtenemos la Fórmula de Abel:
Derecho de autor. Prohibida su reproducción
80 3. EDO LINEALES DE ORDEN SUPERIOR
Resumiremos algunas propiedades del Wronskiano que tienen que ver con in-
dependencia lineal de funciones en el siguiente resultado:
Note que esta factorización permite de hecho resolver la EDO de orden n como
una serie de EDO de orden 1, En efecto, es claro que
P (D)y = 0 ⇔ (D − λ1 )m1 ◦ . . . (D − λl )ml y = 0,
Derecho de autor. Prohibida su reproducción
82 3. EDO LINEALES DE ORDEN SUPERIOR
Lo que quiere decir que eλ1 x , eλ2 x , . . . , eλl x son soluciones de (H). Veamos que son
l.i.:
eλ1 x ... eλl x
λ1 eλ1 x . . . λl eλl x
W (x, eλ1 x , . . . , eλl x ) = .. .. ..
. . .
λl−1 eλ1 x . . . λl−1 eλl x
1 l
1 ... 1
!
X l λ1 . . . λl
= exp λi x .. .. ..
. . .
i=1
λl−1 . . . λl−1
1 l
l
!
X Y
= exp λi x C (λi − λj )
i=1 i6=j
Derecho de autor. Prohibida su reproducción
6. MÉTODOS DE RESOLUCIÓN DE EDO DE ORDEN n A COEFICIENTES CONSTANTES83
Ejemplo 3.3.
(D − 1)2 (D + 3)5 D3 y = 0
En este caso n = 10. Veamos todas las soluciones que se encuentran con este método:
λ1 = 1, m1 = 2 ⇒ {ex , xex }
λ2 = −3, m2 = 5 ⇒ {e−3x , xe−3x , x2 e−3x , x3 e−3x , x4 e−3x }
λ3 = 0, m3 = 3 ⇒ {1, x, x2 }
λi x
e αi,mi −1 (mi − 1)!p̃(λi ) = 0
Se verifica que p̃(λi ) 6= 0 y (mi − 1)! 6= 0, por lo tanto αi,mi −1 = 0. Repitiendo la
aplicacion progresivamente, se verifica que todos los αi,mi −1 son iguales a cero, con
lo cual las n funciones son l.i..
Antes de pasar a los ejemplos resumiremos la discusión anterior en el siguiente
resultado:
Teorema 3.7. Supongamos que λ1 , . . . , λl son los valores caracterı́sticos, con
multiplicidades m1 , . . . , ml , de la ecuación homogénea
(D − λ1 )m1 . . . (D − λl )ml y = 0.
Una base para H es la generada por las siguientes funciones linealmente indepen-
dientes:
Para cada valor caracterı́stico λ real
eλx , xeλx , ..., xm−1 eλx .
Para cada par de valores caracterı́sticos complejos α ± iβ,
eαx cos(βx), xeαx cos(βx), ... xm−1 eαx cos(βx),
eαx sen(βx), xeαx sen(βx), ... xm−1 eαx sen(βx).
Ejemplo 3.4.
(D − 1)2 (D + 3)5 D3 y = 0.
En este caso n = 10. Veamos todas las soluciones que se encuentran con este método:
λ1 = 1, m1 = 2 ⇒ {ex , xex },
λ2 = −3, m2 = 5 ⇒ {e−3x , xe−3x , x2 e−3x , x3 e−3x , x4 e−3x },
λ3 = 0, m3 = 3 ⇒ {1, x, x2 }.
con λ0 ∈ C y qk0 eλ0 x una funcion que puede ser un polinomio por una exponencial,
un polinomio por seno o coseno por una exponencial o cualquier combinación lineal
de estos, donde gr(qk0 ) ≤ k0 en caso de ser un polinomio.
Este método tiene la desventaja aparente de ser muy restrictivo en cuanto al
tipo de lado derecho para el cual es útil. Sin embargo, estas funciones son sumamente
frecuentes en la aplicaciones.
El operador diferencial que anula el lado derecho es (D − λ0 )k0 +1 , pues
Z Z
1
y2 = Cy1 exp − a1 (x)dx du.
y12
Ejemplo 3.9. Se tiene la EDO x2 y ′′ + xy ′ − y = 0 y una solución y1 = x. Se
pide encontrar otra solución linealmente independiente y2 . Para ello escribimos la
ecuación normalizada:
1 1
y ′′ + y ′ − 2 y = 0.
x x
La fórmula de Abel nos da
Z
dx
xy2′ − y2 = C exp − ,
x
de donde
1 ′ 1 C
y2 − 2 y2 =
x x x3
y ′ C
2
=
x x3
y2 C
= + D.
x −2x2
Tomamos C = −2, D = 0 y obtenemos y2 = x1 .
5Se puede usar la regla de Cramer en lugar de calcular la inversa. Esto da origen a la fórmula
de representación de Green, que estudiaremos más adelante.
Derecho de autor. Prohibida su reproducción
92 3. EDO LINEALES DE ORDEN SUPERIOR
Capı́tulo 4
Transformada de Laplace
1. Definiciones y ejemplos
Definición 4.1. Dada f : [0, +∞) → R se llama transformada de Laplace de
f a la función
Z +∞
(28) L[f ](s) = e−st f (t)dt
0
Ejemplo 4.3. Recordemos ahora que cos wt + i sen wt = eiwt , de manera que
Re(eiwt ) = cos wt e Im(eiwt ) = sen wt son la parte real e imaginaria de eiwt , res-
pectivamente. De ejemplo anterior tenemos que
1 s + iw
L[eiwt ](s) = = 2 para s > 0.
s − iw s + w2
95
Derecho de autor. Prohibida su reproducción
96 4. TRANSFORMADA DE LAPLACE
0 1 2 3 4 5 6
−1
0 1 2 3 4 5 6
Propiedad 4.2. Si f ∈ Cα entonces para todo s > α, existe L[f ](s) (y converge
absolutamente). Además
M
|L[f ](s)| ≤
s−α
para todo s > α. En particular, lı́ms→+∞ L[f (t)](s) = 0.
Demostración. Recordemos que si la integral converge absolutamente, en-
tonces converge. Con esto basta probar que
Z ∞
M
|e−st f (t)|dt ≤
0 s −α
para todo s > α. En efecto,
Z ∞ Z ∞
|e−st f (t)|dt = e−st |f (t)|dt
0 0
Z ∞
≤ M e−st eαt dt
0
Z ∞
≤ M e−(s−α)t dt
0
∞
1
(−s+α)t
≤ M e
−s + α 0
C
≤ .
s−α
Ejercicio Propuesto 4.1. Demostrar que tk , k ∈ N, tiene transformada de
k!
Laplace y que L[tk ](s) = k+1 , s > 0.
s
La funcı́on escalón de Heaviside se define como
0 t<a
Ha (t) =
1 t≥a
a
Función escalón de Heaviside Ha (t).
Derecho de autor. Prohibida su reproducción
2. PROPIEDADES BÁSICAS DE LA TRANSFORMADA DE LAPLACE 99
a b
Pulso entre a y b.
Notar que Pab (t) = Ha (t) − Hb (t). Luego, su transformada de Laplace para
0 ≤ a < b será
L[Pab (t)](s) = L[Ha (t) − Hb (t)](s) = L[Ha (t)](s) − L[Hb (t)](s).
Por lo tanto
1 −as
L[Pab (t)](s) = (e − e−bs ) para s > 0.
s
Z +∞
L[f ′ ](s) = e−st f ′ (t)dt
0
Z ∞ ∞
= s e−st f (t)dt + e−st f (t)
0 0
= sL[f ](s) + lı́m e−st f (t) − lı́m e−st f (t).
t→∞ t→0+
para s > α.
Derecho de autor. Prohibida su reproducción
2. PROPIEDADES BÁSICAS DE LA TRANSFORMADA DE LAPLACE 101
Z tZ t Z t Z u
Si denotamos f (u)du = f (v)dv du. La transformada de esta
a a a a
expresión es
Z tZ t Z t Z Z
1 1 a t
L f (u)du (s) = L f (u)du (s) − f (u)du
a a s a s 0 a
Z a Z a
1 1 1
= L[f ](s) − 2 f (u)du − F (t)dt.
s2 s 0 s 0
repitiendo el proceso, se obtiene la fórmula para la transformada de la n-ésima
integral:
Z t Z t Xn Z a Z t Z t
L . . . f (u)du (s) = 1 L[f ](s) − 1
. . . f (u)du.
sn sk 0 a
a a
| {z } k=1 | {z a}
n veces n − k veces
Z ∞ Z t
L[(f ∗ g)](s) = e−st f (u)g(t − u)du dt
0 0
Z ∞ Z t
−st
= e f (u)g(t − u)dudt.
0 0
Z ∞ Z t Z ∞ Z ∞
e−st f (u)g(t − u)dudt = e−st f (u)g(t − u)dtdu
0 0
Z0 ∞ Zu∞
= e−s(w+u) f (u)g(w)dwdu
0 0
Z ∞ Z ∞
= e−su f (u)du e−sw g(w)dw
0 0
= L[f ](s)L[g](s).
Luego notemos que ([kT, (k + 1)T ])k∈N es una partición disjunta de [0, ∞), luego:
Z ∞ ∞ Z
X (k+1)T
e−su f (t)du = e−su f (u)du u = kT + t du = dt
0 k=0 kT
X∞ Z T
= e−s(kT +t) f (kT + t)du
k=0 0
X∞ Z T
= e−s(kT +t) f (kT + t)du
k=0 0
X∞ Z T
= e−skT e−st f (t)dt
k=0 0
∞
!Z
X T
−sT k
= (e ) e−st f (t)dt
k=0 0
P
Luego como |e−sT | ≤ 1 ∀s ≥ 0, entonces ∞ k=0 (e ) = 1−e1−sT (serie geométrica)
−sT k
Concluyendo lo pedido.
C −(s−α)n
kΦn − ΦkI ≤ sup e
s>α s − α
C
= e−(α0 −α)n .
α0 − α
Por lo tanto kΦ − ΦnkI → 0 cuando n → ∞, con lo que Φn converge uniformemente
a Φ en I.
Derecho de autor. Prohibida su reproducción
104 4. TRANSFORMADA DE LAPLACE
p(x) A1 An
= + ··· + .
q(x) x − a1 x − an
2. q(x) tiene una raı́z real de multiplicidad n: q(x) = (x − a)n . Hacemos
p(x) A1 A2 An
= + + ··· + .
q(x) x − a (x − a)2 (x − a)n
3. q(x) tiene 1 raı́z compleja simples. Entonces su conjugado también es raı́z simple,
de modo que q(x) = ax2 + bx + c = (x − z)(x − z), z ∈ C \ R. Hacemos
p(x) Ax + B
= 2 .
q(x) ax + bx + c
4. Si q(x) tiene 1 raı́z compleja de multiplicidad n: q(x) = (ax2 + bx + c)n =
(x − z)n (x − z)n con z ∈ C \ R. Hacemos
p(x) A1 x + B1 A2 x + B2 An x + Bn
= 2 + + ··· + .
q(x) ax + bx + c (ax2 + bx + c)2 (ax2 + bx + c)n
5. Cualquier combinación de estos casos.
Escribimos entonces
R(s) L[Q(t)](s)
L[y(t)](s) = + ,
P (s) P (s)
donde
n
X n
X
P (s) = aj sj y R(s) = aj Rj (s).
j=1 j=1
1
= L[H(t)](s − λ) = L[eλt H(t)](s)
s−λ k−1
1 t tk−1
= L (s − λ) = L eλt (s)
(s − λ)k (k − 1)! σt (k − 1)!
1 1 e
= L sen(wt) (s − σ) = L sen(wt) (s)
(s − σ)2 + w2 w w
s−σ
= L [cos(wt)] (s − σ) = L [eσt cos(wt)] (s).
(s − σ)2 + w2
−1 s −1 1
1. L (t) y 2. L (t)
(s2 + a2 )n (s2 + a2 )n
Z t
1
2. Recordando que L f = L[f ], se tiene:
0 s
1 −1 1 s
L−1 (t) = L (t)
(s2 + a2 )2 s (s2 + a2 )2
1 1
= L−1 L t sen(at) (s) (t)
s 2a
Z t
−1 1
= L L t sen(at) (s) (t)
0 2a
Z t
1
= t sen(at)
0 2a
1 1
= sen(at) − t cos(at) .
2a2 a
Resumimos los ejemplos más relevantes de funciones y sus transformadas en la
siguiente tabla:
f (t) L[f ](s)
1
eλt
s−λ
1 −as
Ha (t) e
s
1 −as
Pab (t) (e − e−bs )
s
tk−1 1
eλt
(k − 1)! (s − λ)k
eσt 1
sen(wt)
w (s − σ)2 + w2
s−σ
eσt cos(wt)
(s − σ)2 + w2
2as
t sen(at)
(s2 + a2 )2
1 2a2
sen(at) − t cos(at)
a (s2 + a2 )2
Ejemplo 4.12. Continuando el Ejemplo 4.9 del cohete, hay que encontrar las
funciones que al ser transformadas por L, equivalen a cada uno de los sumandos
del miembro derecho de la ecuación
v0 a a g
L[y](s) = 2 + 3 e−t1 t − 3 e−t2 t − 3 .
s s s s
Utilizando las conclusiones de la discusión precedente vemos que
a 2 a 2 gt2
L[y](s) = L v0 t + (t − t1 ) H(t − t1 ) − (t − t2 ) H(t − t2 ) + .
2 2 2
Derecho de autor. Prohibida su reproducción
110 4. TRANSFORMADA DE LAPLACE
s+1 1
Como L [e−x cos x] (s) = 2
y L [e−x sen x] (s) = , se tiene
(s + 1) + 1 (s + 1)2 + 1
que
L[y](s) = L e−x cos x (s) + 2L e−x sen x (s)
−x
L[y](s) = L e cos x + 2e−x sen x (s)
y(x) = e−x (cos x + 2 sen x).
cada función continua f le asigna un número real, que es su valor en cero f (0).
Las funciones generalizadas se definen en términos de su acción sobre cierto tipo
de funciones. En estricto rigor, definimos la acción de la delta de Dirac sobre una
función continua f mediante
Z ∞
δ[f ] = δ(t)f (t) dt = f (0).
−∞
La notación de integral es bastante sugerente pero se debe recordar que la delta
de Dirac no es una función, aunque comparta ciertas propiedades con las funciones
integrables. Siguiendo la notación presentada arriba, escribiremos
Z ∞
δa [f ] = δ(t − a)f (t) dt = f (a).
−∞
La delta de Dirac modela un impulso, golpe o chispazo. Se interpreta como un
gran estı́mulo aplicado instantáneamente a un sitema.
12
10
para s ∈ R.
Capı́tulo 5
1. Introducción
Hasta aquı́ se han estudiado ecuaciones lineales con una sola función incógnita.
Estudiaremos ahora sistemas de ecuaciones diferenciales lineales, tanto porque per-
miten modelar problemas fı́sicos que involucran la interacción de diversas variables
dependientes, como también porque a ellos pueden reducirse las ecuaciones lineales
de grado superior.
Lema 5.1. Si A(t) ∈ Mn×m (R) y B(t) ∈ Mm×p (R), con n, m, p ∈ N, entonces
113
Derecho de autor. Prohibida su reproducción
114 5. SISTEMAS LINEALES DE PRIMER ORDEN
bλ
b V V b
b (1+λ)
Consideremos el sistema
(31) x′1 = a11 x1 + a12 x2 + b1
(32) x′2 = a21 x1 + a22 x2 + b2
con las condiciones iniciales x1 (0) = x01 y x2 (0) = x02 .
2.2.1. Método de sustitución. Suponiendo a21 6= 0, despejamos x1 de (32)
para obtener
x′2 − a22 x2 − b2
(33) x1 =
a21
y luego derivamos
x′′2 − a22 x′2 − b′2
(34) x′1 = .
a21
Reemplazando (33) y (34) en (31) resulta
x′′2 − (a11 + a22 )x′2 + (a11 a22 − a12 a21 )x2 = a21 b1 − a11 b2 + b′2 ,
que es una ecuación de segundo orden con condiciones iniciales
x2 (0) = x02
x′2 (0) = a21 x01 + a22 x02 + b2 (0).
2.2.2. Método de reducción. Utilizando la notación de operador diferencial
Dx = x′ el sistema (31)-(32) se reescribe
(35) (D − a11 )x1 − a12 x2 = b1
(36) −a21 x1 + (D − a22 )x2 = b2 .
Aplicando el operador D − a11 a la segunda ecuación, multiplicando por a21 la
primera y sumando, se obtiene
(37) (D − a11 )(D − a22 )x2 − a12 a21 x2 = a21 b1 + (D − a11 )b2 ,
que es la misma EDO de orden 2 para x2 que la obtenida con el método anterior.
Observación: Puede verificar que en todos los casos, despejando x1 o x2 , se
llega a una EDO del tipo
x′′i − (a11 + a22 )x′i + (a11 a22 − a12 a21 )xi = fi , i = 1, 2
| {z } | {z }
tr(A) det(A)
4. Si A o B son continuas por pedazos, el Teorema 5.1 sigue siendo cierto, salvo
porque X ′ es sólo continua por pedazos (se suele decir que X es C 1 por pedazos).
Observe que coincide con la definición usual de norma euclideana para vectores
(pensados aquı́ como matriz fila o columna).
Derecho de autor. Prohibida su reproducción
118 5. SISTEMAS LINEALES DE PRIMER ORDEN
Observemos primero que la función
F (t, X) = A(t)X + B(t)
es Lipschitz con respecto a su segunda variable. Es decir, existe L > 0 tal que
kF (t, X) − F (t, Y )k ≤ LkX − Y k
para todo t ∈ I y para todos los X, Y ∈ Rn . En efecto,
|F (t, X) − F (t, Y )| ≤ kA(t)(X − Y )k
1/2
Xn X n
≤ a2ij (t) kX − Y k.
i=1 j=1
Para la demostración del Teorema 5.2 primero observemos que, integrando entre
t0 y t la ecuación diferencial, obtenemos
Z t
X(t) = X0 + F (s, X(s))ds.
t0
Definimos el operador
T : C 0 (I)m → C 0 (I)m
X 7 → T (X),
donde T (X) es la función definida por
Z t
(38) T (X)(t) = X0 + F (s, X(s))ds para t ∈ I.
t0
Derecho de autor. Prohibida su reproducción
3. LOS TEOREMAS DE EXISTENCIA Y UNICIDAD. DEMOSTRACIÓN. 119
para todo t ∈ I. Además, todo punto fijo tiene que ser de clase C 1 , por ser primitiva
de una función continua. Al derivar tenemos que todo punto fijo de T es solución
del problema de Cauchy
X ′ (t) = F (t, X(t)) para todo t ∈ I
X(t0 ) = X0 .
Recı́procamente, toda solución del problema de Cauchy debe ser un punto fijo de
T . Por lo tanto, para demostrar el Teorema 5.2 basta verificar que T tiene un único
punto fijo.
1. Una sucesión {xn } es de Cauchy si para todo ε > 0 existe N ∈ N tal que
|xn − xk | < ε si k, n ≥ N . En R toda sucesión de Cauchy es convergente. Es decir,
existe x ∈ R tal que lı́mn→∞ xn = x.
Veremos que estas dos propiedades también son válidas en el conjunto de las
funciones continuas de I en Rm . Comencemos por definir lo que es una sucesión
de Cauchy en este espacio. En primer lugar, la norma del supremo de una función
f ∈ C 0 (I)m es
kf k∞ = sup kf (t)k.
t∈I
Con esto, decimos que una sucesión de {fn } de funciones en C 0 (I)m es de Cauchy
si para todo ε > 0 existe N ∈ N tal que
kfn − fk k∞ < ε si k, n ≥ N.
Diremos también que una sucesión {fn } en C 0 (I)m converge a f ∈ C 0 (I)m si para
todo ε > 0 existe N ∈ N tal que kfn − f k∞ < ε si n ≥ N .
kfn − f k∞ < ε
Corolario 5.1 (Teorema del Punto Fijo de Banach). Cada operador contrac-
tante de C 0 (I)m en C 0 (I)m tiene un único punto fijo.
f0 = f y fn = T (fn−1 ) = · · · = T n (f ) para n ≥ 1.
Kk − Kn
≤ kT (f ) − f k∞
1−K
kT (f ) − f k∞ k
≤ K .
1−K
Tomemos ε > 0. Dado que K k → 0 cuando k → ∞ (pues K < 1), existe N ∈ N tal
que K k < kTε(1−K)
(f )−f k∞ para todo k ≥ N . Tenemos entonces que kfn − fk k∞ ≤ ε si
n, k ≥ N . Al ser una sucesión de Cauchy, el Teorema 5.3 nos dice que es convergente.
Como mencionamos arriba, el lı́mite es el único punto fijo de T .
3.2. Demostración del Teorema 5.2. Recordemos que sólo nos resta pro-
bar que el operador T definido en (38) tiene un único punto fijo. Para ello usaremos
el Teorema del Punto Fijo de Banach, pero primero demostraremos un resultado
auxiliar:
Lema 5.3. Supongamos que para cierto N ∈ N el operador
TN =T | ◦ ·{z
· · ◦ T}
N n veces
Lema 5.4. Bajo las hipótesis del Teorema 5.2, existe N ∈ N tal que T N es
contractante.
Demostración. Sean X, Y ∈ C 0 (I)m y sea L la constante de Lipschitz de la
función F con respecto a su segunda variable. Tenemos que
Z t
kT (X)(t) − T (Y )(t)k = kF (s, X(s)) − F (s, Y (s))k ds
t0
Z t
≤ L kX(s) − Y (s)k ds
t0
≤ L(t − t0 )kX − Y k∞ .
Derecho de autor. Prohibida su reproducción
122 5. SISTEMAS LINEALES DE PRIMER ORDEN
De manera similar,
Z t
kT 2 (X)(t) − T 2 (Y )(t)k ≤ L kT (X)(s) − T (Y )(s)k ds
t0
Z t Z s
≤ L2 kX(σ) − Y (σ)k dσ ds
t0 t0
L2 (t − t0 )2
≤ kX − Y k∞ .
2
Si repetimos este argumento vemos que para cada n ∈ N se tiene que
(t − t0 )n Ln
kT n (X)(t) − T n (Y )(t)k ≤ kX − Y k∞ .
n!
Denotamos por M la longitud del intervalo I. Tenemos entonces que
(M L)n
kT n (X) − T n (Y )k∞ ≤ kX − Y k∞ .
n!
an
Como n! tiende a cero cuando n → ∞ para cada a ∈ R, podemos escoger N ∈ N
N
tal que (ML)
N! < 1, de donde T N resulta ser contractante.
Esto termina la demostración del Teorema 5.2. Observemos que, dada la fun-
e
ción constante X(t) e converge a la
≡ X0 , la sucesión de funciones Xn = T n (X)
única solución del problema de Cauchy. Estas funciones se llaman aproximaciones
sucesivas de Picard.
Algunas observaciones:
Ejemplo 5.3 (Cadena con resortes). Se tiene una cadena con extremos fijos
formada por n cuentas conectadas por resortes idénticos de constante elástica k y
largo natural l0 . Cada cuenta de masa mi , i = 1 . . . n, tiene libertad para oscilar
horizontalmente, pero presenta un rozamiento lineal con el medio de constante ci ,
i = 1 . . . n. Para cada cuenta tenemos la ecuación de movimiento
mi x′′i = −ci x′i − k[(xi − xi+1 ) + (xi − xi−1 )], i ∈ {1, . . . , n},
donde x0 = xn+1 = 0. Matricialmente tenemos
′′ ′
m1 x1 c1 x1
.. .. . .. .
. . = − ..
mn xn cn xn
2 −1
x1
..
−1 2 . x2
.. .. .. ..
−k . . . . .
. xn−1
.. 2 −1 xn
−1 2
El sistema tiene la forma M X ′′ + CX ′ + KX = 0, con X ∈ Rn . Si hacemos el
cambio de variables
X ′ X′
Z= , Z = ,
X′ X ′′
Derecho de autor. Prohibida su reproducción
126 5. SISTEMAS LINEALES DE PRIMER ORDEN
Teorema 5.6. El conjunto {φk }nk=1 es una base de H, llamada base fundamen-
tal canónica. En consecuencia, dim(H) = n y la solución del sistema homogéneo
′
Xh (t) = A(t)Xh (t) t∈I
Xh (t0 ) = X0
está dada por
Xh = Φ(t)X0 = x10 φ1 (t) + · · · xn0 φn (t).
Demostración. Debemos probar que todo elemento de H se escribe como
combinación lineal de φ1 , . . . , φn y que estas funciones son linealmente indepen-
dientes.
Partamos por la primera afirmación. Sea Xh ∈ H la solución del sistema ho-
mogéneo
Xh′ (t) = A(t)Xh (t), para t ∈ I
Xh (t0 ) = X0 .
Escibimos X0 = x10 e1 + · · · + xn0 en con x10 , . . . , xn0 ∈ R. Demostraremos que Xh =
x10 φ1 + · · · xn0 φn . Para esto definimos
Z(t) = x10 φ1 (t) + · · · xn0 φn (t) = Φ(t)X0 .
Derivando obtenemos
Z ′ (t) = x10 φ′1 (t) + · · · xn0 φ′n (t),
de donde
Z ′ (t) = A(t) x10 φ1 (t) + · · · xn0 φn (t)
= A(t)Z(t)
pues φ′k = Aφk para cada k. Además, por la definición de los φk tenemos que
Z(t0 ) = x10 e1 + · · · + xn0 en = X0 .
Tenemos entonces que Z es solución del mismo sistema de ecuaciones diferenciales
que Xh , con las mismas condiciones iniciales. En virtud del Teorema de Existencia
y Unicidad,
Xh (t) = Z(t) = Φ(t)X0 para todo t ∈ I.
n
Pn Veamos ahora que {φk }k=1 son linealmente independientes. Suponemos que
k=1 αk φk (t) = 0 para todo t ∈ I. Debemos probar que, necesariamente, αk = 0
para todo k = 1, . . . , n. Observemos primero que
.. .. α1
X n
. . .
αk φk (t) =
φ1 (t) · · · φn (t) .. = Φ(t)α.
k=1 .. .. αn
. .
Derecho de autor. Prohibida su reproducción
4. ESTRUCTURA GEOMÉTRICA DE LAS SOLUCIONES 127
En virtud del Lema 5.5, la matriz Φ(t) es invertible para cada t y por lo tanto
α = 0.
Para generalizar los resultados anteriores, definamos una matriz fundamental
(no necesariamente la canónica) como
.. ..
. .
M (t) =
ψ1 (t) · · · ψn (t) ,
.. ..
. .
con t ∈ I, donde ψk es solución del sistema
′
ψk (t) = A(t)ψk (t) para t ∈ I
ψk (t0 ) = vk
y {vk }nk=1 es una base de Rn .
Corolario 5.2. Sea M (t) una matriz fundamental del sistema. Entonces det(M (t)) 6=
0 para todo t ∈ I y
X(t) = M (t)M −1 (t0 )X0 .
Demostración. Observemos que la matriz M (t0 ) es invertible y satisface
M (t0 )ek = vk para k = 1, . . . , n. Luego M (t0 )−1 vk = ek para k = 1, . . . , n. La
matriz
Ψ(t) = M (t)M (t0 )−1
satisface (39) pues
Ψ(t0 ) = M (t0 )M (t0 )−1 = In
Ψ′ (t) = M ′ (t)M (t0 )−1 = A(t)M (t)M (t0 )−1 = A(t)Ψ(t).
De acuerdo con el Lema 5.5, Φ(t) = Ψ(t) = M (t)M (t0 )−1 , de donde se sigue el
resultado.
Para el sistema no homogéneo aplicaremos el método de variación de paráme-
tros estudiado en el capı́tulo 3 para ecuaciones de orden superior. Recordemos la
regla para derivar un producto contenida en el Lema 5.1 y la propiedad de Φ dada
por la ecuación (39) en el Lema 5.5. Buscamos una solución particular de la forma
Xp (t) = Φ(t)F (t),
donde Φ es la matriz fundamental canónica y F es una incógnita. Xp debe satisfacer
B(t) = Xp′ (t) − A(t)Xp (t)
= (Φ(t)F (t))′ − A(t)Φ(t)F (t)
= Φ′ (t)F (t) + Φ(t)F ′ (t) − A(t)Φ(t)F (t)
= Φ(t)F ′ (t),
para lo cual basta que
F ′ (t) = Φ(t)−1 B(t)
Z t
F (t) = Φ−1 (s)B(s)ds.
t0
Como la serie del lado derecho es convergente a cero, la sucesión {sn (i, j)}∞
n=0 es
de Cauchy. Concluimos que la esponencial de una matriz está bien definida.
1. e0·t = eA·0 = I.
d At
2. dt e = AeAt .
A(t+s)
3. e = eAt eAs . En particular, eAt es invertible y su inversa es (eAt )−1 =
−At
e .
4. AB = BA si, y sólo si, BeAt = eAt B.
5. AB = BA si, y sólo si, eAt eBt = e(A+B)t .
Demostración. La propiedad 1. es inmediata de la definición.
2. En cada intervalo [−b, b], b ∈ R, las componentes de la matriz At están uni-
formemente acotadas por | máx(aij )||b|, luego tenemos la convergencia uniforme
Pp (Ak )ij tk
de hp (t) = k=0 k! . Como hp converge uniformemente a h y la sucesión
′
Pp (Ak )ij tk−1 Pp (Ak )ij tk−1
hp (t) = k=1 k k! converge uniformemente en (−b, b) a k=1 k k! ,
tenemos que h(t) es derivable y
∞
X
′ (Ak )ij tk−1
h (t) = k
k!
k=1
de donde
d At
(e )ij = Aij (eAt )ij para todo t ∈ (−b, b).
dt
Dado que b era arbitrario, el resultado es válido en todo R.
4. Por inducción se prueba que AB = BA si, y sólo si, Ak B = BAk para cada
k ≥ 1. Luego
∞
X ∞
X ∞
X
Ak tk BAk tk Ak tk
BeAt = B = = · B = eAt B.
k! k! k!
k=0 k=0 k=0
5. Se prueba como 3.
0 ··· λn .. .. ..
. . .
El cálculo de la exponencial de una matriz diagonalizable es sumamente sencillo.
∞ ∞
!
X (P DP −1 )k tk X D k tk
At P DP −1 t
e =e = =P P −1 = P eDt P −1 ,
k! k!
k=0 k=0
donde
eλ1 t ··· 0
eDt = ... ..
.
..
. .
0 ··· eλn t
De esta forma la solución del sistema homogéneo se puede escribir como:
Xh (t) = eA(t−t0 ) X0 = P eDt e−Dt0 P −1 X0 = P eDt C,
| {z }
C
C1
.
donde C = .. es un vector constante que depende de las condiciones iniciales.
Cn
Desarrollando el producto podemos escribir
Xh (t) = C1 eλ1 t v1 + C2 eλ2 t v2 + · · · + Cn eλn t vn .
Derecho de autor. Prohibida su reproducción
5. RESOLUCIÓN DE SISTEMAS LINEALES 131
Observación. Esta última fórmula puede obtenerse fácilmente viendo que ca-
da función de la forma eλt v es solución del sistema lineal siempre que v sea un vector
propio asociado al valor propio λ. Como toda matriz diagonalizable de tamaño n×n
tiene n vectores linealmente independientes, es evidente que las funciones de esta
forma generan todo el espacio de soluciones homogéneas. Sin embargo, preferimos
construir las soluciones a partir de la exponencial para que el lector se familiarice
con este método, imprescindible en el caso no diagonalizable.
Ejemplo 5.4 (Confort de un auto). Un modelo para estudiar el confort de un
auto está dado por
x′′ = −(k1 + k2 )x + (k1 L1 − k2 L2 )θ
G G
k1 k2 x(t) θ (t)
L2 L1
donde a = −(1 + µ)k1 y b = −µ(1 + µ)k1 L22 . Dado que tenemos la forma X ′ = AX,
podemos calcular la exponencial de la matriz. Para esto analicemos cómo son las
potencias de A:
0 0 1 0 a 0 0 0
0 0 0 1 0 b 0 0
A= 2
a 0 0 0 ⇒ A = 0 0 a 0 ;
0 b 0 0 0 0 0 b
luego, elevando a k,
ak 0 0 0
0 bk 0 0
A2k =
0
,
0 ak 0
0 0 0 bk
y si multiplicamos por A queda
0 0 ak 0
0 0 0 bk
A2k+1 =
ak+1
.
0 0 0
0 bk+1 0 0
Teniendo todas las potencias calculadas, procedemos a determinar la matriz expo-
nencial
k
a 0 0 0
X∞
t2k 0 b
k
0 1
eAt =
(2k)! 0 0 ak 0
k=0
0 0 0 bk
0 0 ak 0
∞
X t2k+1 0
0 0 bk .
+
(2k + 1)! a k+1
0 0 0
k=0
0 bk+1 0 0
Si tomamos a = −α2 y b = −β 2 obtenemos
∞
X ∞
X
ak t2k (−1)k (αt)2k
= = cos(αt)
(2k)! (2k)!
k=0 k=0
X∞ X∞
bk t2k (−1)k (βt)2k
= = cos(βt),
(2k)! (2k)!
k=0 k=0
y
∞
X ∞
ak t2k+1 1 X (−1)k (αt)2k+1 1
= = sen(αt)
(2k + 1)! α (2k + 1)! α
k=0 k=0
∞
X ∞
bk t2k+1 1 X (−1)k (βt)2k+1 1
= = sen(βt).
(2k + 1)! β (2k + 1)! β
k=0 k=0
Por lo tanto:
Derecho de autor. Prohibida su reproducción
5. RESOLUCIÓN DE SISTEMAS LINEALES 133
1
cos(αt) 0 αsen(αt) 0
0 cos(βt) 0 1
eAt = β sen(βt) .
−α sen(αt) 0 cos(αt) 0
0 −β sen(βt) 0 cos(βt)
Si ponemos alguna condición inicial, como por ejemplo una frenada, represen-
tada por
x(0) = 0, θ(0) = 0, x′ (0) = −1, θ′ (0) = −1,
obtenemos
x − sen(αt)
α
θ − sen(βt)
′ =
.
x β
− cos(αt)
θ′ − cos(βt)
Es decir, la solución del movimiento es:
− sen(αt) − sen(βt)
x(t) = , θ(t) = ,
α β
p p
con α = (1 + µ)k y β = L µ(1 + µ)k, que fı́sicamente representan las frecuen-
cias de cada modo.
Otro método para solucionar este problema es analizar los valores y vectores
propios del sistema.
−λ 0 1 0
0 −λ 0 1
det(A − λI) = det
a
= (λ2 − a)(λ2 − b).
0 −λ 0
0 b 0 −λ
En términos de las variables antes definidas quedan los valores propios imaginarios
puros:
λ1 = iα, λ2 = −iα, λ3 = iβ, λ4 = −iβ,
que tienen asociados respectivamente los vectores propios:
−i/α i/α 0 0
0 0 −i/β i/β
v1 =
1 , v2 = 1 , v3 = 0 ,
v4
0 .
0 0 1 1
Luego la solución general será:
x −i/α i/α
θ
′ = c1 eiαt 0 + c2 e−iαt 0
x 1 1
θ′ 0 0
0 0
−i/β i/β
+c3 eiβt
0 + c4 e
−iβt
0
.
1 1
Derecho de autor. Prohibida su reproducción
134 5. SISTEMAS LINEALES DE PRIMER ORDEN
x(t) (ic2 e−iαt − ic1 eiαt ) 1 (ic4 e−iβt − ic3 eiβt ) 0
= + .
θ(t) α 0 β 1
−c1 + c2 =0
−c3 + c4 =0
−i −i
α c1 + c2 α = −1
−i
β c3 + c4 −i
β = −1 ,
x(t) i 1/α i 0
= (−e−iαt + eiαt ) + (−e−iβt + ieiβt ) ;
θ(t) 2 0 2 1/β
8 0 0 0 0 0
0 8 0 0 0 0
0 0 3 0 0 0
.
0 0 0 3 1 0
0 0 0 0 3 0
0 0 0 0 0 4
Para λ2 calculamos:
l1 = dim ker(A − I) = 2
l2 = dim ker(A − I)2 = 4
ls = dim ker(A − I)s = 4, ∀s ≥ 2.
Como l0 = 0, tenemos que k1 = 0 y k2 = 2. En virtud de la Proposición 5.6 no
hay cadenas de largo 1, pero hay 2 cadenas de largo 2. Cada cadena de largo 2
determina un bloque de 2 × 2, por lo que tenemos completamente determinada la
matriz J:
1 1 0 0 0
0 1 0 0 0
J = 0 0 1 1 0 .
0 0 0 1 0
0 0 0 0 −1
Buscamos ahora los vectores propios generalizados. Para λ2 tenemos
(A − I)v1 = 0 ⇒ v1 = (α, α, 0, β, β)
(A − I)v2 = v1 ⇒ v2 = (β + γ, γ, α, α + δ, δ).
Si α = 0, β = 1 ⇒ v1 = (0, 0, 0, 1, 1) ⇒ v2 = (1 + γ, γ, 0, δ, δ). Si γ = δ = 0 ⇒ v2 =
(1, 0, 0, 0, 0).
Si α = 1, β = 0 ⇒ v1 = (1, 1, 0, 0, 0) ⇒ v2 = (γ, γ, 1, 1 + δ, δ). Si γ = δ = 0 ⇒ v2 =
(0, 0, 1, 1, 0).
La matriz P de cambio de base tiene por columnas los vectores propios generaliza-
dos:
0 1 1 0 0
0 0 1 0 1
P = 0 0 0 1 1 .
1 0 0 1 0
1 0 0 0 0
N0 = I
N1 = N
0 1 0 0 0 1 0 0
.. .. .. .. .. . . .. ..
. . . . . . . .
N 2
=
.. .. .. .. . . ..
. . . 1 . . . 1
0 ··· ··· 0 0 ··· ··· 0
0 0 1 0
.. .. ..
. . . 1
=
.
.. .. ..
. . . 0
0 ··· ··· 0
∞
X m−1
X
N k tk N k tk
eN t = =
k! k!
k=0 k=0
2
t tm−1
= I + tN + N 2 + · · · + N m−1
2! (m − 1)!
t2 tm−1
1 t 2! · · · (m−1)!
. ..
. .. .. ..
. . . . .
.
= . ... ... ... t2
.
.
2!
. . . .
.. .. .. .. t
0 ··· ··· ··· 1
Finalmente,
t2 tm−1
1 t 2! · · · (m−1)!
. ..
. .. .. ..
. . . . .
eJt = eλt
.
.. . .
.
..
.
..
. t 2 .
2!
. . .. ..
.. .. . . t
0 ··· ··· ··· 1
De esta forma hemos demostrado el siguiente resultado:
Proposición 5.8. Sea A ∈ Mn×n (C), y sea A = P JP −1 su representación en
forma canónica de Jordan. Entonces
J1 t
e ... 0
.. P −1 .
eAt = P ... ..
. .
Jn t
0 ... e
Ejemplo 5.7. Consideremos el sistema:
X ′ = AX,
con condición inicial X(0) = X0 , donde la matriz A está dada en el Ejemplo 5.6.
Ya habı́amos calculado su forma canónica de Jordan. En virtud de la proposición
anterior se tiene que la solución es:
t
e tet 0 0 0
0 et 0 0 0
−1
X(t) = P 0 t
0 e te t
0
P X0 ,
0 0 0 et 0
0 0 0 0 e−t
donde P es la matriz de cambio de base descrita en el Ejemplo 5.6.
Probemos por inducción que efectivamente la solución del sistema está dada por:
n
X
n
Xn = A X0 + An−j Bj−1 , n≥1
j=1
Para n = 1 resulta
X1 = AX0 + B0 ,
Pn
que es solución del sistema. Suponemos ahora que Xn = An X0 + j=1 An−j Bj−1
satisface el sistema propuesto. Hay que probar que
n+1
X
Xn+1 = An+1 X0 + An+1−j Bj−1
j=1
= AXn + Bn .
Estudiemos ahora la solución del sistema discreto homogéneo, para el caso en
que A sea diagonalizable, es decir A = P DP −1 , con
λ1 · · · 0
.. , donde λ , . . . λ son los valores propios de A.
D = ... ..
. . 1 d
0 · · · λd
Derecho de autor. Prohibida su reproducción
142 5. SISTEMAS LINEALES DE PRIMER ORDEN
De esta forma vemos que se repite la solución que tenı́amos para el caso diagonali-
zable: independiente de la condiciones iniciales, nuestro modelo indica que cuando
n → ∞ , si |λ| > 1, el sistema diverge, es decir las poblaciónes de ballenas y
plancton se expanden indefinidamente, al igual que la temperatura, mientras que
si |λ| < 1, las poblaciones se extinguen y la temperatura del mar desciende a 0.
En el caso crı́tico |λ| = 1 se tiene que las poblaciones animales divergen, pero con
temperatura constante del mar T0 .
4 3 2 1
entonces
1 2 2 1 3 3
1 Kt 2K t 6K t
−Kt 0 1 Kt 1 2 2
eAt =e 2K t .
0 0 1 Kt
0 0 0 1
Se deja al lector escribir la solución general usando la fórmula de variación de
parámetros y analizar los casos en que hay algún volumen (y por lo tanto algún
valor propio) repetido.
Si λ ≥ 0 entonces θ ∈ [0, 1] y
−b(1 + λ) −b(1 + λ)
s1 = (1 + θ) < 0 y s2 = (1 − θ) < 0.
V V
Si introducimos los parámetros
b(1 + λ) bλ
α= , β= ,
V V
vemos que
bσ + αx01 + βx02 sx01
Lx1 = +
(s + α(1 + θ))(s + α(1 − θ)) (s + α(1 + θ))(s + α(1 − θ))
αbσ
+
s(s + α(1 + θ))(s + α(1 − θ))
α(x01 + x02 ) sx02
Lx2 = +
(s + α(1 + θ))(s + α(1 − θ)) (s + α(1 + θ))(s + α(1 − θ))
αbσ
+ .
s(s + α(1 + θ))(s + α(1 − θ))
Los tres sumandos de cada ecuación anterior tienen una antitransformada conocida:
1 eat − ebt
L−1 =
(s − a)(s − b) a−b
s aeat − bebt
L−1 =
(s − a)(s − b) a−b
1 (c − b)eat + (a − c)ebt + (b − a)ect
L−1 = .
(s − a)(s − b)(s − c) (a − b)(b − c)(c − a)
Por lo tanto las soluciones son:
e−α(1+θ)t − e−α(1−θ)t
x1 = −(bσ + αx01 + βx02 )
2αθ
−α(1+θ)t
(1 + θ)e − (1 − θ)e−α(1−θ)t
+x01
2θ
(1 − θ)e−α(1+θ)t − (1 + θ)e−α(1−θ)t + 2θ
+αbσ
2α2 θ(1 − θ2 )
e−α(1+θ)t − e−α(1−θ)t
x2 = −α(x01 + x02 )
2αθ
−α(1+θ)t
(1 + θ)e − (1 − θ)e−α(1−θ)t
+x02
2θ
(1 − θ)e−α(1+θ)t − (1 + θ)e−α(1−θ)t + 2θ
αbσ .
2α2 θ(1 − θ2 )
Luego
x1 = C1 es1 t + C2 es1 t + C3
x2 e1 es1 t + C
= C e2 es1 t + C
e3 ,
donde C1 , Ce1 , C2 , C
e2 , C3 , C
e3 son constantes dadas de la agrupación de términos
semejantes en las expresiones anteriores, y que dependen sólo de las condiciones
iniciales y de la geometrı́a del sistema.
Derecho de autor. Prohibida su reproducción
146 5. SISTEMAS LINEALES DE PRIMER ORDEN
λ=0
λ>0
0 10 20 30 40 50 60
c1 N
ecuador
c2
S
Introduzcamos la masa media entre los dos hemisferios como c(t) = 12 (c1 (t) +
c2 (t)) y la emisión media como f = 12 (f1 + f2 ). Si derivamos la expresión de la masa
media con respecto al tiempo obtenemos
1 ′
(c (t) + c′2 (t)).
c′ (t) =
2 1
Luego, si sumamos las EDO que tenemos para c1 y c2 , nos queda
c′1 + c′2 = 2c′ = (f1 + f2 ) − α(c1 − c2 + c2 − c1 ) − β(c1 + c2 )
c′ = f − βc
Derecho de autor. Prohibida su reproducción
148 5. SISTEMAS LINEALES DE PRIMER ORDEN
que es una EDO lineal de primer orden para c. La condición inicial está dada por
1 1
c(0) = (c1 (0) + c2 (0)) = (84 + 60)[kton] = 72[kton].
2 2
La solución del problema de Cauchy asociado es
f
c(t) = 72e−βt + (1 − e−βt ).
β
Se estima que el lı́mite de c(t) cuando t → +∞ es de 100 [kton]. De allı́ podemos
encontrar una estimación para el tiempo de vida del contaminante
f f
lı́m c(t) = 100[kton] = lı́m 72e−βt + (1 − e−βt ) = .
t→∞ t→∞ β β
100[kton]
β −1 = = 5 [años].
20[kton/año]
El tiempo estimado de vida del contaminante es de 5 años.
k1 φ1 (s) + k2 φ2 (s) + k3 φ3 (s)
= e ,
k1 φ1 (s) + e
k2 φ2 (s) + e
k3 φ3 (s)
donde
1 1
L−1 [φ1 (s)] = L−1 −
2α(s + β) 2α(s + 2α + β)
1 −1 1 1 −1 1
= L − L
2α s+β 2α s + 2α + β
1 −βt 1 −(2α+β)t
= e − e ,
2α 2α
1 −1 s 1 −1 s
L−1 [φ2 (s)] = L − L
2α s+β 2α s + 2α + β
1
= (δ(t) − βe−βt − δ(t) + (2α + β)e−(2α+β)t )
2α
1
= (2α + β)e−(2α+β)t − βe−βt ,
2α
1 1 1 −1 1
L−1 [φ3 (s)] = L−1 − L
β(2α + β) s 2αβ s+β
1 1
+ L−1
2α(2α + β) s + 2α + β
1 1 −βt 1
= − e + e−(2α+β)t .
β(2α + β) 2αβ 2α(2α + β)
Finalmente, agrupando constantes, la cantidad de contaminante en cada hemisferio
resulta:
Capı́tulo 6
151
Derecho de autor. Prohibida su reproducción
152 6. SISTEMAS NO LINEALES
x′ = y
′ c g
y = − y − sen(x) + f (t).
m L
La no linealidad se aprecia claramente en el término sen(θ) pues
x′ = y
c g
y′ = − y − sen(x).
m L
Decimos que (x̄, ȳ) es un punto crı́tico o de equilibrio del (SN LA) si
F (x̄, ȳ) = 0
G(x̄, ȳ) = 0.
Se llama nulclina en x a la curva definida por F (x, y) = 0 y nulclina en y a la
curva definida por G(x, y) = 0.
Los puntos crı́ticos son las intersecciones de las nulclinas. El conjunto de puntos
crı́ticos del sistema se denota por
C = {(x̄, ȳ) ∈ R2 | F (x̄, ȳ) = G(x̄, ȳ) = 0}.
Si hacemos una expansión de Taylor en el (SN LA) en torno a un punto crı́tico
(x̄, ȳ), resulta:
∂F ∂F
F (x, y) = (x̄, ȳ)(x − x̄) + (x̄, ȳ)(y − ȳ) + hF (r)
∂x ∂y
∂G ∂G
G(x, y) = (x̄, ȳ)(x − x̄) + (x̄, ȳ)(y − ȳ) + hG (r),
∂x ∂y
donde r = (x, y) − (x̄, ȳ) (recordemos que F (x̄, ȳ) = G(x̄, ȳ) = 0). Además
hF (r) hG (r)
lı́m =0 y lı́m = 0.
r→0 krk r→0 krk
2. Si |J(x̄, ȳ)| = 0 el sistema tiene infinitas soluciones, que forman una recta si
J(x̄, ȳ) 6= 0 y todo el plano si J(x̄, ȳ) = 0.
3. Por contradicción, supongamos que para cada δ > 0 existe otro punto crı́tico
(w̄, z̄) de (SN LA) tal que |(x̄, ȳ) − (w̄, z̄)| < δ. Por definición,
a(w̄ − x̄) + b(z̄ − ȳ) + hF (r̄) = 0
c(w̄ − x̄) + d(z̄ − ȳ) + hG (r̄) = 0,
donde r̄ = (w̄, z̄) − (x̄, ȳ). De manera equivalente,
a b w̄ − x̄ hF (r̄)
=− .
c d z̄ − ȳ hG (r̄)
Si |J(x̄, ȳ)| 6= 0 tenemos que
−1
w̄ − x̄ a b hF (r̄)
=− .
z̄ − ȳ c d hG (r̄)
Tomando norma y acotando obtenemos
p
||r̄|| ≤ ||J(x̄, ȳ)−1 || hF (r̄)2 + hF (r̄)2
s 2 2
−1 hF (r̄) hF (r̄)
1 ≤ ||J(x̄, ȳ) || + .
r̄ r̄
Derecho de autor. Prohibida su reproducción
1. SISTEMAS NO LINEALES Y SISTEMAS LINEALIZADOS 155
Pero el lado derecho tiende a cero cuando r̄ → 0, con lo que llegamos a una contra-
dicción.
Tenemos entonces que C = {(0, 0), (0, 14), (20, 0), (12, 6)}, conjunto que muestra
las posibilidades de equilibrio entre ambas poblaciones. El único caso en que pueden
coexistir ambas especies es (x̄, ȳ) = (12, 6). En los demás casos, al menos una de
las especies se extingue. Aun ası́, todas son soluciones de equilibrio. El Jacobiano
y
(0,14)
(12,6)
(20,0)
(0,0) x
Figura 1. Soluciones de equilibrio para el modelo de conejos y ovejas
Derecho de autor. Prohibida su reproducción
156 6. SISTEMAS NO LINEALES
Ejemplo 6.4. Habı́amos visto que las ecuaciones de un péndulo sin forzamiento
están dadas por ′
x = y
y′ = − mc
y − Lg sen(x).
Sus puntos crı́ticos satisfacen
0 = ȳ
0 = −m c
ȳ − Lg sen(x̄),
de donde
C = {(kπ, 0) | k ∈ Z}.
x
(− 2π,0)(−π,0) (0,0) (π,0) (2π,0)
Figura 2. Soluciones de equilibrio del péndulo no lineal
Derecho de autor. Prohibida su reproducción
1. SISTEMAS NO LINEALES Y SISTEMAS LINEALIZADOS 157
El Jacobiano es
0 1
J(x̄, ȳ) = ,
− Lg (−1)k −mc
y 2
–4 –2 x
2 4 6
–2
–4
Figura 3. Intersección de Cónicas
(xo,yo)
Denotamos este punto por (x0 , y0 ). La trayectoria T3 definida por (x3 (t), y3 (t)) =
(x2 (t + t2 − t1 ), y2 (t + t2 − t1 )) tiene el mismo recorrido que T2 porque es una tras-
lación en tiempo y el sistema es autónomo. Pero también tiene el mismo recorrido
que T1 porque es solución del mismo problema de Cauchy.
Esto nos dice que las curvas en el diagrama de fase no se intersectan, de modo
que una situación como la presentada en la Figura 4 es imposible.
Ejemplo 6.7. Lo que sı́ se puede dar para dos trayectorias distintas, es que
sus recorridos se confundan. Por ejemplo, tomemos la trayectorias:
(42) (x1 (t), y1 (t)) = (sen(t), cos(t)) con (x(0), y(0)) = (0, 1)
(43) (x2 (t), y2 (t)) = (cos(t), sen(t)) con (x(0), y(0)) = (1, 0).
Elevando al cuadrado cada componente y sumando vemos que ambas trayectorias
satisfacen x2 + y 2 = 1. También es fácil ver que recorren toda la circunferencia. Sin
embargo lo hacen en sentidos opuestos, como se ilustra en la figura 5. Si queremos
hacer que las dos trayectorias recorran la curva en el mismo sentido basta modificar
(43) a
(e
x2 (t), ye2 (t)) = (cos(t), − sen(t)) con (x(0), y(0)) = (1, 0).
y y
x x
Por otra parte, un punto crı́tico (x̄, ȳ) es asintóticamente estable si es estable
y además existe δ > 0 tal que
lı́m (x(t), y(t)) = (x̄, ȳ), siempre que k(x0 , y0 ) − (x̄, ȳ)k < δ.
t→∞
Notemos que no todo punto estable es asintóticamente estable.
xo
x
yo
y
que representan rectas de pendiente xy00 , cuyo gráfico en el plano XY queda ilustra-
do en la Figura 7, donde las flechas indican el sentido positivo del tiempo. Cuando
t → ∞ tenemos que (x(t), y(t)) → 0.
t=0
yo
y
t=1
t=2
xo
x
yo
y
xo
x
Se dice que una trayectoria t 7→ x(t), y(t) converge o tiende al punto crı́tico
(x̄, ȳ) cuando t → ∞ si
lı́m x(t) = x̄ y lı́m y(t) = ȳ.
t→∞ t→∞
Análogamente se define la convergencia cuando t → −∞.
entonces se dice que la trayectoria entra al punto crı́tico (x̄, ȳ) tangente a una
semirrecta de pendiente l cuando t → ∞. De la misma forma, si además de converger
al punto (x̄, ȳ) cuando t → −∞ la trayectoria es tal que existe
y(t) − ȳ
l = lı́m
t→−∞ x(t) − x̄
se dice que la trayectoria sale del punto crı́tico (x̄, ȳ) tangente a una semirrecta de
pendiente l cuando t → −∞.
Un punto crı́tco aislado (x̄, ȳ) se llama nodo si todas las trayectorias vecinas,
o bien entran al punto (x̄, ȳ), o bien salen de él.
Derecho de autor. Prohibida su reproducción
164 6. SISTEMAS NO LINEALES
yo
y
x
xo
x
Un punto crı́tco aislado (x̄, ȳ) se llama punto silla si existen dos trayectorias
que entran a él tangentes a semirrectas opuestas y dos que salen de la misma for-
ma. Las restantes trayectorias no convergen a (x̄, ȳ) y tienen como ası́ntotas a las
semirrectas antes mencionadas.
Diremos que un punto crı́tico es punto espiral si todas las trayectorias en una
vecindad del punto convergen pero no entran cuando t → ∞, o convergen pero no
salen cuando t → −∞. Un punto crı́tico es un centro si todas las trayectorias en
Derecho de autor. Prohibida su reproducción
3. CLASIFICACIÓN DE LOS PUNTOS CRÍTICOS 165
una vecindad del punto son cerradas y existen trayectorias arbitrariamente cerca
del punto.
y y
x x
y y
x x
Las Figuras 13, 14 y 15 nos dicen también las caracterı́sticas de estabilidad del
punto crı́tico (0, 0): si α < 0 es un punto espiral asintóticamente estable; si α > 0
es un punto espiral inestable; si α = 0 es un centro estable. Si β = 0 es un nodo
(estable o inestable según el signo de α).
Derecho de autor. Prohibida su reproducción
4. PUNTOS CRÍTICOS DE SISTEMAS LINEALES 167
y y
x x
A. Casos Principales:
1. λ1 , λ2 reales distintos pero de igual signo.
2. λ1 , λ2 reales distintos y de distinto signo.
3. λ1 , λ2 complejos (conjugados), con parte real no nula.
B. Casos Frontera:
1. λ1 , λ2 reales e iguales.
2. λ1 , λ2 imaginarios puros (conjugados).
A. Casos Principales:
En todos estos casos los valores propios son distintos, por ende la matriz es
diagonalizable, es decir, se puede escribir de la forma A = P DP −1 , con
λ1 0
D= P = v1 | v2 ,
0 λ2
donde v1 , v2 son los vectores propios asociados a λ1 y λ2 , respectivamente. De esta
forma la solución del sistema lineal es
x At x0
= e
y y0
λt
x e 1 0 −1 x0
= P P
y 0 eλ2 t y0
λt
x e 1 0 x0
P −1 = P −1
y 0 eλ2 t y0
| {z } | {z }
λt
x
e e 1
0 x
e0
= .
ye 0 eλ2 t ye0
De esta forma, el problema se desacopla en estas nuevas coordenadas x e, ye intro-
ducidas por las direcciones de los vectores propios (llamadas también direcciones
propias). Resulta simplemente
e = eλ1 t x
x e0
ye = eλ2 t ye0 .
Observemos que si, por ejemplo, λ1 < 0, entonces
e = lı́m eλ1 t x
lı́m x e0 = 0.
t→∞ t→∞
Si λ1 > 0, entonces
e = lı́m eλ1 t x
lı́m x e0 = ±∞,
t→∞ t→∞
donde el signo de ±∞ está dado por el signo de la condición inicial. El mismo
argumento es válido para λ2 .
Sólo hay dos casos: λ2 /λ1 > 1 o λ2 /λ1 < 1. El bosquejo del diagrama de fases
para el caso λ2 /λ1 > 1 se tiene en la Figura 16. En el gráfico de la izquierda, para
λ2 < λ1 < 0 se tiene estabilidad asintótica, mientras que en el de la derecha, para
0 < λ1 < λ2 , el origen es inestable.
Para el caso λ2 /λ1 < 1, se tiene la Figura 17. En el gráfico de la izquierda, para
λ1 < λ2 < 0 se tiene estabilidad asintótica, mientras que en el de la derecha, para
0 < λ2 < λ1 , el origen es inestable.
Derecho de autor. Prohibida su reproducción
4. PUNTOS CRÍTICOS DE SISTEMAS LINEALES 169
y y
x x
y y
x x
Como conclusión tenemos que en ambos casos el punto (0, 0) es un nodo. Si los
valores propios son negativos, es asintóticamente estable; si son positivos, entonces
es inestable. Además, vemos que en general para el caso de valores propios reales
distintos de igual signo, los recorridos de las trayectorias son siempre tangentes a
la recta dada por la dirección del vector propio asociado al valor propio de menor
módulo.
En efecto:
λ2
|λ1 | < |λ2 | ⇒ λ1 > 1 ⇒ recorrido tangente a v1 .
λ2
|λ2 | < |λ1 | ⇒ λ1 < 1 ⇒ recorrido tangente a v2 .
Los diagramas de fase serán los de la Figura 18. El punto crı́tico es un punto
silla, que por supuesto es inestable. Las trayectorias entran en la dirección propia
asociada al valor propio negativo y salen en la dirección propia asociada al valor
propio positivo.
Ejemplo 6.13 (Continuación de conejos y ovejas). Volvamos al ejemplo de los
conejos compitendo con las ovejas:
′
x = 60x − 3x2 − 4xy ≡ F (x, y)
y ′ = 42y − 3y 2 − 2xy ≡ G(x, y).
Sabemos que
C = {(0, 0), (0, 14), (20, 0), (12, 6)}.
Derecho de autor. Prohibida su reproducción
170 6. SISTEMAS NO LINEALES
y y
x x
Vimos que este sistema no era degenerado entorno a ninguno de estos cuatro puntos
crı́ticos, que ahora clasificaremos con respecto al sistema linealizado. Ya habı́amos
evaluado el jacobiano en cada punto crı́tico:
60 0
1. J(0, 0) = . Sus valores propios son λ1 = 60 y λ2 = 42, por lo
0 42
que se trata de un nodo inestable. Los vectores propios asociados son
1 0
v1 = v2 = ,
0 1
respectivamente. Los recorridos de las trayectorias son tangentes a la di-
rección v2 (salvo
los 2 queson tangentes a v1 ).
−60 −80
2. J(20, 0) = . Sus valores propios son λ1 = −60 y λ2 = 2,
0 2
por lo que éste es un punto silla. Los vectores propios asociados son
1 80
v1 = v2 = ,
0 −62
respectivamente.
El eje delas trayectorias convergentes es la dirección v1 .
4 0
3. J(0, 14) = . Sus valores propios son λ1 = 4 y λ2 = −42,
−28 −42
por lo que éste también es un punto silla. Los vectores propios asociados
son
46 0
v1 = v2 = ,
−28 1
respectivamente.
El eje delas trayectorias convergentes es la dirección v2 .
−36 −48 √
4. J(12, 6) = . Sus valores propios son λ1 = −27 − 3 73 ≈
−12 −18
√
−52,6 y λ2 = −27 + 3 73 ≈ −1,4 por lo que el punto crı́tico (12, 6) es un
nodo asintóticamente estable. Los vectores propios asociados son
16√ 16√
v1 = v2 = ,
−3 + 73 −3 − 73
respectivamente. Los recorridos de trayectorias son tangentes a la dirección
v2 (salvo los 2 que son tangentes a v1 ).
Derecho de autor. Prohibida su reproducción
4. PUNTOS CRÍTICOS DE SISTEMAS LINEALES 171
b b
c g
k impar: Tenemos que λ2 + mλ −= 0, de donde
L
r
c c2 g
λ=− ± 2
+ .
2m 4m L
Puesto que todas las contantes son positivas tenemos dos valores propios reales de
distinto signo. Los puntos crı́ticos resultan ser puntos sillas inestables.
c g
k par: Aquı́ λ2 + mλ + L = 0 y luego
r
c c2 g
λ=− ± − .
2m 4m2 L
Analicemos por casos el signo de la cantidad subradical:
2
c g
1. Sobreamortiguado: El caso 4m 2 > L nos entrega dos valores propios reales
distintos, ambos negativos, por lo que los puntos crı́ticos (kπ, 0) resultan
nodos asintóticamente estables.
c2 g
2. Subamortiguado: El caso 4m 2 < L nos entrega dos valores propios complejos
conjugados, por lo que los puntos crı́ticos (kπ, 0) resultan puntos espirales,
c
que además son asintóticamente estables pues α = − 2m < 0.
Este caso corresponde al más común, que ocurre cuando el coeficiente de
c2 g
pg
roce es pequeño (notar que 4m 2 < L si, y sólo si, c < 2m L ).
El sentido de giro de las espirales se puede determinar calculando el flujo
el algunos puntos vecinos (ver flechas en la Figura 19).
Notemos que si el péndulo está en posición vertical hacia arriba, que
corresponde a los puntos k impar, está en un equilibrio inestable. Intuiti-
vamente, si en esta posición lo perturbamos ligeramente, el péndulo oscila
hasta la posición vertical hacia abajo, que son los puntos espirales asintóti-
camente estables con k par (ver situación cerca de (π, 0) en la Figura 19).
4
y
2
−π 00 π x
− 3π − 2π 2π 3π
–2 x
–4
2
c g
3. Crı́ticamente amortiguado: El caso 4m 2 = L nos entrega un sólo valor
B. Casos Frontera.
Tenemos el sistema
x′ = ax + by
y′ = cx + dy,
a b
con la matriz A = no es necesariamente diagonalizable (pues tiene el
c d
valor propio repetido), pero siempre se puede llevar a su forma canónica de Jordan
A = P JP −1 , donde J puede tener uno o dos bloques. En el caso de tener dos
bloques, la matriz es diagonalizable, con
λt
λ 0 Jt e 0
J= y de aquı́ e = .
0 λ 0 eλt
Luego, la ecuación en las nuevas coordenadas dadas por los vectores propios gene-
ralizados v1 , v2 es:
λt
x
e e 0 x
e0 e = eλt x
x e0
= λt ⇒ .
ye 0 e ye0 ye = eλt ye0
El diagrama de fase de este caso ya es bien conocido del comienzo del capı́tulo (ver
Figura 8 para caso λ < 0). El punto crı́tico será un nodo asintóticamente estable si
el valor propio es negativo y un nodo inestable si es positivo. Este tipo de nodos se
denomina nodos estrella.
más adelante sobre este mismo punto cuando analicemos la energı́a de sistemas
conservativos.
4 det(A) = tr(A)2
det(A)
1 2
4 3
tr(A)
5
Plano traza - determinante.
En el semiplano det(A) < 0 (región 5) las raı́ces son reales y tienen distinto
signo, por lo que tendremos sólo puntos sillas inestables.
Sobre la parábola tenemos una raı́z real con multiplicidad 2. Se tienen nodos
inestables si tr(A) < 0 y nodos asintóticamente estables si tr(A) > 0.
Hemos concluido ası́ una forma que nos permite tener una idea general del
punto crı́tico de un sistema sólo viendo la traza y determinante de la matriz, sin
tener que calcular valores y vectores propios explı́citamente. Sin embargo, si que-
remos un gráfico preciso del diagrama de fases es inevitable calcular los vectores
propios, pues indican las direcciones de tangencia de las trayectorias, para el caso
que corresponda.
entonces
|x′ (t1 )|2 |x′ (t1 )|2
+ F (x1 ) = + F (x2 ).
2 2
Esto quiere decir que la siguiente cantidad (que es exactamente la energı́a salvo una
constante multiplicativa y/o aditiva) es conservada a lo largo de las trayectorias:
y2
E(x, y) = + F (x).
2
Las curvas de iso-energı́a quedan ası́ parametrizadas por E de la forma
p
y = ± 2(E − F (x)).
Derecho de autor. Prohibida su reproducción
5. ENERGÍA, FUNCIONES DE LIAPUNOV Y ESTABILIDAD 177
−1
−2
−3
−4
−15 −10 −5 0 5 10 15
2
Notar que la energı́a cinética de este sistema es K = 12 mL2 θ′ = 21 mL2 y 2 , y
la energı́a potencial U = mgL(1 − cos(θ)) = mgL(1 − cos(x)), entonces la energı́a
total es
1
V (x, y) = K + U = mL2 y 2 + mgL(1 − cos(x))
2
y podemos verificar que V y E se diferencian solamente en una constante multipli-
cativa y una constante aditiva:
V (x, y) = mL2 E(x, y) + mgL
de manera que las curvas de nivel de E y V son cualitativamente iguales (ver Figura
20) y da lo mismo cuál de las dos funciones energı́a considerar.
Si queremos verificar directamente que V es conservada a lo largo de las tra-
yectorias, derivemos la función V con respecto al tiempo usando la regla de la
cadena:
dV ∂V dx ∂V dy
= + = mgL sen(x)x′ + mL2 yy ′ .
dt ∂x dt ∂y dt
Derecho de autor. Prohibida su reproducción
178 6. SISTEMAS NO LINEALES
2. Para la establidad asintótica debemos δ > 0 de manera que si (x0 , y0 ) ∈ Bδ (x̄, ȳ)
entonces lı́mt→∞ (x(t), y(t)) = (x̄, ȳ). Como la función
t 7→ V (x(t), y(t))
es positiva y decreciente para cualquier condición inicial, ella tiene un lı́mite L
cuando t → ∞. Si L = 0, la continuidad de V asegura que lı́mt→∞ (x(t), y(t)) =
Derecho de autor. Prohibida su reproducción
180 6. SISTEMAS NO LINEALES
(x̄, ȳ) que es el único punto donde V se anula. Probaremos ahora que L no puede
ser positivo. En efecto, si L > 0 entonces V (x(t), y(t)) ≥ L para todo t. Por otro
lado, siempre se tiene que V (x(t), y(t)) ≤ V (x0 , y0 ). Además, el punto 1 dice que
(x̄, ȳ) es estable, de modo que existe δ > 0 tal que si (x0 , y0 ) ∈ Bδ (x̄, ȳ) entonces
(x(t), y(t)) ∈ Br/2 (x̄, ȳ). En virtud de la continuidad de V el conjunto
K = { (x, y) ∈ B r/2 (x̄, ȳ) | L ≤ V (x, y) ≤ V (x0 , y0 ) }
es cerrado y acotado. Además contiene a la trayectoria y no contiene a (x̄, ȳ). Luego
d
máx V (x(t), y(t)) | (x(t), y(t)) ∈ K = −k < 0.
dt
Tenemos entonces que
V (x(t), y(t)) ≤ V (x0 , y0 ) − kt,
cantidad que tiende a −∞ cuando t → ∞. Esto es imposible puesto que V no toma
valores negativos.