UNIVERSIDAD DE PAMPLONA

TRABAJO FINAL DE METODOS NUMERICOS

PROFESOR: FRANCISCO CABRERA

ALUMNO: CARLOS ALBERTO HERNANDEZ
CODIGO: 1094264563
PROGRAMA: FISICA

PAMPLONA 09 DE AGOSTO DE 2019

 1) Para el modelo depredador-presa
a) Hacer un código que llame la matriz de cofactores A y resuelva el
sistema no lineal por runge-kutta
b) Dados los parámetros alfa, beta, gama y teta, obtener el sistema
linealizado y comparar con el inciso a).
Solución
Ecuaciones de Lotka-Volterra: Modelo depredador presa

Resulta intuitivo pensar que las poblaciones de un animal depredador y su presa

están relacionadas de algún modo en el que, si una aumenta, la otra lo hace
también. Utilizaremos como ejemplo en este artículo un ecosistema aislado y
formado por leones y cebras que viven en armonía, es decir, los leones
comiéndose a las cebras.
- Ecuaciones:

Las ecuaciones consisten en dos ecuaciones de primer orden, acopladas,

autónomas y no lineales:
dx
=αx−βxy
dt
dy
=−γy+ δyx
dt
Donde:
x=¿ numero de presas.
y=¿ numero de depredadores.
α : tasa de crecimiento de las presas.
β : éxito en la casa del depredador.
γ :tasa de crecimiento en los depredadores.
δ : éxito en la casa y cuanto alienta cazar una presa al depredador.
Asignándole valores a los anteriores parámetros se obtuvieron los siguientes
gráficos:
α =0.4 , β=0.37 , γ =0.3 , δ =0.05
α =0.4 , β=0.37 , γ =0.3 , δ =0.05
  Linealización del problema del depredador-presa:

Buscamos los puntos de equilibrio derivando con respecto al tiempo la población

de los depredadores y de la presa y las igualamos a cero:
dx
=αx−βxy =0
dt
dy
=−γy + δyx=0
dt
De aquí obtenemos las coordenadas del punto de equilibrio:
α
y=
β
θ
x=
γ

( θγ , αβ ).
El punto de equilibrio ( x 0 , y 0 ) =

La matriz A a linealizar corresponde a la matriz jacobiana:

∂f ∂f θ
0

( )( )
−β
∂x ∂y l
A= =
∂g ∂g α
γ 0
∂x ∂y β
Utilizando la ecuación secular, para calcular los valores y vectores propios para la
solución:
| A−λI |=0

Llegamos al sistema de ecuaciones lineales a resolver:

θ
λx− β y=0

{ γ
α
β
x + λy=0

donde las soluciones de λ=±

l

√ γ∝θ
l

Donde λ=± √ a1 b 1 ≡iω, esto da lugar a soluciones periódicas.

( xy (t)(t ))=C ⃗V e
1 1
i λ1 t
V2e
+C 2 ⃗
−i λ t
2

Las soluciones la podemos escribir de la siguiente forma:

x (t) cos ⁡(ωt ) sen (ωt)

( ) (
y (t )
=C1
Asen(ωt )
+C 2 ) (
− Acos( ωt) )
Para los parámetros α =0.4 , β=0.37 , γ =0.3 , δ =0.05
C 1=100 y C 2=50

(C 1 y C2 son la condisiones iniciales para la

poblacion de depredadores y de las presa)
2) Realizar el mismo procedimiento que en 1) para el péndulo simple.
Solución
 Péndulo simple:

Un péndulo simple se define como una partícula de masa m suspendida del punto
O por un hilo inextensible de longitud l y de masa despreciable.
Un péndulo simple es un ejemplo de oscilador no lineal. Se puede aproximar a un
oscilador lineal cuando su amplitud es pequeña.
Ecuaciones de movimiento:
Si la partícula se desplaza a una posición θ0 (ángulo que hace el hilo con la
vertical) y luego se suelta, el péndulo comienza a oscilar.
El péndulo describe una trayectoria circular, un arco de una circunferencia de radio
l. Estudiaremos su movimiento en la dirección tangencial y en la dirección normal.
Las fuerzas que actúan sobre la partícula de masa m son dos:
- el peso mg
- La tensión T del hilo.

Esta grafica:
Representación de las
componentes del
movimiento péndulo
simple en un determinado
sistema de referencia

La ecuación diferencial la cual la podemos obtener de un análisis del péndulo con

la segunda ley de Newton o por medio de un análisis de energías con la mecánica
lagrangiana, esta ecuación diferencial es:
g
θ̈+ sin ( θ ) =0
l
Esta es una ecuación diferencial de segundo orden no lineal, por lo tanto, es muy
complicado hallar una solución analítica. Aun así, podemos reducir el orden de
esta ecuación, para formar un sistema de dos ecuaciones lineales de primer
orden, el cual podemos solucionar por el método de los valores y vectores propios.
θ̇= y
g
ẏ + sen ( θ )=0
l

1 0
θ̇ =
()
ẏ 0 ( )(
−g
l
y
sen (θ) )
Para linealizar el problema consideramos que θ → 0entonces sen ( θ ) ≈θ, pequeñas
oscilaciones, y el sistema de ecuaciones nos quedaría de la siguiente forma:
 1 0
θ̇ =
()
ẏ 0(−g
l
() θy )
La siguiente grafica pertenece a una simulación hecha del péndulo simple
elaborada en Matlab.

3) Resolver el sistema de ecuaciones (θ̇1=w1 , θ˙2=w 2 , ẇ 1=f ( θ 1 , θ2 ) , ẇ2 =g ( θ1 , θ2 ))

por el método de runge-kutta.
a) Dadas las condiciones iniciales, mostrar una animación.
b) Mostrar graficas de teta1 y teta2 vs t, w1 y w1 vs t, teta1 vs teta2, w1 vs
w2.
x 1=l 1 sen θ1 y 1=−l 1 cos θ 1

x 2=x 1+ l2 sen θ 2 y 2= y 1−l 2 cos θ 2

Derivando con respecto al tiempo obtenemos:

ẋ 1=θ̇1 l 1 cos θ1 ẏ 1=θ̇1 l 1 sen θ1

ẋ 2= ẋ 1+ θ̇2 l 2 cos θ2

ẏ 2= ẏ 1 + θ̇2 l 2 sen θ2

Ecuaciones de movimiento
 A partir de las ecuaciones anteriores, tras realizar numerosas operaciones
algebraicas con la finalidad de encontrar las expresiones de θ̈1 , θ̈2 en
términos de θ1 , θ̇1 , θ2 , θ̇2, llegaríamos a las ecuaciones de movimiento para el
péndulo doble:

−g ( 2 m1+ m2 ) sen θ1−m2 gsen ( θ 1−2θ 2 )−2 sen ( θ 1−2 θ2 ) m2 (θ̇ 22 l 2+ θ̇21 l 1 cos ( θ1 −θ2 ) )
θ̈1=
l 1 ( 2 m1 +m2−m2 cos ( 2 θ1−2 θ2 ) )

2 sen ( θ 1−2 θ2 ) ( θ̇21 l 1 ( m1−m2 ) + g ( m1−m2 ) cos θ1 + θ̇22 l 2 m2 cos ( θ1−θ 2) )

θ̈2=
l 2 ( 2m1 +m2−m2 cos ( 2θ 1−2θ 2 ) )

Energía
La energía cinética viene expresada por:
1 1
T = m 1 ( ẋ 21 + ẏ 21 ) + m 2 ( ẋ 22+ ẏ 22 )
2 2
1 1
T = m 1 l 21 θ̇ 21+ m2 [ l 21 θ̇ 21+l 22 θ̇22 +2l 1 l 2 θ̇1 θ̇2 cos ( θ1−θ 2) ]
2 2
La energía potencial:
V =m1 g y 1+ m2 g y 2

V =−( m1+ m2 ) g l 1 cos θ 1−m 2 g l 2 cos θ 2

Por tanto, el movimiento se regirá por la lagrangiana

L=T −V
1 1
L= ( m1 +m 2 ) l 21 θ̇21 + m 2 l 22 θ̇ 22+ m2 l 1 l 2 θ̇ 1 θ̇2 cos ( θ 1−θ2 ) + ( m 1 +m 2 ) g l 1 cos θ1+ m 2 g l 2 cos θ 2
2 2
Ecuaciones de movimiento de Lagrange
Usando las ecuaciones de Lagrange en este caso particular son:
d ∂L ∂L
( ) −
dt ∂ θ̇1 ∂ θ1
=0

d ∂L ∂L
dt ( ∂ θ̇ ) ∂ θ
− =0
2 2

Calculando explícitamente las derivadas de la expresión anterior se llega a:

2
( m1 +m2 ) l 1 θ̈1 +m2 θ̈2 l2 cos ( θ 1−θ2 ) +m 2 θ̇ 2 l 2 sen ( θ1 −θ2 ) + ( m 1 +m 2 ) gsen θ1=0
m2 l 2 θ̈2 +m2 θ̈1 l 1 cos ( θ 1−θ2 ) −m2 θ̇21 l 1 sen ( θ1−θ2 ) +m2 gsen θ 2=0
 4) Obtener la versión linealizada del sistema del punto 3) y hacer lo mismo
que 3a y 3b. comparar 3 y 4.
Solución
Este sistema de EDOs no lineal, puede ser simplificado considerando que los giros
θ1 y θ2 son pequeños, por cuanto las aproximaciones cos ( θ 1−θ2 ) =1, sen ( θ1−θ 2) =0,
sen ( θ1 )=θ 1 y sen ( θ2 )=θ 2, permiten transformar el sistema mediante la linealización
indicada a continuación:

( m1 +m2 ) l 1 θ̈1 +m2 θ̈2 l2+ ( m1 +m2 ) g θ1=0

m 2 l 2 θ̈2 +m 2 θ̈1 l 1+ m2 gθ 2=0

Reformulando el problema por reducción de orden en 4 ecuaciones de primer

orden con condiciones iniciales adaptadas al nuevo sistema, pero partiendo del
sistema original.
d θ1 d θ2
=w1 ; =w 2
dt dt
d w1 [ m2 g θ2−( m1 +m2 ) g θ 1 ]
=
dt [ ( m1+ m2 ) l 1−m2 l 1 ]

d w1
d w2
=
[ −( m1 +m2) g θ1−( m1 +m2 ) l 1
dt ]
dt m2 l 2

5) Resolver el sistema de lorentz (ver diapositivas).

a) Obtener las graficas
b) Resolver la versión linealizada y comparar
SOLUCION
Modelo de Lorenz
dx
=σ ( y−x )
dt

dy
=ρx− y−xz
dt

dz
=xy−βz
dt

σ: Número de Prandtl 1
ρ: Número Rayleigh2
β: razón entre la longitud y altura del sistema

La variable x representa la velocidad y la dirección de circulación del fluido. Si x >

0 el fluido circula en sentido horario mientras que si x < 0 el fluido circula en
sentido anti-horario. La variable y representa la variación de temperatura vertical, y
z es la desviación del gradiente vertical de temperatura de la linealidad (es decir, si
z>0, los gradientes de temperatura más acusados ocurren en las fronteras).
Debiendo pertenecer x, y, z a los números reales

ρ, σ y β son constantes que dependen de las propiedades fisicoquímicas y

geométricas de la capa de fluido. Se toman los parámetros σ = 10, ρ = 28 y β =
8/3, donde el modelo de Lorenz exhibe un comportamiento caótico.

1 El número de Prandlt viene dado por la razón entre la viscosidad y la

conductividad térmica.
2 El número de Rayleight relaciona la diferencia de temperatura entre la base y el
tope del sistema, la conductividad térmica, la capacidad calorífica y la viscosidad.
Estas ecuaciones son no lineales y su solución es una ecuación imposible de
resolver con exactitud. Por consiguiente, aplicaremos métodos numéricos para
obtener soluciones aproximadas de estas funciones a partir de determinar distintos
valores iniciales para las variables.

3. ANÁLISIS DE LOS PUNTOS CRÍTICOS

Las ecuaciones de Lorenz son un modelo no lineal de 3 dimensiones con

ecuaciones diferenciales de primer orden:

F (x , y , z )=σ ( y −x )

G ( x , y , z )= ρx− y−xz

H (x , y , z)=xy−βz

Para poder describir el comportamiento de este sistema, se evaluará la existencia

de puntos críticos. Un punto crítico se define como el punto ( x*,y*,z*) donde las
tres funciones se hacen cero. Es decir,
F¿

( x*,y*,z*) es un punto crítico del sistema, entonces las funciones constantes:

x ( t )=x∗y ( t )= y∗z ( t )=z∗¿

dx dy dz
tienen derivadas, =0 , =0 , y =0 . Estas soluciones constantes son
dt dt dt
llamadas solución de equilibrio del sistema.

Para encontrar los puntos críticos de las ecuaciones de Lorenz podemos utilizar el
método de sustitución. Al igualar a cero (1), podemos despejar x,

x = -y

Sustituyendo en (3) obtenemos:

x2
z=
β

Así, reemplazando en la ecuación (2):

 x3
σx−x− =0
β

x2
x (σ −1− )=0
β

Para lo cual, si x= 0, entonces y=0 y z=0.

x2
σ −1− =0
β
Entonces

x 1,2 ¿−¿+ ¿ √ β (σ−1 ) ¿

¿

De esta forma, el sistema tiene 3 puntos críticos: (0,0,0), (

( √ β ( σ−1 ) ) , ( √ β ( σ−1 ) ) ,(σ−1)¿ y (− √ β ( σ−1 ) ) , (− √ β ( σ−1 ) ) ,( σ−1)¿

ESTUDIO DE LA ESTABILIDAD DE LOS PUNTOS CRITICOS MEDIANTE

SISTEMA LINEALIZADO
Para comprender el comportamiento del sistema es de gran importancia el análisis
de la estabilidad de los puntos críticos. Para esto, diremos que un punto crítico x =
c de una ecuación de primer orden es estable siempre que el valor inicial x0 esté
suficientemente cercano a c, entonces x(t) permanece cercano a c para todo t>0.
Un punto crítico se encuentra aislado si en la vecindad de éste no contiene otro
punto crítico. Entonces, es posible determinar las trayectorias del sistema en las
zonas cercanas a este. Para esto utilizaremos un método de linealización.
La fórmula de Taylor para funciones de 3 variables implica que si la función f(x, y,
z ) es continua y derivable cerca de un punto fijo entonces:

f (x¿ ¿ ¿+u , y¿ +v , z ¿ +u)=f ( x ¿ , y ¿ , z ¿ ) + f x ( x¿ , y ¿ , z ¿ ) + f y ( x¿ , y ¿ , z ¿ ) + f z ( x ¿ , y ¿ , z ¿ ) r (u , v , w)¿

Si se aplica la fórmula de Taylor para F, G y H, y se asume que ( x*,y*,z*) s un
punto crítico aislado tal que F(x*,y*,z*)=G(x*,y*,z*)=H(x*,y*,z*)=0, el resultado es

du
=f x ( x¿ , y ¿ , z ¿ ) + f y ( x¿ , y ¿ , z ¿ ) + f z ( x ¿ , y ¿ , z ¿ ) r (u , v , w)
dt

dv
=gx ( x ¿ , y ¿ , z ¿ )+ g y ( x ¿ , y ¿ , z¿ ) + g z ( x ¿ , y ¿ , z ¿ ) +r (u , v , w)
dt

dw
=h x ( x ¿ , y ¿ , z ¿ ) +h y ( x ¿ , y ¿ , z¿ ) +h z ( x ¿ , y ¿ , z ¿ ) +r (u , v , w)
dt

Donde los residuos r(u,v,w), s(u,v,w) y t(u,v,w) satisfacen la condición:

r (u , v , w) , s(u , v , w), t (u , v , w),
lim 2 2 2
= = =0
(u ,v , w)→(0,0,0) √ u +v +w √ u + v + w √u2 + v 2+ w2
2 2 2

y por tanto, pueden despreciarse. Convirtiéndose el sistema, en un sistema lineal:

du
=f x ( x¿ , y ¿ , z ¿ ) u+ f y ( x ¿ , y ¿ , z ¿ ) v + f z ( x ¿ , y ¿ , z¿ ) w
dt

dv
=gx ( x ¿ , y ¿ , z ¿ ) u+ g y ( x ¿ , y ¿ , z ¿ ) v +g z ( x ¿ , y ¿ , z¿ ) w
dt

dw
=h x ( x ¿ , y ¿ , z ¿ ) u+h y ( x ¿ , y ¿ , z ¿ ) v+ h z ( x ¿ , y ¿ , z ¿ ) w
dt

el cual puede escribirse como u’=Ju (siendo u=[u v w]T), cuya matriz de
coeficientes se conoce como la matriz jacobiana de las funciones F, G y H
evaluadas en el punto

F x (x¿ , y ¿ , z ¿ ) F y (x ¿ , y ¿ , z ¿ ) F z ( x¿ , y ¿ , z ¿ )

(
J ( x ¿ , y ¿ , z ¿ )= G x (x ¿ , y ¿ , z ¿ ) G y (x ¿ , y ¿ , z ¿ ) G y ( x ¿ , y ¿ , z ¿ )
H x ( x ¿ , y ¿ , z ¿ ) H y ( x ¿ , y ¿ , z ¿ ) H z (x ¿ , y ¿ , z ¿ ) )
Para las ecuaciones de Lorenz,
F (x , y , z )=σ ( y −x )

G ( x , y , z )= ρx− y−xz

H (x , y , z)=xy−βz

la matriz jacobiana resulta:

−σ δ 0

[
J ( x , yz )= (ρ−z) 1 −x
y x −β ]
Así, el sistema lineal izado resulta:

F( x , y , z ) −σ σ 0 x

[ ][ ][ ]
G( x , y , z ) ≅ ( ρ−z ) 1 −x y
H (x , y , z ) y x −β z

5. CARACTERIZACIÓN DE LOS PUNTOS CRÍTICOS

Para investigar el comportamiento de los puntos críticos de un sistema lineal se
puede utilizar el método del autovalor-autovector:

x' a b c x

[ ] [ ][ ]
y' = d e f y
z' g h y z

Donde los autovalores λ de la matriz de coeficientes A son soluciones de la

ecuación característica:

a− λ b c
det ( A−λI )=¿
[ d
g
e− λ
h y −λ ]
f =( ( a−λ ) ( e−λ )−bd ) ( y−λ )=0

5.A. PUNTO CRÍTICO (0, 0, 0)

La matriz jacobiana para el punto crítico (0,0,0) resulta:
 −σ σ 0

[
J (0,0,0)= ρ −1 0
0 0 −β ]
Aplicando el método antes descripto,

(−σ−λ) σ 0
det ( A−λI )= ρ
0 [
(−1−λ)
0
0
(−β−λ)
=0
]
Por tanto, el polinomio característico es:
( (−σ −λ )(−1−λ )− ρσ)(−β− λ)=0

Buscando las soluciones para este polinomio, si(−β−λ )=0 , entonces λ 1=−β . Y si,

(−σ −λ )(−1− λ )− ρσ=λ2 + λ ( 1+σ )+ ( σ− ρβ )=0

Entonces,
2
−(1+σ ) ± √ ( 1+σ ) −4 ( σ−ρβ )
λ 2,3=
2
Para β=8 /3, ρ=28 y σ =10 los autovalores de la matriz de coeficientes constantes
en el punto crítico (0,0,0) pertenecen a los R, son diferentes y tienen signos
opuestos. El punto ( 0,0,0 ) es un punto silla.

5. B. PUNTO CRÍTICO ( √2 β ( σ−1 ) , √2 β ( σ −1 ) , ( σ −1 ) )

Continuando el análisis para el punto ( √2 β ( σ−1 ) , √2 β ( σ −1 ) , ( σ −1 ) ) , construimos la

matriz jacobiana

−σ σ 0
2 2

[
J ( √ β ( σ−1 ), √ β ( σ −1 ) , (σ −1) )=
2
−1
√ β ( σ−1 ) 2
−1
√ β ( σ−1 )
2
−√ β ( σ −1 )
−β ]
para aplicar el método del auto-vector

(−σ−λ ) σ 0
2 2

[
det ( J ( √ β ( σ−1 ), √ β ( σ −1 ) , (σ −1) )−λI )=
2
−1
√ β ( σ −1 ) 2
2
(−1−λ ) −√ β ( σ −1 ) =0
√ β ( σ−1 ) (−β−λ ) ]
obteniendo como polinomio característico:
( (−σ −λ )(−1−λ ) +σ )(−β− λ)=0
Las soluciones para los autovalores son:
2
−(1+ σ) ± √( 1+σ ) −8 σ
λ 1=−β , λ 2,3=
2
Para β=8 /3, ρ=28 y σ =10 los autovalores para el punto crítico
( √2 β ( σ−1 ) , √2 β ( σ −1 ) , ( σ −1 ) ) resultan λ 1=−8/3 , λ 2=−2.2984 , λ 3=−4.5969
Como los tres autovalores pertenecen a los R, son distintos y con signos iguales,
el punto ( √2 β ( σ−1 ) , √2 β ( σ −1 ) , ( σ −1 ) ) es un nodo asintóticamente estable.

5. C. PUNTO CRÍTICO (− √2 β ( σ−1 ) ,−√2 β ( σ −1 ) , ( σ −1 ) )

Continuando el análisis para el punto (− √2 β ( σ−1 ) ,−√2 β ( σ −1 ) , ( σ −1 ) ) , construimos la

matriz jacobiana

−σ σ 0
2 2

[
J ( √ β ( σ−1 ), √ β ( σ −1 ) , (σ −1) )=
2
−1
2
−1
− √ β ( σ−1 ) −√ β ( σ −1 )
2
√ β ( σ −1 )
−β ]
Y como z se repite con el punto crítico antes evaluado, obtendremos el mismo
polinomio característico ( (−σ −λ )(−1−λ ) +σ )(−β− λ)=0
Cuyas soluciones como vimos son:
2
−(1−σ )± √ ( 1+ σ ) −8 σ
λ 1=−β , λ 2,3=
2
que para β=8 /3, ρ=28 y σ =10 los autovalores resultan λ 1=−8/3 , λ 2=−2.2984 ,
λ 3=−4.5969

Como los tres autovalores pertenecen a los ℝ, son distintos y con signos iguales,
el punto (− √2 β ( σ−1 ) ,−√2 β ( σ −1 ) , ( σ −1 ) ) también es un nodo asintóticamente
estable.

6. PLANOS DE FASE

Cuando el punto inicial ( x 0 , y 0 , z 0 ) no es un punto crítico, entonces la trayectoria

correspondiente es una superficie en el plano xyz, dentro del cual se mueve el
punto ( x (t) , y (t) , z ( t ) )conforme t incrementa.
El comportamiento de las soluciones de un sistema puede mostrarse
construyendo un diagrama que muestre sus puntos críticos junto con una
colección de superficies de solución típicas.

CODIGO

sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) +
a(1)*a(2)];
[t,a] = ode45(f,[0 100],[1 1 1]);% Runge-Kutta 4th/5th order ODE solver
plot3(a(:,1),a(:,2),a(:,3))
ylabel('y');
xlabel('x');
zlabel('z');
title('ECUACIÓN DE LORENTZ','FontSize',18,'FontWeight','Bold','Color',[1 0 0]);

clc;
clear all;
close all;

theta1 = 15;
theta1_ini= theta1*(pi/180);
theta1_pto_ini=0;
theta2 = 5;
theta2_ini=theta2*(pi/180);
theta2_pto_ini=0;

options = odeset('RelTol',1e-3,'AbsTol',[1e-3 1e-6]);

tf=30;
[T,Y] = ode45(@odePDobleL,[0 tf],[theta1_ini theta1_pto_ini theta2_ini
theta2_pto_ini]);

writerObj = VideoWriter('video_pendu_doble.avi');
open(writerObj);

figure1=figure('color','white');

%L1= input('Inserte el valor de L1 : ');

%L2= input('Inserte el valor de L2 : ');
%m1= input('Inserte el valor de m1 : ');
%m2= input('Inserte el valor de m2 : ');

L1=1.5;
L2=1.5;
m1=1.5;
m2=0.3;
g=9.8;
M=(m2/(m1+m2));
theta1=Y(:,1);
theta2=Y(:,3);

x1=L1*sin(theta1);
y1=-L1*cos(theta1);
x2=L1*sin(theta1)+L2*sin(theta2);
y2=-L1*cos(theta1)-L2*cos(theta2);
x=x1+x2;
y=y1+y2;
set(figure1,'Position',[80 80 700 600]);
axes1 = axes('Parent',figure1);
for t=1:100 %size(T,1);
cla(axes1,'reset');

plot([0 0],[0 -6*1],'-k','Parent',axes1);

% plot([0 -0.5],[0 0.5*1],'-k','Parent',axes2);

%
hold on;
plot([0 x1(t)],[0 y1(t)],'-b','LineWidth',1,'Parent',axes1);

hold on;
plot(x1(t),y1(t),'ok','MarkerFacecolor','r','MarkerSize',10,'Parent',axes1);
hold on;

plot([x1(t) x(t)],[y1(t) y(t)],'-b','LineWidth',1,'Parent',axes1);

hold on;
plot(x(1:t),y(1:t),'-k','MarkerFacecolor','k','MarkerSize',5);
hold on;
plot(x(t),y(t),'ok','MarkerFacecolor','r','MarkerSize',10,'Parent',axes1);

axis([-4 4 -8 2]);
text(-3.5,1.2,['t= ', num2str(t)],'FontSize',15,'FontWeight','Bold');
title('PÉNDULO DOBLE (LINEAL) ','FontSize',18,'FontWeight','Bold','Color',[1 0
0]);
xlabel('Y');
ylabel('X');
grid on
box(axes1,'on');
hold(axes1,'on');
grid(axes1,'on');

F=getframe(gcf);
writeVideo(writerObj,F);

hold off;
end

close(writerObj);
a=0.4;
b=0.37;
c=0.3;
d=0.05;
t=100;
i=1;
Z=[2.2;0.32];
f1(i)=Z(1);
f2(i)=Z(2);
x(i)=t;
h=-0.01
while t>0
i=i+1;
D=Z;
t=t+h;
F=[a*Z(1)-c*Z(1)*Z(2);d*Z(1)*D(2)-b*Z(2)];
Z=D+h*F;
f1(i)=Z(1);f2(i)=Z(2);
x(i)=t;
end
figure(1)
hold on
plot(x,f1,'b',x,f2,'r')
title('SISTEMA DEPREDADOR-PRESA','FontSize',18,'FontWeight','Bold','Color',
[1 0 0]);
figure(2)
hold on
plot(f1,f2,'g')
Z;

title('SISTEMA DEPREDADOR-PRESA','FontSize',18,'FontWeight','Bold','Color',[1
0 0]);

1. Modelo Depredador-Presa No Lineal

N=1000;
t0=0; tN=100; h=(tN-t0)/N;
a=2;b=1;c=1.2;d=0.9; %Valor de las constantes del modelo.
R(1)=3; F(1)=1;
for i = 1:N
k1=[a*R(i)-c*F(i)*R(i), -b*F(i)+d*F(i)*R(i)];
k2=[a*(R(i)+h*k1(1)/2)-c*(F(i)+h*k1(2)/2)*(R(i)
+h*k1(1)/2), -b*(F(i)+h*k1(2)/2)+d*(F(i)+h*k1(2)/2)*(R(i)
+h*k1(1)/2)];
 k3=[a*(R(i)+h*k2(1)/2)-c*(F(i)+h*k2(2)/2)*(R(i)
+h*k2(1)/2), -b*(F(i)+h*k2(2)/2)+d*(F(i)+h*k2(2)/2)*(R(i)
+h*k2(1)/2)];
k4=[a*(R(i)+h*k3(1))-c*(F(i)+h*k3(2))*(R(i)+h*k3(1)),
-b*(F(i)+h*k3(2))+d*(F(i)+h*k3(2))*(R(i)+h*k3(1))];
R(i+1)=R(i)+h*(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6;
F(i+1)=F(i)+h*(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6;
end;
T=t0:h:tN;
t=t0:h:tN;
figure(1)
hold on
plot(t,R,'g')%Población de conejos en función del tiempo
plot(t,F,'r')%Población de lobos en función del tiempo
title('Población de conejos(verde) y zorros(rojo) en función
del tiempo');
ylabel('t(años)');
hold off
figure(2)
plot(R,F,'k')
xlabel('Población de conejos');
ylabel('Población de lobos');

2. Péndulo simple
function [] = call_pend()
tspan=[0 50];
y0=[pi/4,0];
[t,y]=ode23(@pend,tspan,y0);
figure(1)
plot(t,y(:,1),'r')
xlabel('tiempo');
ylabel('theta');
title('Péndulo No lineal')

function dz = pend(t,y)
dz = zeros(2,1);
G=9.8; L=10; B=0.1; M=4;
dz(1)=y(2);
dz(2)=-B/M*y(2)-G/L*sin(y(1)); %para que sea lineal quite sin
y deje y(1)

end

end

3. Péndulo Doble
clear all;
close all;
clc;
theta1 =15;
theta1_ini=theta1*(pi/180);
theta1_pto_ini=0;
theta2 =20;
theta2_ini=theta2*(pi/180);
theta2_pto_ini=0;

options = odeset('RelTol',1e-3,'AbsTol',[1e-3 1e-6]);

tf=30;
[T,Y] = ode45(@odePendulumDoble,[0 tf],[theta1_ini
theta1_pto_ini theta2_ini theta2_pto_ini]);

writerObj = VideoWriter('video_Pendu_DobleNL.avi'); %avi

open(writerObj);

figure1=figure('color','white');

l1=1;
l2=1;
m1=1;
m2=1;
g=9.8;

theta1=Y(:,1);
theta2=Y(:,3);

x1=l1*sin(theta1);
y1=-l1*cos(theta1);
x2=l1*sin(theta1)+l2*sin(theta2);
y2=-l1*cos(theta1)-l2*cos(theta2);

x=x1+x2;
y=y1+y2;

% plot(theta1,theta2) %retrato de fases

% xlabel ('angulo1')
% ylabel ('angulo2')
%%%%% tiempo vs theta1, theta 2
% t2 = size(theta1);
% t1 = t2(:,1);
% t=0:t1-1;
%
% plot(t,theta1,'r',t,theta2,'k')
% %%%%%% t
% theta1p = Y(:,2); derivada de theta
% theta2p = Y(:,4);
% plot(theta1,theta1p,theta2,theta2p)
% plot(t,theta1p,t,theta2p)
%%%%

%VIDEO
set(figure1,'Position',[80 80 700 600]);
axes1 = axes('Parent',figure1);
for t=1:100 %size(T,1);
cla(axes1,'reset');

plot([0 0],[0 -6*1],'-k','Parent',axes1);

% DIBUJAMOS LOS OBJ A REPRESENTAR

hold on;
plot([0 x1(t)],[0 y1(t)],'-
k','LineWidth',1,'Parent',axes1);
hold on;

plot(x1(t),y1(t),'ok','MarkerFacecolor','r','MarkerSize',8,'P
arent',axes1;
hold on;
plot([x1(t) x(t)],[y1(t) y(t)],'-
k','LineWidth',1,'Parent',axes1); hold on;
plot(x(1:t),y(1:t),'-
k','MarkerFacecolor','k','MarkerSize',5); %linea
hold on;

plot(x(t),y(t),'ok','MarkerFacecolor','r','MarkerSize',8,'Par
ent',axes1);

title('PÉNDULO NO
LINEAL','FontSize',18,'FontWeight','Bold','Color',[1 0 0]);
axis([-4 4 -6 4]);
text(-4,4.2,['t= ',
num2str(t)],'FontSize',15,'FontWeight','Bold');
xlabel('x');
ylabel('x');
grid on
box(axes1,'on');
hold(axes1,'on');
grid(axes1,'on');

F=getframe(gcf);
writeVideo(writerObj,F);
 hold off;
end

close(writerObj);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%otro script
function dy = odePendulumDoble(t,y) %METODOS NUMERICOS 2019
dy = zeros(4,1);
l1=1;
l2=1;
m1=1;
m2=1;
g=9.8;

dy(1) = y(2);
dy(2) = -((dy(4))*m2*l2*cos(y(1)-y(3))
+m2*l2*(y(4)^2)*sin(y(1)-y(3))+(m1+m2)*g*sin(y(1)))/
((m1+m2)*l1);

dy(3) = y(4);
dy(4) = (m2*l1*(y(2)^2)*sin(y(1)-y(3))-
m2*l1*(dy(2))*cos(y(1)-y(3))-m2*g*sin(y(3)))/(m2*l2);

4. Lorenz

function ELorenz
clear all
t=[0 50]; % hasta 50 cuando es retrato fase 5 si son las
funciones
xinit=[0.1;0.1;0.1]; % Initial condition x=y=z=0.1
[t,x]=ode45(@lorenz,t,xinit); % Integrate in time
plot3(x(:,1),x(:,2),x(:,3) % Plot solution retrato fase
%plot(t,x(:,1),t,x(:,2),t,x(:,3)) %funciones hasta t=5
%title('Sistema de Lorenz');
%legend 'x(t)' 'y(t)' 'z(t)'
%xlabel('t');
%ylabel('x,y,z');
%zlabel('z');

function dx=lorenz(t,x)
% Parameters
sigma=10; s=28; b=8/3;
%Right hand sides
dx1=sigma*(x(2)-x(1));
dx2=s*x(1)-x(2)-x(1)*x(3);
dx3=x(1)*x(2)-b*x(3);
dx=[dx1;dx2;dx3];