MODELOS MATEMÁTICOS DE

DINÁMICA ECONÓMICA
INTRODUCCIÓN A LOS
SISTEMAS DINÁMICOS CON MATLAB.

GUÍA DE PRÁCTICAS

Última modificación: 27 de Enero de 2005

Fernando Fernández Rodríguez

Departamento de Métodos Cuantitativos en Economía y Gestión
Universidad de Las Palmas de Gran Canaria

1
CONTENIDO
1 El entorno MATLAB
1.1 Introducción de matrices en MATLAB
1.2 Algebra matricial
1.3 Matrices especiales
1.4 Rescate de elementos a partir de una matriz
1.5 Resolución de sistemas lineales

2 Funciones en MATLAB
2.1 Funciones matemáticas elementales
2.2 Funciones estadísticas elementales
2.3 La función plot
2.4 Diversas funciones de cálculo numérico (roots, poly, polyval,
fval, polyfit,quad)

3 Creación de archivos . m
3.1 Archivos “script”
3.2 Archivos “function”

4 Algoritmos
4.1 El bucle for
4.2 El bucle while
4.3 La bifurcación if

5 Sistemas dinámicos en Economía

5.1 Ecuaciones en diferencias y sistemas
5.2 Ecuaciones diferenciales
5.3 Aplicaciones a la Economía

2
INTRODUCCION DE MATRICES EN MATLAB

 1 2 3
 
Dada una matriz A =  4 5 6  puede ser escrita en MATLAB de cualquiera de las formas
7 8 9
 
alternativas:

>>A=[1 2 3;4 5 6 ;7 8 9 ]

o bien

>>A=[1 , 2 , 3;4 , 5 , 6 ;7 , 8 , 9]

Obsérvese que con la introducción de punto y coma al final de una variable matricial , esta quedará
guardada en la memoria de trabajo pero no aparecerá explícitamente en la pantalla.

Ejercicio.
Escribir las siguientes matrices:

 0 −1 3  3 
  1 3  
B =  2 − 1 2  , C =   , D = (2 3 9 − 1) , E =  π  .
 2 4 3  2 4  − e
   
El número π se escribe en MATLAB simplemente como pi, mientras que el número e se representa por
medio de exp(1).

OPERACIONES ELEMENTALES CON MATRICES

Una vez introducidas las matrices en el espacio de trabajo es posible realizar operaciones elementales
entre ellas, siempre que se verifiquen las condiciones de dimensionalidad adecuadas.
A continuación mostramos las principales operaciones entre matrices.

Suma de matrices :
>> A+B

Producto de matrices :
>> A*B

Potencia n-ésima :
>> A^n

Sumar un escalar a todos los elementos de una matriz:

>> A+2

Multiplicar un escalar por los elementos de una matriz:

>> k*A

Transpuesta:
>> A’

Determinante :
>> det(A)

3
Inversa :
>> inv(A)

División izquierda :
>> A\ B; (similar a inv(A)*B)

División derecha:
>> B/A; (similar a B*inv(A))

Producto, elemento a elemento, de dos matrices:

>> A . * B

Cociente, elemento a elemento, de dos matrices:

>> A . / B

Potencia n-ésima de los elementos de una matriz:

>> A . ^ n

Potencia de los elementos de una matriz a los exponentes dados por otra:
>> A . ^B

Cálculo de la dimensión de una matriz :

>> size(A)

Cálculo de la longitud de un vector :

>> length(V)

Valores propios o autovalores:

>> eig(A)

Vectores propios o autovectores :

>> [ V , D ]=eig(A)

Esta última instrucción da lugar a una matriz diagonal D formada por los autovalores y a una matriz V
−1
cuyas columnas son los autovectores, de modo que se verifica V XV =D.

Diagonal de una matriz :

>> diag(A)

proporciona un vector columna cuyos elementos son los de la diagonal principal de una matriz A.

Diagonal k-ésima de una matriz :

>> diag(A,k)

proporciona un vector columna formado por los elementos de la k-ésima diagonal de la matriz A (k es un
entero positivo o negativo).

Matrices diagonales :
>>diag(v)

proporciona una matriz cuadrada, diagonal, cuya diagonal principal son los elementos del vector v. Para
obtener una matriz cuadrada cuya diagonal coincida con la de determinada matriz A emplearemos la
instrucción
>> diag(diag(A))

Cuando v es un vector y k un número entero,

>> diag(v,k)

es la matriz cuadrada de orden n+abs(k) con los elementos de v en la k-ésima diagonal.

4
Matriz triangular inferior :
>> tril(A)

Elementos debajo de la k_ésima diagonal inferior de la matriz A :

>> tril(A,k)

Matriz triangular superior :

>> triu(A)

Elementos encima de la k_ésima diagonal superior de A :

>> triu(A,k)

Obsérvese que a partir de las operaciones elementales anteriores es posible calcular otros conjuntos
elementos más complejos de una matriz. Por ejemplo, la norma de las columnas de una matriz A puede
calcularse como
>> sqrt( sum ( A.^2 ) )

Ejercicio.
Hallar la norma de las filas de una matriz A haciendo uso de las operaciones anteriores.

Ejercicio.
Calcular las matrices que se especifican a continuación a partir de la A, B, C, D.

3 1   2 − 1
   2 − 2  
A=  0 2  B=   C=  − 1 3  D= (2 1)
 2 −1 1 3  0 1 
  

a) DB 2
b) BC T
T
c) (CB) D
T
d) det( AC )
T −1
e) ( A C )

Ejercicio.
Dadas las matrices A y B calcular las matrices que se especifican en los diferentes apartados

1 3 2 3 1 1
   
A = 0 1 4 B = 1 − 2 0
1 − 2 0 2 1 1
   

a) A + B T , AT B , det( A 2 + B −1 )
b) 2 A + ( A + B )
2

−1 −1
c) A + 3I + B .
d) Autovalores de A+B
e) Autovectores de B
f) Matriz cuyos elementos se obtienen multiplicando los de A y B, elemento a elemento, es decir,
( ) (
cij = aij bij . )
g) Matriz cuyos elementos se obtienen dividiendo los de A y B, elemento a elemento, es decir,
( ) (
cij = aij / bij . )

5
h) Matriz cuyos elementos se obtienen multiplicando los cuadrados de los elementos A y por los
cuadrados de los elementos de B, elemento a elemento, es decir (c ) = (a
ij ij
2 2
bij . )
i) Matriz cuyos elementos se obtienen elevando los elementos de A a los exponentes de la matriz B,
( )
elemento a elemento, es decir (cij) = aij ij .
b

OPERACIONES ELEMENTALES CON VECTORES

Si a y b son vectores de la misma dimensión pueden definirse, de forma específica, las siguientes
operaciones entre vectores:

a. * b = (a1b1 a2b2 .......anbn ) (Producto elemento a elemento)

a * b' = a1b1 + a 2 b2 + ....... + a n bn (Producto escalar)

 a1   a1b1 .... a1bn 

   
a '*b =  .... (b1 .... bn ) =  ... .... ...  (Producto de Kroneker)
a   a b .... a b 
 n  n 1 n n 

MATRICES ESPECIALES
El MATLAB es capaz de generar numerosas matrices especiales con gran facilidad, como se muestra a
continuación:

Números enteros entre 1 y 20:

>> X1=1:20

Números enteros entre 1 y 200:

>> X2=(1:200)

Números entre 1 y 5 tomados de 0.1 en 0.1:

>> X3=(1:0.1:5)

Matriz 3x4 de números aleatorios gaussianos:

>> X5=randn(3,4)

Matriz 2x10 de números aleatorios distribuidos uniformemente:

>> X6=rand(2,10)

Matriz identidad 4x4:

>> X7=eye(4)

Matriz de ceros 3x3:

>> X8=zeros(3)

Matriz de ceros de orden 3x4:

>> X9=zeros(3,4)

Matriz de unos de orden 4x4:

>> X10=ones(4)

6
Matriz de unos de orden 4x2:
>> X11=ones(4,2)

Igualmente, siendo A la matriz definida inicialmente, MATLAB es capaz de generar otras matrices del
tamaño de la matriz A en cualquiera de las siguientes formas:

Matriz de números aleatorios gaussianos de las dimensiones de A:

>> X11=randn(size(A))

Matriz de números aleatorios uniformes de las dimensiones de A:

>> X12=rand(size(A))

Matriz identidad de las dimensiones de A:

>> X12=eye(size(A))

Matriz de unos de las dimensiones de A:

>> X13=ones(size(A))

Matriz de ceros de las dimensiones de A:

>> X14=zeros(size(A))

Ejemplo
Basándonos en el generador randn de números aleatorios Gaussianos N(0,1), es posible muestras
aleatorias de cualquier distribución normal. Por ejemplo, para generar números aleatorios muestreados de
una normal N ( µ , σ 2 ) bastará con realizar sucesivamente la operación µ + σ * randn , es decir,

>> mu+sigma*randn

donde mu = µ y sigma = σ .

Ejercicio.
Escribir una sucesión de 1000 números en los intervalos que se citan a continuación. Empleando la
función hist(x) , que veremos posteriormente, hallar sus correspondientes histogramas.
a) Numeros aleatorios uniformemente distribuidos entre –2 y 2.
b) Números aleatorios uniformemente distribuidos entre 0 y 5.
c) Números aleatorios gaussianos con media 2 y varianza 3.
d) Números aleatorios gaussianos con media 1 y varianza 0.1.

Ejemplo
Veamos una ilustración del teorema central del límite: Si tipificamos una suma de distribuciones
idénticas obtenemos, asintóticamente, una distribución normal. Para ello procedamos a sumar
distribuciones uniformes en el intervalo [0,1], tipifiquemos, y hallemos su histograma.

RESCATE DE DETERMINADOS ELEMENTOS A PARTIR DE

UNA MATRIZ
Dada una determinada matriz A = ( aij ) es posible obtener, de forma sistemática, numerosos
subconjuntos de sus elementos :

>> A(i, j) proporciona el elemento (i,j) de la matriz A.

>> A(i, :) proporciona la fila i.

>> A(:, j ) proporciona la columna j.

7
Si X es un determinado vector e I un conjunto de índices, X (I ) proporciona un nuevo vector que
contiene los elementos del vector X correspondientes a los índices I.
Por ejemplo, sea X=randn(1,30) y sea I=(5:2:20); X ( I ) proporciona los elementos de X comprendidos
entre el 5 y el 20 tomados de dos en dos.

Igualmente pueden considerarse subconjuntos de elementos de un vector X que verifiquen alguna relación
de desigualdad:

>> Y = X ( X <= 1);

da lugar a un vector Y formado con los elementos menores o iguales que uno del vector X.

Si deseamos el conjunto de índices de los elementos de un vector que verifican determinada relación
acudiremos al comando:

>> I=find(X<1)

que proporciona el conjunto de índices donde la variable X toma valores menores que 1.

También es posible eliminar los elementos de un vector que cumplan determinada relación:

>> X ( X < 0) = [ ];
hace desaparecer los elementos negativos del vector X.

Ejercicio:
1. Dado el vector X=rand(20,1), comprobar que find ( X < 2 ) = X ( X < 2 ).
2. Hallar el conjunto de índices donde X<0.

Igualmente se pueden crear variables binarias que toman los valores 0 ó 1 según sea verdadera o falsa
determinada relación.

>> Y=X<4;

Crea un vector Y que toma el valor 1 donde el correspondiente elemento de X sea menos que 4 y 0 en
caso contrario.

Por último, es posible cambiar el valor de determinados elementos de una matriz.

Si, por ejemplo, tenemos que
 1 2 3
 
A =  4 5 6
7 8 9
 
y hacemos

>> A(1,3)=0.5

el elemento A(1,3) de la matriz tomará el valor 0.5.

Si hacemos

>> A(:,3)=[2 3 –1]'

sustituiremos la tercera columna de la matriz A por (2 3 –1)'. Análogamente

>> A(1,:)=[1 0 1]

sustituye la primera fila de la matriz A por (1 0 1).

8
RESOLUCION DE SISTEMA LINEALES DE ECUACIONES
Un sistema de ecuaciones lineales puede ser resuelto con MATLAB de varias formas.
→ → → →
Si escribimos el sistema en la forma matricial A X = b , donde A es una matriz nxn y X y b son
vectores columna n-dimensionales, la solución puede obtenerse de dos formas alternativas:
i) Empleando la matriz inversa de A:
→ →
−1
>> X = A b .
ii) Empleando un método más eficiente, de eliminación gaussiana, por medio de la
división izquierda:
→ →
>> X = A \ b .

Por ejemplo, para resolver el sistema

2x-y=1   2 −1  x   1 
 , es decir,    =  
x+y=2  1 1  y   2
procederemos de la forma,

>> inv ( [2 -1;1 1] ) * [1;2]

ans = 1 , 1

→ →
Si escribimos el sistema en la forma matricial X A = b , donde A es una matriz nxn y
→ →
X y b representan ahora vectores fila n-dimensionales, la solución puede obtenerse de dos formas
alternativas:

iii) Empleando la matriz inversa de A:

→ →
−1
>> X = b A .
iv) Empleando un método más eficiente, de eliminación gaussiana, por medio de la
división derecha :
→ →
>> X = b / A .

Ejercicio.
Hallar las soluciones de los siguientes conjuntos de ecuaciones lineales. Describir cada conjunto
→ → → →
empleando las dos formas de notación matricial A X = b y X A = b . Comparar los resultados
obtenidos según empleemos en el método de resolución la matriz inversa, la división izquierda o la
división derecha.

x+3y=1 
a) 
2x-y=3 
3x+2y-z=10 

b) -x+3y+2z=5 
x-y-z=-1 


9
 x+y+z+t=4 
2x-y+t=2 

c) 
x-2y-z-t=2 
x-2y-3z+t=-3

SAVE y LOAD: GUARDAR Y RESCATAR RESULTADOS

La instrucción save, junto con sus diversas variantes, permite salvar los resultados obtenidos.

Si A es una matriz, save A guarda, en el directorio desde donde ha sido llamado el MATLAB, una matriz
denominada A.mat.

Se pueden guardar igualmente diversas matrices previamente generadas:

save datos A B C;

guarda las matrices A, B y C en un archivo denominado datos.mat.

Es posible, igualmente, guardar ficheros en formato ascii:

save dastos.dat A -ascii -double

de modo que la matriz A será guardada en formato ascii con el nombre de datos.dat, en doble densidad.

La instrucción load permite cargar al espacio de trabajo los ficheros previamente guardados en un
determinado directorio, siempre que dicho directorio se encuentre abierto:

load a y load datos

permitirán recuperar los archivos a.mat y datos.dat, respectivamente.

10
ALGUNAS OPERACIONES DE INTERFACE CON EL ESPACIO
DE TRABAJO
La instrucción who ofrece un listado de las variables definidas en el espacio de trabajo.

La instrucción whos ofrece un listado de las variables definidas en el espacio de trabajo que incluye
información adicional sobre su tamaño y otras características.

La instrucción clear borra todas las variables previamente creadas en el espacio de trabajo.

El formato numérico que MATLAB muestra por defecto es un formato de cuatro cifras decimales.
Existen otros formatos numéricos disponibles. El más importante es un formato numérico de dieciséis
cifras decimales que es el número de cifras con el que MATLAB realiza las operaciones de cálculo. Para
acceder a este formato es suficiente con teclear en la pantalla format long. Para regresar al formato
habitual debe teclearse format short.

La instrucción clc limpia la pantalla.

FUNCIONES MATEMATICAS ELEMENTALES

El MATLAB contiene numerosas funciones matemáticas elementales, algunas de las cuales se numeran a
continuación.

Trigonométricas

sin : seno.
asin : arco seno.

cos : coseno.
acos : arco coseno.

tan : tangente.
atan : arco tangente.

Exponenciales

exp : exponencial .

log : logaritmo natural.

log10 : logaritmo decimal.

sqrt : raiz cuadrada.

Complejas

abs : valor absoluto.

angle : ángulo de fase.

conj : complejo conjugado.

imag : parte imaginaria.

real : parte real.

11
Numéricas

fix : parte entera.

sign : función signo.

La característica fundamental de cada una de estas funciones es que no solamente son aplicables a un
número real sino a un vector o a una matriz. En tal caso la función actuará sobre cada uno de los
elementos de la matriz o del vector.
 a11 a12 a13 
 
Por ejemplo, sea X =  a 21 a 22 a 23  . Si calculamos el seno de la matriz X, el resultado será:
a a 32 a 33 
 31

 sin( a11 ) sin( a12 ) sin( a13 ) 

 
sin( X ) =  sin( a 21 ) sin( a 22 ) sin( a 23 )  .
 sin( a ) sin( a ) sin( a ) 
 31 32 33 

Ejercicio.
Realizar las siguientes operaciones entre números reales.

5 3
a) 1+ +
e π2
23
b) + 2− 5
e
15
c)
1 e sen(π / 8)
+ + + ln(23)
3 2π 4
π + e2 + 2
d)
sen(3π / 4) + log10 (14)
cos(π 2 ) + tg(π / 5) 5
e) ( )
ln( 3 + 2π )
120 1
f) 3 π2 − +
π − 14 π

FUNCIONES ESTADISTICAS ELEMENTALES

El MATLAB dispone de diversas funciones estadísticas elementales. Si éstas se aplican a un vector X
darán lugar al correspondiente estadístico sobre sus observaciones. Si X es una matriz estas funciones
actúan sobre cada una de las columnas proporcionando un vector de resultados, que corresponde al valor
del estadístico aplicado sobre cada uno de los vectores columna de dicha matriz.

Media aritmética:
>> mean(X)

Desviación estándar:
>> std(X)

12
Mediana:
>> median(X)

Producto:
>> prod(X)

Suma:
>> sum(X)

Máximo:
>> max(X)

Mínimo:
>> min(X)

Otras funciones estadísticas importantes son cov y corrcoef que tienen una forma similar de actuación.

Si X es un vector

>> cov(X)

da lugar a la varianza del vector X.

Si X es una matriz
>> cov(X)

da lugar a la matriz de covarianzas entre las columnas.

Es inmediato que

>> diag(cov(X))

proporciona un vector de varianzas de cada columna, mientras que

>> sqrt( diag ( cov ( X ) ) )

proporciona un vector de desviaciones standard entra las columnas.

Cuando X e Y son vectores columna

>> cov(X,Y)

es igual a

>> cov([X Y]).

Resulta, también, de importancia el operador primeras diferencias diff.

Para un vector X

diff(X) = [X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)].

Para una matriz X diff(X) es la matriz de diferencia entre las columnas

diff(X) =[X (2 : n , : ) – X ( 1: n-1 , : ) ]

Por último,

>> diff(X,n)

es el operador de diferencias n-ésimas.

La confección de histogramas se realiza mediante la instrucción

13
>> hist(X)

donde se supone que X es un vector numérico que representa una determinada muestra. Si queremos
precisar el número k de intervalos de clase que se consideran en el histograma, utilizaremos la instrucción

>> hist(X,k).

Para reordenar, de manera ascendente, los elementos de un determinado vector X emplearemos la función

>> sort(X)

Esta función permite igualmente recuperar el orden en que estaban situados inicialmente los elementos
en un determinado conjunto de índices I, de la forma:

>> [ Y , I ]=sort(X)

Ejercicio.
Para comprobar el funcionamiento de sort, realizar el siguiente conjunto de instrucciones:

>> X = rand(1:10);
>> sort(X');
>> [ Y , I ] = sort(X');
>> [ Y’ I’ X' ]

Obsérvese que Y=X( I ).

LA FUNCIÓN PLOT
Dados dos vectores X e Y de la misma longitud,

>> plot(X,Y)

dibuja el vector X frente al vector Y.

Cuando X es un vector

>> plot(X)

dibuja los elementos de X frente al conjunto de sus índices.

Cuando X es una matriz

>> plot(X)

dibuja cada una de sus columnas frente al conjunto de sus índices.

Por defecto la función plot dibuja uniendo los puntos de la gráfica en forma poligonal. Existen otras
posibilidades como las siguientes:

>> plot(X , ' . ' ) que dibuja solamente los puntos de la gráfica.

>> plot(X , ' + ' ) que dibuja cruces en lugar de puntos en la gráfica.

>> plot(X , 'r - - ' ) que dibuja una gráfica en rojo con trazo discontinuo.

MATLAB permite dibujar simultáneamente varios vectores de la misma longitud pero con trazos y color
diferenciados. Para ello hay que especificar en cada caso el conjunto de índices, el vector y el formato de
color y unión entre puntos, tal como se muestra en el siguiente ejemplo:

>> X=rand(1,50); Y=rand(1,50);

14
>> plot( (1,50) , X, ' r - ' , (1,50) , Y, ' b - - ' )

Si escribimos previamente a una instrucción plot la instrucción grid la gráfica quedará dibujada dentro de
una rejilla.

Ejercicio.
Mírese en el sistema de ayudas del MATLAB las diferentes posibilidades de la función plot:

>> help plot

Ejemplo.
sen( x) 1
Dibujar las gráficas de las funciones y = sen( x) , y = x sen( x ) , y = , y = sen( ) .
x x
MATLAB dibuja la gráfica de una función obteniendo los valores que toma sobre un determinado
conjunto de índices. El primer paso, por lo tanto, para dibujar una gráfica es definir el conjunto de
índices. Para estas funciones trigonométricas resulta adecuado definir el siguiente conjunto de índices:

>> a = ( 1 : 200 ) ; b = a / pi ;

las funciones anteriormente señaladas pueden ahora dibujarse simplemente con

>> plot ( sin ( b) )

>> plot (b.*sin ( b) )

>> plot (sin(b). / b)

>> plot (sin(1. / b) )

Ejemplo.
Dibujar la gráfica del polinomio y = x − 3x + 2 x − 1 .
3 2

>> x=-10:0.1:10;
>> y=x.^3-3*x.^2+2*x-1;
>>plot(x,y);

Ejercicios
x2 x2
• Dibujar la gráfica de las funciones y = , e y = en el intervalo [−10,10] .
1 − x2 1 + x2
• Dibujar la gráfica de la función y = tan(sen( x ) − sen(tan( x )) en el intervalo [− π , π ] .
Dibujar conjuntamente las gráficas de las funciones y = x − 3x + 2 x − 1 e y = sen (x ) .
3 2
•

MATLAB posee una función inteligente de representación gráfica denominada fplot.m.

fplot tiene como inputs el nombre de la función en una cadena de caracteres y el rango en que se desea
pintar la función.
Por ejemplo fplot ( ' sin ' , [0 pi] ) .

Otro ejemplo más sofisticado sería.

>> f='2*exp(-x). * sin (x)';
>> fplot ( f , [0 8] )
>> title ('seno') , xlabel( 'eje x' ) , ylabel( 'eje y' )

Ejercicio

15
 ax0
Dibujar la curva de crecimiento logístico x(t ) = , que corresponde a las
bx0 + (a − bx0 )e− (t −t0 )
dx
soluciones de la ecuación diferencial = ax − bx 2 , para diversos valores de los parámetros
dt
a, b, x0 , t0 . Dibujar conjuntamente varias soluciones.

DIVERSAS FUNCIONES DE CALCULO NUMÉRICO

MATLAB posee numerosas funciones específicas que permiten realizar complejas operaciones de cálculo
numérico. A continuación se muestran algunas de ellas.

roots :

Permite calcular las raíces de un polinomio.

roots(X) es un vector fila formado por las raíces de un polinomio cuyos coeficientes son los elementos
del vector X.
Por ejemplo, para calcular las raíces del polinomio x 3 + 3x 2 − 2 x + 5 , formaremos previamente el
vector

>> X=[1 3 –2 5];

las raíces se obtendrán de la forma:

>> roots(X) = -3.8552 0.4276 + 1.0555i 0.4276 - 1.0555i

Ejercicio
Hallar la tasa interna de rendimientos (TIR) de una inversión (un bono) cuyo precio es de 110 € y cuyo
valor facial es de 100 € y tiene tres periodos hasta la amortización, produciendo unos pagos de cupón de
10€ durante sus tres años de vida.
El precio de dicho bono viene dado por

c c c+ A
p= + +
1 + r (1 + r ) (1 + r )3
2

La tasa interna de rendimiento se obtendrá resolviendo la ecuación

p (1 + r )3 − c(1 + r ) 2 − c(1 + r ) − c − A = 0

es decir, 110(1 + r ) − 10(1 + r ) 2 − 10(1 + r ) − 110 = 0

3

>> roots ( [ 110 -10 -10 -100 ] )

obtiene las raices 1.0322 , -0.4707 + 0.8119i , -0.4707 - 0.8119i ,

luego 1 + r = 1.0322 , r = 0.0322 .

fzero :
Busca las raíces de cualquier función de una variable en torno a un determinado punto x0. Su sintaxis es
fzero( 'función' , x0 ) .

16
Ejemplos
Hallar el cero de la función y = sen( x) más cercano al punto x=2.
fzero( ' sin ', 2 )=3.1416
Cuando más adelante aprendamos a construir funciones arbitrarias en MATLAB podremos, por ejemplo,
1 1
hallar un cero de la función y = + − 6 próximo al valor 1.4.
( x − 0.3) + 0.01 ( x − 0.9) 2 + 0.04
2

>> x_cero=fzero ( ' funcion ' , 1.4 );

>>x_cero=1.2995
>>y_cero=3.5527e-15 % valor obtenido cuando evaluemos el valor de la función en la raíz.

poly :

Proporciona el polinomio característico de una matriz cuadrada.

Si A es una matriz nxn poly(A) es un vector fila de dimensión n+1 cuyos elementos son los coeficientes
del polinomio característico.

Ejemplo.
1 4 − 2
 
Sea la matriz cuadrada A =  4 0 7  . Calcularemos su polinomio característico por aplicación
9 −1 5 
 
de la instrucción poly, resultando :

>> poly(A) = 1 –6 14 -187

Por tanto, el polinomio característico será λ3 − 6λ2 + 14λ − 187 .

Cuando V es un vector, poly(V) es un vector cuyos elementos son los coeficientes del polinomio que
tiene por raíces los elementos de V.

Ejemplo: Construir el polinomio que tiene por raíces los números 2 , 5, 3 , 8 .

Para ello formamos el vector

>> V=[ 2 5 3 8 ];

El polinomio buscado será λ4 − 18λ3 + 111λ2 − 278λ + 240 , según se obtiene mediante la
instrucción

>> poly(V)= 1 -18 111 -278 240

Ejercicio.
Encontrar las raíces del polinomio cuyos coeficientes son los números naturales comprendidos entre 1 y
20.

Ejercicio.
1 4 − 2
 
Dada la matriz A =  4 0 7  calcular sus autovalores utilizando las instrucciones poly y roots.
9 −1 5 
 

eig:

17
Calcula los valores propios o autovalores de una matriz A.
2 1 1 
 
Sea la matriz A =  2 3 4 
 −1 −1 − 2
 
>> eig(A) = 3 1 -1

Igualmente es posible obtener los vectores propios o autovectores guardando las salidas de eig(A) en dos
nuevas matrices D y V de la forma :

>> [ V , D ]=eig(A) .

Esta instrucción da lugar a una matriz diagonal D formada por los autovalores y a una matriz V cuyas
−1
columnas son los autovectores, de modo que se verifica V X V = D . Concretamente:
>> [ V , D ]=eig(A)

 − 0.5345 − 0.7071 0.0000 

 
>> V =  − 0.8018 0.7071 − 0.7071
 0.2673 0.0000 0.7071 

3 0 0 
 
>> D =  0 1 0 
 0 0 − 1
 

Obsérvese que roots ( poly ( A ) ) = eig (A) .

polyval:
Permite evaluar el valor numérico de un polinomio.
Si V es un vector que contiene los coeficientes de un polinomio la función polyval(V, s) proporciona el
valor numérico de dicho polinomio en el número s.

Ejemplo.
Dado el polinomio p ( x ) = x − 4 x − 2 x + 1 , calcular el valor numérico de dicho polinomio en 5.
3 2

Para ello formamos el vector

>> V=[ 1 -4 -2 1];

el valor numérico p(5) se calcula simplemente como

>> polyval(V,5) = 16

Si S es una matriz o vector el polinomio se evalúa en todos los números que forman S.

Ejercicio.
1 4 − 2
 
Consideremos una matriz cuadrada A, como por ejemplo A =  4 0 7  . Calcular su polinomio
9 −1 5 
 
característico y hacer uso de la instrucción polyval para verificar el teorema de Cayley-Hamilton con la
matriz A: Toda matriz es raíz de su polinomio característico.

18
feval :
Evalúa el valor numérico de una función en un punto x0 . Su sintaxis es la siguiente

>>feval ( ' función ' , x0 )

Ejemplo.
>> feval ( ' sin ', pi / 2 ) = 1

Observación:
En general, x0 puede ser un vector multidimensional.

Polyfit :

Permite realizar ajuste polinomial a una nube de puntos por el método de los mínimos cuadráticos.
Dados dos vectores numéricos x e y de la misma longitud

>> polyfit(x,y,n)

encuentra los coeficientes de un polinomio p(x) de grado n que ajusta los datos en el sentido de los
mínimos cuadrados.

Ejemplo.
Vamos a construir, primeramente, una nube de puntos (x,y) perturbando aleatoriamente la función
polinómica y = x :
2

>> x=0:0.1:3;
>> y=x.^2 + 2*randn ( size ( x ) ) ;
>> plot ( x , y )

Podemos ahora realizar diversos ajustes polinomiales, como por ejemplo:

>> polyfit(x,y,1) = 31 -165.3333


que daría lugar a la resta de regresión y1 = 31x − 165.3 .

Permitiendo dibujar, simultáneamente, la nube de puntos y la recta de regresión mediante

>> plot ( x ,y ,’ . ’ , x , y1 , ’ - ‘ )

De forma análoga,

>> polyfit(x,y,3) = 0 -0.0001 0.0738 2.8073

daría lugar al polinomio de regresión y = −0.0001x + 0.0738 x + 2.8073

2

Ejercicio.
Consideremos la nube de puntos (x,y) dada por

>> x=-2 : 0.1 : 2

mientras que las abscisas vienen dadas por yt = xt − 2 xt + xt + 1 + ε t , siendo ε t una familia de
3 2

variables aleatorias independientes y N( 0 , 1 ) para cada instante t.

1. Hallar los ajustes polinómicos de grados 1, 2 y 3 en el sentido de los mínimos cuadrados.
2. Calcular, en cada caso, el error cuadrático medio que se produce en cada uno de los ajustes
realizados.

19
 3. Dibujar en la misma gráfica la nuve de puntos y los ajustes polinomiales del apartado 1.

quad:

Permite realizar integración numérica (el término quad viene de cuadratura, sinónimo con que se designa
a la integral).

>> quad( 'nombre_de_función' , a , b )

calcula el valor aproximado de la integral de la función f(x) entre a y b. Para realizar el cálculo
aproximado MATLAB realiza una partición adecuada del intervalo de integración [a,b] y utiliza la regla
de Simpson. Dentro de las comillas ' ' debemos colocar una función previamente definida por el
MATLAB o que nosotros creemos en un fichero .m, como veremos más adelante.

Si deseamos emplear como procedimiento de integración aproximada el método de Newton-Cotes,

MATLAB dispone del comando

>> quad8( 'nombre_de_función ', a , b ) .

Ejemplo.
3
Calcular ∫ sen( x) dx .
0
Como el MATLAB dispone de una función sin que representa la función seno, la integral se calculará
como :

>> quad( ' sin ' , 1 , 3 ) = 1.5303

Ejercicio.
π
Obténgase el valor de la integral ∫ sen( x) dx por tres procedimientos de integración numérica:
0
• Regla de Simpson empleando quad( ' sin ', 0 , pi ) .
• Método de Newton-Cotes empleando quad8( ' sin ' , 0 , pi ) .
• Un burdo procedimiento de integración aproximada consistente en dividir el intervalo de integración
[ 0 , pi ] en segmentos de longitud 0.1, formar sobre cada uno de ellos un rectángulo y finalmente
sumar el área de todos los rectángulos, es decir : x = 0 : 0.1: pi ;
π

∫ sen( x) dx ≅ sum(sin( x)) * 0.1 .

0
• Comparar cada uno de los resultados anteriores con el valor exacto de dicha integral:

π π

∫ sen( x) dx = [− cos( x)]

0
= 2.
0
Emplear formato numérico de dieciséis cifras ( format long ) para realizar las comparaciones entre las
diferentes técnicas.

20
DERIVACIÓN
La derivación numérica es una operación mucho más inestable que la integración numérica. Mientras que
la integración numérica da cuenta de una propiedad macroscópica o global, la derivación describe una
operación microscópica o local. MATLAB puede tratar el tema de la derivación de diversas formas.

• Cuando estamos en presencia de polinomios la función derivada puede encontrarse utilizando la

función polyder.m .

Ejemplo.
Dado el polinomio y = x − 5 x + 2 x + 1 , con el fin de calcular su derivada formamos el vector de
3 2

coeficientes y aplicamos el comando polyder.

>> p=[ 1 , -5 , 2, 1 ]
>> polyder ( p ) =[ 3 , -10 , 2 ]

que son los coeficientes del polinomio derivada.

• Cuando estamos en presencia de una función definida numéricamente

( x1 , y1 )
( x2 , y2 )
................
( xn , yn )

la derivada se aproxima por medio por medio de la función diff.m .

Si x = ( x1 , x 2 ,....., x n ) , la función diff actúa de la forma diff( x ) = ( x 2 − x1 ,......, x n − x n −1 ) .

La derivada puede obtenerse y dibujarse así:

>> dy = diff ( y ) . / diff( x ) ;

>> xx = x ( 1 : length ( x ) – 1 ) ;
>> plot ( xx , dy )

• Cuando se trata de una función arbitraria y = f(x ) la derivada en un punto x0 puede

aproximarse numéricamente como la media de las derivadas laterales:

+ f ( x 0 + h) − f( x 0 ) − f ( x 0 ) − f ( x 0 − h)
Siendo f ' ( x 0 ) = lim , f ' ( x 0 ) = lim
h →0 h h →0 h

f ' ( x0+ ) + f ' ( x0− ) f( x0 + h) − f( x0 − h)

resulta que f ( x 0 ) = ≅
'

2 2h

• MATLAB dispone igualmente de una toolbox de matemática simbólica que permite la

determinación de la función derivada.

fmin , fmins :
Calcula los mínimos relativos de funciones de una variable. La sintaxis es la siguiente:

>> fmin ( ' función ' , x1 , x2 )

donde busca el mínimo de la función en el intervalo [ x1 , x2 ].

21
Con el formato

>>fmin ( ' función ' , x1 , x2 , options)

utiliza un vector de parámetros de control.

Ejercicio
Consultar el menú de ayuda fmin haciendo

>> help fmin

Para encontrar un máximo relativo es suficiente observar que max {f( x )} = min {− f( x )} .

Ejemplo.
−x
Encontrar los máximos y mínimos de la función y = 2e sen( x ) . Podemos proceder como sigue:

>> f=' 2 * exp (-x) * sin(x) ';

>> x_min=fmin ( f , 0 , 8 ); (3.9270)
>> y_min =feval ( f , x_min ); (-0.0279)
>> g='- 2*exp(-x)*sin(x)';
>> x_max=fmin ( g , 0 , 8 ); (0.7854)
>> y_max =feval ( g , x_max ); ( 0.6448)

De forma similar fmins calcula mínimos de funciones de varias variables en el entorno de un punto x0.

La sintaxis es ahora:

>> fmins ( ' función ' , x0)

y también

>> fmins ( ' función ' , x0, options) .

Ejemplo.
Calcular el mínimo local de la función f( x1 , x 2 ) = 100( x 2 − x1 ) + (1 − x1 ) .
2 2 2

Para ello construimos el archivo.m

% banana.m

function y=banana(x)

y=100*( x (2) – x (1). ^2). ^2 + ( 1 – x (1) ). ^2 ;

Para obtener el mínimo es suficiente escribir

>> fmins ( 'banana' , [ -1.2 , 1] );

Resultando

>> x = ( 1 , 1 ).

22
CREACION DE UN GUIÓN EN FORMA DE ARCHIVO .M
Una de las capacidades más importantes de MATLAB es la de poder escribir una sucesión de
instrucciones en un determinado archivo que pueden, posteriormente, ser ejecutadas al llamar a ese
archivo. Tal sucesión de instrucciones juega un papel similar al de un guión de una película
cinematográfica y se almacenarán en un archivo .m del bloc de notas. Si estos archivos son
conjuntamente almacenados dentro de una determinada carpeta para poder ejecutarlos es preciso que
desde MATLAB este abierta dicha carpeta.

Ejemplo.
Crear en el bloc de notas un archivo con el nombre de grado2.m que calcule las soluciones de una
ecuación de segundo grado.

Procederemos escribiendo en un archivo que guardaremos con el nombre de grado2.m el siguiente

conjunto de instrucciones:

%calcula las raíces de una ecuación de segundo grado

a=1; b=2; c=3;

x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);

[x1,x2];

Obsérvese que se pueden añadir comentarios dentro de un archivo .m siempre que vayan precedidos de
la instrucción % .

Para cambiar el valor de los coeficientes a b y c es preciso introducir nuevos valores numéricos dentro de
grado2.m .
De forma alternativa MATLAB permite la introducción de estos valores por medio de la instrucción
input que producirá al llamar el programa un mensaje en la pantalla demandando los valores de a, b y c.
Concretamente el programa anterior podría escribirse en la forma:

% calcula las raíces de una ecuación de segundo grado

a = input(‘coeficiente del término x^2 = ‘);

b = input(‘coeficiente del término en x = ‘);

c = input(‘coeficiente del término independiente = ‘);

x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);

[x1,x2];

Alternativamente los valores numéricos de los coeficientes a, b y c pueden introducirse simplemente

escribiendo su valor en la ventana de comandos. Por ejemplo escribiendo

>> a=3 ; b=5 ; c=-2 ;

23
también somos capaces de ejecutar la siguiente versión del programa grado2.m , donde los valores de a,
b y c se introducen desde la ventana de comandos, sin tener que manipular el programa en el bloc de
notas:

% calcula las raíces de una ecuación de segundo grado

x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);

[x1,x2];

Ejemplo.
Escribir un archivo .m que dibuje la gráfica de un polinomio de tercer grado.
Crearemos un archivo llamado, por ejemplo, graf1.m de la forma

% gráfica de un polinomio de grado tres

% y=ax^3+bx^2+cx+d

a=input('coeficiente a = ');
b=input('coeficiente b = ');
c=input('coeficiente c = ');
d=input('coeficiente d = ');

x=-10:0.1:10;
y=a*x.^3+b*x.^2+c*x+d;

plot(x,y)

Ejercicio
Construir un archivo.m con el fin de dibujar la curva de crecimiento logístico
ax0
x(t ) = , que corresponde a las soluciones de la ecuación diferencial
bx0 + (a − bx0 )e− (t −t0 )
dx
= ax − bx 2 , para diversos valores de los parámetros a, b, x0 , t0 .
dt

CREACION DE NUEVAS FUNCIONES

El MATLAB permite la creación de nuevas funciones. Tales funciones deben ser escritas en archivos .m
que llevarán idéntico nombre que la función que deseamos crear. La estructura sintáctica de un
archivo.m, para una función, contiene en la primera línea, como primera instrucción, la palabra function
seguido de la expresión y= f ( x ) , donde f(x) puede tener el nombre alusivo a la nueva función que
deseamos crear. Una característica deseable de las funciones que vamos a crear, con el fin de que posean
las mismas propiedades que las funciones originarias del MATLAB, es que puedan actuar sobre una
variable matricial. Tal es el caso del ejemplo que citamos a continuación.

Ejemplo.
Definir, en el MATLAB, la función cúbica y = x + 2 x − x + 4 .
3 2

Construiremos un archivo llamado cubica.m que tendrá la siguiente estructura

24
function y=cubica(x)

% define una función cubica

a=1; b=2; c=-1; d=4;

y=a*x.^3+b*x.^2+c*x+d;

Una característica esencial de las funciones es que sus variables son locales, lo que quiere decir que no
dejan ninguna señal en el espacio de trabajo después de ser ejecutadas. Como consecuencia de esto no
será posible alterar sus parámetros desde la ventana de comandos a menos que dichas variables sean
globalizadas por medio del comando global.
Así, en el programa cubica.m podemos conseguir que los parámetros a, b, c y d se introduzcan desde el
espacio de trabajo globalizando las variables a, b, c y d, de la forma:

function y=cubica(x)

% define una función cubica

global a b c d

y=a*x.^3+b*x.^2+c*x+d;

Una vez globalizadas, para dar valores concretos a los parámetros a, b, c y d habremos de escribir en la
ventana de comandos:

>> global a b c d

>> a=1 ; b=2 ; c=-1 ; d = 4 ;

Una vez definida la función cubica.m es posible, por ejemplo, calcular integrales numéricas.

Ejemplo.
3

∫ (2 x + 5 x 2 − 4 x + 6)dx .
3
Calcular una integral
1
Para ello será suficiente escribir en el espacio de trabajo las siguientes instrucciones:

>> global a b c d
>> a=2 ; b=5 ; c=-4 ; d = 6 ;
>> quad( 'cubica', 1 , 3 )

Ejercicio.
Obtener el mínimo de la función f( x) = 2 x 3 + 5 x 2 − 4 x + 6 mediante la instrucción

>> fmin ( 'cubica' , 0 )

Ejemplo

25
Construir una función de distribución N (0,1) , es decir, a construir la función
2
x 1 − s2
F ( x) = ∫ e ds
−∞
2π
2
1 − s2
Para ello construimos, previamente, la función f ( x ) = e , que denominamos gauss.m
2π

% Campana de Gauss

function y=gauss(x)

y=(1/sqrt(2*pi))*exp(-x.^2/2);

A continuación debemos construir otra función llamada normal.m , desde llamaremos a gauss.m para
realizar la integral.

% Funcion de distribucion normal

function y=normal(x)

y=quad('gauss', -10 , x );

Así,

>> normal ( 1 ) = 0.8413

Ejercicios.
Construir un archivo .m para definir cada una de las funciones que se representan a continuación.
Representar gráficamente cada una de estas funciones en el intervalo [-10,10].

a) f ( x) = 2 x 3 − 4 x 2 + 6 x − 3
b) g ( x) = x sen( x) + e x
x 3 + 5x + 1
c) h( x ) =
x2 +1
d) t ( x ) = x x + 1
2

x2
1 −2
e) s ( x ) = e
2π
1
f) c ( x ) =
π (1 + x 2 )
g) rect(x) = 1 si x ≤ 0.5 ; 0 en otro caso.

h) step(x) = 0 si x < 0; 1 en otro caso.

i) ramp(x) = o si x < 0; x en otro caso.

j) q(x)= 0 si x < 0; sen(πx / 2) si 0 ≤ x ≤ 1 ; 1 si x > 1.

26
 0.01t
k) Crear un archivo .m que dibuje conjuntamente las gráficas de las funciones 100 e y
100 (1 + 0.01) t que representan el crecimiento, con el tiempo, de un capital de 100 euros
colocado a un interés r = 0.01 , continuo y compuesto, respectivamente. Comparar el
crecimiento del capital en ambos tipos de capitalización.

Ejercicios.
Calcular el valor de las siguientes integrales definidas.

3
a) ∫−1
(2 x 3 − 4 x 2 + 6 x − 3)dx
3

∫ ( x sen( x) + e
x
b) )dx
1
1
x 3 + 5x + 1
c) ∫(
−1 x2 +1
)dx
2
d) ∫( x
0
x 2 + 1)dx

10 x2
1 −2
e) ∫
−10 2π
e dx

∞
dx
f) ∫ π (1 + x
−∞
2
)
g) Teniendo en cuenta que la función de densidad para una distribución de Cauchy viene dada
1
por la función f ( x ) = , calcular la probabilidad de que dicha variable aleatoria
π (1 + x 2 )
1
dx
este comprendida entre -1 y 1. ( ∫ π (1 + x
−1
2
)
).

Las funciones definidas hasta el momento tienen una única entrada y una única salida. Resulta muy
sencillo construir funciones con más de una entrada y más de una salida.

En primer lugar, veamos la sencilla función vectorial banana.m

f( x1 , x 2 ) = 100( x 2 − x12 ) 2 + (1 − x1 ) 2 .

% banana.m

function y=banana(x)

y=100*( x (2) – x (1). ^2). ^2 + ( 1 – x (1) ). ^2 ;

Como segundo ejemplo vamos a construir la siguiente transformación vectorial

y1 = 2 x12 − x2
y2 = x1 − x22

27
Para ello construimos el m-archivo de función, que llamaremos trans.m

function [y1,y2]=trans (x1,x2)

y1=2*x1.^2-x2; y2=x1-2*x2.^2;

Alternativamente esta función también puede representarse mediante el archivo trans1.m , de la forma

function y=trans1(x)

y=[2*x(1).^2-x(2); x(1)-2*x(2).^2];

Ejemplo.
Vamos a construir una función, que llamaremos stat.m , que tiene como entrada un vector y tres salidas
como son la longitud del vector, su media aritmética y su desviación típica.

function [ n media desv_stand ]=stat(x)

n=length(x);
media =sum ( x ) / n ;
desv_stand = sqrt ( sum ( ( x – media ) . ^ 2 ) / n ) ;

Ejercicio.
Escribir una función cuya entrada sea un vector y cuya salida sean los todos los momentos centrales,
hasta el orden 3.

Ejemplo.
Vamos a escribir una función cuya entrada sean los coeficientes de una ecuación de segundo grado y cuya
salida sean sus raíces. Pudiera ser lógico darle por nombre grado2f.m y su estructura como archivo .m es
la siguiente:

function [x1,x2]=grado2f(a,b,c)

%calcula las raíces de una ecuación de segundo grado

x1=(-b+sqrt(b^2-4*a*c))/(2*a);
x2=(-b-sqrt(b^2-4*a*c))/(2*a);

28
ALGORITMOS

Un algoritmo constituye una secuencia de operaciones aritméticas y lógicas que deben ser realizadas en
un determinado orden. En cada uno de los pasos de la secuencia de operaciones, se altera el valor de
alguna determinada variable empleando, para ello, el valor de dicha variable obtenido en la etapa anterior.
El proceso algorítmico continúa hasta que se verifiquen determinadas condiciones respecto a las variables
involucradas, aunque en ciertos casos no es posible saber de antemano el número de veces que se ha de
repetir el algoritmo. Los ordenadores sólo son capaces de llevar a cabo, de forma primaria, las reglas
elementales de la aritmética. Cualquier otro tipo de operación matemática más compleja ha de reducirse a
las reglas aritméticas elementales mediante el uso de los algoritmos.

Veamos, como ejemplo, un algoritmo destinado a obtener la raíz cuadrada de un número positivo “a”.
Partiendo de la ecuación x
2
= a , la transformamos en x 2 − 1 = a − 1 y por tanto en
a −1
( x + 1)( x − 1) = a − 1 , de donde se obtiene x − 1 = , es decir
x +1
a −1
x = 1+ (1) .
1+ x
La raíz cuadrada del número a=2 se obtendrá partiendo de un valor inicial de x, por ejemplo x=1,
colocándolo en el segundo miembro de la expresión (1) , obteniendo entonces un nuevo valor de x=3/2.
Este nuevo valor volverá a ser llevado al segundo miembro de la expresión 1 y así sucesivamente.

Este algoritmo puede, igualmente visualizarse , con el empleo de fracciones continuas como se muestra a
continuación.
a −1 a −1
x = 1+ = 1+
a −1 a −1
1+1+ 2+
1+ x 1+ x
y así sucesivamente se obtiene

a −1
x = 1+
a −1
2+
a −1
2+
a −1
2+
a −1
2+
2 + .....

La estructura algorímica que acabamos de describir puede ser tratada, de modo muy eficiente, mediante el
bucle for, que mostramos a continuación.

EL BUCLE FOR

Un bucle es una estructura computacional que permite repetir sucesivamente una serie de instrucciones.
El bucle for tiene la siguiente estructura:

for i =1 : n

INSTRUCCIONES

end

29
De modo que cuando un determinado índice i vaya desde el valor 1 hasta el valor n las instrucciones se
repetirán una vez por cada valor que avanza el índice.
Veamos algunos ejemplos elementales del uso del bucle for.

Ejemplo
En primer lugar construiremos un programa que permita calcular la raíz cuadrada de un determinado
número.

% Calcula la raíz cuadrada de un número a

x=1

for i=1:n
x=1+(a-1)/(1+x);
end

sprintf(‘la raíz cuadrada del número a es %g’,x)

Ejercicio
Establecer un algoritmo destinado a calcular la raíz n-ésima de un número.
Sugerencia: Partiendo de la ecuación x = a , usar para ello la identidad
n

x n − 1 = ( x − 1)( x n −1 + x n − 2 + ..... + x + 1) = a − 1

a −1
para concluir la expresión x = 1 +
1 + x + x 2 + ...... + x n

Ejemplo.
n
Escribir la tabla de los n primeros números naturales de la forma 2 .
Procederemos creando un archivo.m de la forma:

% Muestra la tabla de las n primeras potencias de 2

n=input( ‘número de potencias de 2 que deseamos encontrar = ‘ ) ;

xx = [ ];

for i = 1 : n

x = 2 ^ i;

xx= [ xx x ] ;

end

xx’

La lógica que hay en esta construcción es la siguiente: para cada valor del índice i el número 2 es elevado
a i y tal resultado se guarda dentro de la variable x. Finalmente los sucesivos valores de la variable x se
almacenan dentro de otra variable xx. Obsérvese que previamente hemos procedido a vaciar la variable
xx mediante la instrucción xx = [ ].

Ejemplo.

30
Escribir una función que sume las componentes de un vector dado. Tal función que podríamos llamar
vectsuma.m es equivalente a la función sum vista previamente, aunque esta última también es capaz de
actuar sobre matrices. Su estructura es la siguiente:

function y=vectsuma(x)

% suma las componentes de un vector x

n=length(x);

s=0;

for i= 1 : n

s=s+x(i);

end

y = s;

Ejercicios.
Escribir de archivos.m donde sea empleado el bucle for en su construcción y que realicen cada una de las
funciones que se detallan a continuación.

a) Suma de los n primeros números naturales.

b) [
Dibujar la curva seno en el intervalo 0,2π . ]
c) Escribir una función llamada medgeo.m que calcule la media geométrica de un vector x.
d) Determinación del producto de los elementos de la diagonal principal de una matriz de dimensiones
arbitrarias.
e) Una tabla de tres entradas en la que figuren los radios , las áreas y los perímetros de veinte círculos
con radios comprendidos entre 1 y 5.
π
f) Una tabla de los valores de la función sen(x) entre 0 y , haciendo uso del desarrollo de Taylor de
5
x3 x5
grado cinco ( sen( x) ≅ x − + ).
3! 5!
g) Escribir una función que sume los elementos pares de un vector.
h) Escribir una matriz de dimensiones arbitrarias cuyas i-ésima fila esté formada por los elementos i,
i+1, i+2, …… , i+n.
i) Escribir una matriz de dimensiones arbitrarias cuyas j-ésima columna esté formada por los elementos
j, j+1, j+2, …… , j+n.
j) Dado un vector x = ( x1 , x2 , ......, xn −1 , xn ) obtener el vector y = ( xn , xn −1 , ..... , x2 , x1 ) .

También es posible escribir diversos bucles for que estén anidados unos dentro de los otros tal como
muestra el siguiente ejemplo.

Ejemplo.
Escribir una matriz de dimensiones arbitrarias de modo que cada uno de sus elementos sea de la forma
1
aij =
i + j −1

31
% construye una matriz cuyos elementos son de la forma a(i,j)=1/(i+j-1)

n=input('número de filas =');

m=input('número de columnas =');

a = [ ];

for i=1:n

for j=1:m

a ( i , j ) = 1 / ( i + j - 1);

end

end

Ejercicios.
Haciendo uso del bucle for doblemente anidado, escribir archivos.m con las siguientes características:

a) Una función que sume todos los elementos de una matriz.

b) Una función cuya entrada sea una matriz y cuya salida sea un vector formado por la suma de cada
una de sus columnas (tal función coincide exactamente con la función sum).

EL BUCLE FOR EN LAS ECUACIONES EN DIFERENCIAS

Las ecuaciones en diferencias proporcionan ejemplos muy característicos e intuitivos del manejo del
bucle for.

Ejemplo.
Escribir un archivo.m que calcule los cincuenta primeros valores del problema de valores iniciales:
X t = 1.1 X t −1 con X 1 = 2 , dibujando finalmente el valor de la variable X frente al tiempo.
Llamaremos a este archivo malthus.m y tendrá la siguiente estructura:

% describe la dinámica del crecimiento malthusiano

xx=[ ]; x=2;

for i=1:50

x = 1.1 * x ;

xx =[ xx x ] ;

end

plot(xx)

32
La lógica que hay en esta construcción es la siguiente: para cada valor del índice i la variable escalar x se
multiplica por 1.1 y su resultado se almacena dentro de otra variable xx. Obsérvese que previamente
hemos procedido a vaciar la variable xx mediante la instrucción xx = [ ] .

MATLAB permite, igualmente, realizar este mismo programa sin el empleo de variables almacén como la
xx, debido a su capacidad de admitir variables vectoriales dentro de los bucles. El archivo malthus.m
puede ser escrito, análogamente, como

%describe la dinámica del crecimiento malthusiano

x=[ ]; x(1)=2;

for i=1:50

x(i+1) = 1.1 * x(i) ;

end

plot(xx)

Obsérvese que la variable x es una variable vectorial.

Las ecuaciones en diferencias de orden superior pueden ser resueltas de forma análoga.

Ejemplo.
Hallar la sucesión de los veinte primeros números de Fibonacci, que vienen dados mediante la ecuación
en diferencias X t + 2 = X t +1 + X t con las condiciones iniciales X 1 = 1 , X ( 2) = 1 .
Para ello crearemos el archivo fibo.m, que tiene la siguiente estructura:

% EL CODIGO DAVINCI

% Calcula los veinte primeros números de Fibonacci

%x=[ ];

x(1)=1 ; x(2)=1 ;

for i=1:20

x(i+2)=x(i+1)+x(i);

end

Ejercicio
xt +1 1 + 5
Comprobar que los términos de la sucesión de Fibonacci verifican lim = ≅ 1.618
n →∞ xt 2

Ejercicio.
Considérese la ecuación en diferencias X t + 2 − 1 / 6 X t +1 − 1 / 6 X t = 2 , con las condiciones iniciales
X (1) = 1 , X (2) = 1 / 2 .

33
a) Determinar sus 25 primeros términos mediante la creación de un archivo .m y un bucle for .
b) Escríbase, igualmente, la solución particular en forma analítica encontrando, previamente, la
solución general de la parte homogénea y la solución de equilibrio.
c) Dibujar la solución obtenida en el apartado a) y comparar con la solución particular obtenida en b).

El MATLAB resulta igualmente muy eficaz a la hora de tratar procesos estocásticos. Veamos a
continuación un ejemplo donde se maneja un paseo aleatorio con deriva.

Ejemplo.
Escribir un archivo .m que describa sucesivas realizaciones del paseo aleatorio con deriva
X t +1 = 0.7 X t + 10 + ε t .

% Paseo aleatorio con deriva

x=[ ] ; x ( 1 ) = 0 ;

for i=1:1000

x ( i + 1 ) = 0.7 * x ( i ) + 10 + randn;

end

plot(x);

Resulta igualmente inmediato el tratamiento de ecuaciones en diferencias no lineales como muestra el

siguiente ejemplo.

Ejemplo.
Construir un archivo .m que permita calcular las cincuenta primeras iteraciones de la ecuación en
π
diferencias logística X t +1 = 4 X t (1 − X t ) , con X (1) = .
4

% ecuación en diferencias logística

x =[ ] ; x ( 1 ) = pi / 4 ;

for i=1:50

x(i+1)=4*x(i)*(1–x(i));

end

plot(x)

Ejercicios.

a) Representar el espacio de fases de la ecuación en diferencias logística, dibujando X t frente a X t +1 .

b) Escribir un archivo.m donde se obtenga las iteraciones de la ecuación logística por medio de una
function que tenga tres entradas: el número de iteraciones, la condición inicial y el parámetro ( dada

34
 la logística X t +1 = mX t (1 − X t ) se denomina parámetro al coeficiente m, que en el ejemplo
anterior tomaba el valor m = 4).
c) Dibujar diversas series temporales, de longitud 250, obtenidas por la logística al dar al parámetro m
π
los sucesivos valores de m = 1 , m = 2 , m = 3 , m = 4 , con la condición inicial X (1) = .
4
d) Hallar el coeficiente de correlación entre las variables X t y X t −1 para diversas logísticas de
longitud 100, 500 y 1000.
e) Comparar las trayectorias de las series temporales de longitud 200 generadas por dos ecuaciones
π π
logísticas con parámetro m=4, con valores iniciales x 0 = y x0 = + 0.0001 ,
4 4
respectivamente.
f) Repetir el ejercicio anterior para series generadas por una logística con parámetro m=3.
π
g) Calcular las 100 primeras iteraciones de la ecuación logística con X (1) = para los diferentes
4
valores del parámetro m que se enumeran a continuación, comprobando que: para m = 3 el sistema
tiene un punto de equilibrio estable; para m = 3.1 hay un ciclo estable de periodo 2; para m = 3.5
hay un ciclo estable de periodo 4; para m = 3.55 hay un ciclo estable de periodo 8; para
m = 3.565 hay un ciclo estable de periodo 16; y para m = 3.57 la trayectoria es aperiódica.

Los sistemas de ecuaciones en diferencias, cuando están escritos en forma normal, pueden ser resueltos
fácilmente empleando un bucle for, tal como muestra el ejemplo que presentamos a continuación.

Ejemplo.
 xt +1   1 2  xt 
Hallar los veinte primeros puntos de la solución particular del sistema   =    , con
 y t +1   − 2 1  y t 
 x1  1
la condición inicial   =   .
 y1  1
Para ello vamos a construir el siguiente archivo que denominaremos discret.m .

% Genera un sistema discreto 2X2

x=[ ]; y=[ ];

x( 1) =1; y(1)=1;

for i=1:10

x(i+1)=x(i)+2*y(i);

y(i+1)=-2*x(i)+y(i);

end

[ x’ y’ ]

Este sistema puede expresarse alternativamente como

35
% discreto 2X2

a=[1 2;-2 1]; x=[ ];

x(: , 1)=[1 ; 1];

for i=1:10

x(: , i+1)=a*x( : , i);

end

Para dibujar la trayectoria en el espacio de fases teclearemos en el espacio de trabajo

>> plot ( x ( 1 , : ) , x ( 2 , : ) )

Para dibujar las trayectorias temporales de cada una de las variables haremos

>> plot ( x ( 1 , : ) )
>> plot( x ( 2 , : ) )

Ejercicio
Construyamos el sistema anterior con coeficientes aleatorios

% Sistema 2x2 con coeficientes aleatorios

a=randn(2,2);

x=[ ];

x(: , 1)=[1 ; 1];

for i=1:10

x(: , i+1)=a*x( : , i);

end

plot ( x ( 1 , : ) , x ( 2 , : ) )

eig (a)

ejecutar el archive sucesivas veces y adivinar a partir de la gráfica del espacio de fases la naturaleza real o
compleja de los autovalores así como si su parte real es mayor o menor que uno en valor absoluto.

Ejercicios.
1) Hallar los veinte primeros puntos de la solución particular del sistema
 xt +1   1 1 xt   3   x1   2 
  =    +   , con la condición inicial   =   .
y
 t +1   4 1  y t   − 2   y1   4 
2) Hallar las mil primeras iteraciones del sistema no lineal:

36
 xt +1 = 1 + yt − 1.4 xt2
yt +1 = 0.3 xt
 x1   0.1
con las condiciones iniciales   =   .
 y1   0.1
3) Hallar las mil primeras iteraciones del sistema no lineal:
xt +1 = 3.9( xt − xt2 )
y t +1 = 4 xt ( y t − y t2 )
 x1   0.1
con las condiciones iniciales   =   .
 y1   0.1
−1
4) Teniendo en cuanta que las matrices MAM y A tienen los mismos autovalores, poner ejemplos
de sistemas lineales de ecuaciones en diferencias 2x2 que tengan como espacio de fase
• Un nodo estable.
• Un nodo inestable.
• Un punto de silla.
• Una estrella estable.
• Una estrella inestable.

Modelo de población por grupos.

Consideremos el siguiente modelo poblacional humana en el que se consideran las siguientes variables
que representan el número de individuos de cada grupo de edad.

X 0 : número de recién nacidos,

X 1 : número de niños,
X 2 : número de jóvenes,
X 3 : número de adultos,
X 4 : número de jubilados.

Las ecuaciones que rigen la dinámica del modelo son

X 0 (t + 1) = a2 X 2 (t ) + a3 X 3 (t )
X 1 (t + 1) = b0 X 0 (t )
X 2 (t + 1) = b1 X 1 (t )
X 3 (t + 1) = b2 X 2 (t )
X 4 (t + 1) = b3 X 3 (t )

Donde los parámetros de fertilidad son a2 = 1.5 y a3 = 1 , y los parámetros de supervivencia son
b0 = 0.99 , b1 = 0.95 , b2 = 0.9 , b3 = 0.8 .
Consideremos ahora dos posibles estados iniciales de los que parte la población:
• Una población con grupos en equilibrio X 0 = X 1 = X 2 = X 3 = X 4 = 1 .
• Una población en la que sólo existen niños, es decir, X 1 = 1 , X 0 = X 2 = X 3 = X 4 = 0 .
En ambos casos las unidades se suponen expresadas en millones.

a) Estudiar la evolución de ambos grupos de población durante los próximos 20 años. Encontrar las
pirámides de población correspondientes. Comprobar que, con independencia de la proporción
inicial de individuos dentro de cada grupo, la distribución por edades de la población tiende al
autovector dominante. Dibujar la gráfica de la evolución de los diferentes grupos de población.

37
 b) Con el avance de la medicina y la falta de nacimientos, las pirámides de población de las
sociedades occidentales tienden a invertirse. Vamos a considerar ahora el caso en que parte de
los jubilados sobrevive pasando a una nueva categoría de jubilados ancianos. Para describir tal
fenómeno añadiremos una nueva ecuación X 41 (t + 1) = b3 X 3 (t ) + b4 X 4 (t ) , b4 = 0.7 ,
donde X 41 representa ahora el número total de jubilados, mientras que X 4 representa ahora el
número de personas recién jubilados. Estudiar ahora la población de jubilados y la nueva
pirámide poblacional.

A continuación vamos a considerar el siguiente modelo de población no lineal con reclutamiento, que
viene dado por el modelo

b
X t +1 = aX t + X t −d
1 + X tk− d
donde a, b, k y d son los parámetros del sistema. Para k = 10 se trata de una versión discreta del
modelo de Mackey-Glass, desarrollado originalmente para modelizar la producción y pérdida de glóbulos
blancos.
Este modelo puede ser interpretado como un modelo de dinámica poblacional. Si 0 < a < 1 y b > 0 y si
X t denota el número de adultos de la población, entonces a es la tasa de supervivencia de adultos y
d es el tiempo de retardo entre nacimiento y maduración de un individuo hasta alcanzar la edad
fértil.
b
El factor puede interpretarse como la tasa de fertilidad de los adultos nacidos hace d periodos
1 + X tk− d
en el pasado, la cual es no lineal como consecuencia de una fecundidad decreciente a altos niveles de la
b
población. Así, el término X t − d da cuenta del reclutamiento de nuevos adultos que han nacido
1 + X tk− d
hace d años.
Estudiar las 400 primeras evoluciones del modelo para los parámetros ( a, b, k , d ) = (0.2, 2, 6, 2) .

MODELO DE EPIDEMIAS
Como un ejemplo ilustrativo de la utilización de sistemas de ecuaciones en diferencias, estudiaremos un
modelo de epidemia bastante conocido: el sarampión.
Todos, de alguna u otra manera, hemos tenido contacto con esta enfermedad, altamente
contagiosa, que afecta especialmente a los niños. También sabemos que una vez superada, el organismo
se vuelve inmune a ella. Al niño no expuesto a la enfermedad (ni la ha sufrido ni ha sido debidamente
vacunado) se le denomina susceptible, y S será el número de niños que en un determinado momento y
lugar se encuentran en este estado.
Un problema que debemos considerar, es la duración de la enfermedad. Existe un periodo de
latencia después del contagio, durante el cual el virus no se ha desarrollado suficientemente, y el
individuo no es contagioso y no tiene síntomas. Pasado este periodo, el individuo se vuelve contagioso. I
representará al numero de individuos infectados que pueden contagiar a aquellos que entren en contacto
con ellos (I es el número de infecciosos). Al cabo de cierto tiempo los enfermos se vuelven inmunes y no
pueden reinfectarse.
Trataremos este problema haciendo algunas simplificaciones. Para ello consideremos que el
periodo de latencia es de una semana, y que permanecen en el estado contagioso, hasta que comienza su
recuperación y por tanto se vuelven inmunes, otra semana.

38
 Sea Ik el número de enfermos, con capacidad de contagiar, presentes en la semana k, mientras
que Sk será el numero de susceptibles esa semana. Ik+1 será el número de susceptibles que fueron
infectados por la enfermedad al principio de la semana k (recordemos que el periodo de latencia es de una
semana).
Sk+1 representa el numero de niños susceptibles de padecer la enfermedad en la semana k, menos
el numero de niños que han sido contagiados al comienzo de esta semana, mas el numero de nacimientos
B (niños que se incorporan al estado de susceptibles) durante la semana en cuestión. Este ultimo sumando
B, lo consideramos constante.
Finalmente, nos queda por señalar cuál es la capacidad de contagio de un niño infectado sobre
los demás. Para ello consideramos que un enfermo contagioso, en promedio, infecta de forma constante
una fracción f del total de susceptibles, y por tanto f Sk es la cantidad de susceptibles infectados por cada
uno de los Ik. El total de infectados, será por tanto, Ik f Sk.
Podemos seguir la evolución de un sistema de este tipo, a partir de unas ciertas condiciones
iniciales, mediante el siguiente sistema de ecuaciones en diferencias
Sk+1=Sk-fSkIk+B
Ik+1=fSkIk.
La situación de equilibrio, es decir, los puntos fijos, se obtienen resolviendo el sistema
Sk=Sk-fSkIk+B
Ik=fSkIk.
1
Dando como resultado los valores S*= ; I*=B.
f
I k +1 1
De la segunda ecuación, dividiendo por Ik, queda = fSk, y por tanto, si Sk < ; se verifica que
Ik f

I k +1 1
<1, luego I decrece. Si por el contrario, Sk> los valores de I crecerán. Lo ideal parece ser
Ik f
1
mantener la población de candidatos a padecer la enfermedad por debajo de , lo cual podría
f
conseguirse con una adecuada campaña de vacunación.
En el siguiente recuadro se propone un programa en MATLAB para simular este proceso.
S0=50000; f=0.00003;I0=15;B=200; % condiciones iniciales
S=[];I=[];S(1)=S0;I(1)=I0;
for i=1:499;
S(i+1)=S(i)-f*S(i)*I(i)+B;
I(i+1)=f*S(i)*I(i);
end
figure(1), plot([S' I'])
Ie = B; Se=1/f;
% espacio de fases
figure(2), plot(S,I,Se,Ie,’+’) % ciclo límite resultante alrededor del punto fijo

39
EL BUCLE WHILE
Mediante este bucle es posible repetir una determinada instrucción un número indeterminado de veces
mientras se verifique determinada condición. Su estructura sintáctica es la siguiente:

while CONDICIONES sobre ciertas variables

INSTRUCCIONES

end

Las condiciones que se emplean en el bucle while se describen por medio de los operadores relacionales
entre variables == , < , > , <= y >= .
Como primer ejemplo veamos como la ecuación en diferencias logística X t +1 = 4 X t (1 − X t ) , con
π
X (1) = , puede ser resuelta con un while.
4

%ecuación en diferencias logística

x = [ ] ; x ( 1 ) = pi / 4 ; i=1 ;

while i < 100

x(i+1)=4*x(i)*(1–x(i));

i=i+1;

end

Existen numerosos algoritmos computacionales donde el uso del while es imprescindible. Se trata de
aquellos procedimientos iterativos que se repetirán hasta haber alcanzado un determinado grado de
aproximación. Tal es el caso de los procedimientos iterativos de resolución de ecuaciones algebraicas o
de sistemas de ecuaciones lineales.

Ejemplo.
Encontrar una solución de la ecuación x + 3 x − 10 = 0 mediante un proceso de aproximaciones
5 2

sucesivas que mejore, sucesivamente, una primera aproximación dada de antemano.

El procedimiento iterativo al que hace referencia el enunciado es el siguiente:
Primeramente escribimos esta ecuación de la forma x = 10 − 3x 2 .
5

Seguidamente consideramos una primera aproximación x0 = 0.1 . Este valor se sustituye en el segundo
miembro de la ecuación para obtener una segunda aproximación x1 y así sucesivamente, es decir,

x0 = 0.1 ,
x1 = 5 10 − 3(0.1) 2 = 1.47577
x2 = 5 10 − 3(1.47577) 2 = 1.28225
x3 = 5 10 − 3(1.28225) 2 = 1.38344
x4 = 5 10 − 3(1.38344) 2 = 1.33613

40
El proceso continuará hasta obtener dos valores sucesivos de las aproximaciones suficientemente
cercanos, en tal caso se dice que el método converge, o hasta sobrepasar un número límite prefijado de
iteraciones en el caso que no haya convergencia.

El siguiente programa realiza el proceso de aproximaciones sucesivas que acabamos de describir. Para
abortar la ejecución del programa hemos empleado el relacional lógico & que es equivalente a un and.

X=[ ];

X(1)=2;x(2)=1;

i=1;

while abs ( x ( i + 1 ) – x ( i ) ) > 0.00000001 & i<1000

x ( i + 2 ) = ( 10 – 3 * x ( i + 1 ) ^ 2 ) ^ 0.2 ;

i=i+1;

end

n = length ( x )

x(n)

Al ejecutar este programa obtenemos un valor aproximado para la raíz de la ecuación de 1.3520 después
de realizar 28 iteraciones con una tolerancia de 0.00000001. No hemos permitido que el número de
iteraciones sea superior a 1000. El uso de la bifurcación if , que veremos a continuación, permitirá
introducir un mensaje advirtiendo cuándo el proceso iterativo es divergente.

LA BIFURCACIÓN CONDICIONAL IF

La bifurcación posibilita ejecutar una determinada instrucción en el caso en que se verifique una
determinada condición. Su forma general es la siguiente:

If CONDICIÓN sobre ciertas variables

INSTRUCCIÓN

end

Las condiciones que se emplean en la bifurcación if se describen por medio de los operadores relacionales
entre variables == , < , > , <= y >= .

Veamos un sencillo ejemplo

if x > 100

y = y + 1;

disp( 'x ha superado el umbral 100' );

end

41
En este ejemplo ha ocurrido lo siguiente:
Si la variable escalar x toma un valor superior a 100 entonces la variable y se incrementa en una unidad y
en la pantalla se muestra el mensaje “x ha superado el umbral 100” ( disp es una nueva instrucción que
permite mostrar mensajes en la pantalla).
Si la variable x no toma un valor superior a 100 se ignoran ambas instrucciones.

Con mucha frecuencia se desea que cuando no se verifican las condiciones señaladas en el if se ejecuten
ciertas instrucciones diferentes al caso en que si se verifican. En este caso la bifurcación condicional tiene
la siguiente sintaxis:

if CONDICION sobre ciertas variables

INSTRUCCIÓN 1

else

INSTRUCCIÓN 2

end

Ejemplo.
Vamos a construir un archivo.m , que llamaremos moneda.m , cuya ejecución simule el lanzamiento de
una moneda. Para ello haremos uso tanto de la bifurcación if como del generador de números aleatorios
gaussianos.

% simula el lanzamiento de una moneda

x=randn;

if x>=0

disp('cara');
disp(x);

else

disp('cruz');
disp(x);

end

En general la bifurcación if admite un conjunto de condiciones alternativas, cada una de ellas asociada
con una determinada instrucción, siguiendo la sintaxis que mostramos a continuación:

42
If CONDICIÓN sobre ciertas variables
INSTRUCCIÓN 1

elseif CONDICIÓN 2
INSTRUCCIÓN 2
................

................

elseif CONDICIÓN n
INSTRUCCIÓN n

else
INSTRUCCIÓN n+1

end

Como primer ejemplo vamos a considerar la función trozos.m que crea una función a trozos:
0.1x 2 + 3sen( x) − 3 x<0

f ( x ) =  2e x + 2 0≤ x≤3
 x2 x>3


function y=trozos(x)

for i=1:length(x)

if x(i)<0
y(i)=0.1*x(i).^2+3*sin(x(i))-3;

elseif 0<=x(i) & x(i)<=3

y(i)=2*exp(x(i))+2;

else %if 3<x(i)

y(i)=-x(i).^2;

end

end

Ejemplo.
Escribir un programa que simule el lanzamiento de dos monedas.
Por analogía con el caso de una moneda llamaremos a este programa moneda2.m . El programa es el
siguiente:

%simula lanzamiento de dos monedas

x=randn(2,1);

if x(1) >=0 & x(2) >=0

disp('c c');
disp(x);

elseif x(1) >=0 & x(2) <0

disp('c +');

43
disp(x);

elseif x(1) <0 & x(2) <0

disp('+ +');
disp(x);

else
disp('+ c ');
disp(x);

end

En el programa que mostramos a continuación, y que llamaremos monedas.m , simulamos el lanzamiento

de un número arbitrario de monedas contando el número de caras y cruces que resultan en cada
lanzamiento.

% simula el lanzamiento de n monedas contando el número de caras y cruces obtenidas

n=input( ' número de tiradas = ' ) ;

cara=0;

for i=1:n

x=randn;

if x>=0
cara=cara+1;
else
end

end

disp('nº de caras =');

disp(cara);
disp('nº de cruces =');
disp(n-cara);

Ejercicios.

1)Construir el anterior programa, de lanzamiento de monedas, como una function cuya entrada sea el
número de lanzamientos y cuya salida sea el número de caras obtenidas en los lanzamientos.
Construir entonces un nuevo programa que permita obtener el histograma que resulta del numero de caras
en el lanzamiento de veinte monedas, cuando dicha veinte monedas han sido lanzadas dos mil veces.

2) Simular el lanzamiento de un dado mediante un archivo.m .

3) Escribir un archivo.m capaz de efectuar la siguiente operación: Dada una matriz cuadrada A
determinar su rango y calcular su inversa. En caso en que al matriz no tenga inversa mostrar un
mensaje en la pantalla que lo indique.
→ →
4) Dado un sistema de ecuaciones A X = b , donde A no es necesariamente una matriz cuadrada
elaborar un archivo.m que discuta la naturaleza de las soluciones del sistema mediante el teorema de
Rouche_Frobenius.

44
5) Programar el método de la bisección de Bolzano para calcular raíces de una ecuación: si f(x) es una
función continua en un intervalo [ a, b ] y el signo de f(a) es distinto del signo de f(b) entonces existe
un número α ∈ ( a, b) de modo que f (α ) = 0.
6) Construir la tabla de los n primeros números primos.
7) Dado un número, averiguar si es primo.
8) Construir un archivo .m para la función
B( x ) = (1 − x 2 ) si x < 1 , mientras que B( x ) = 0 si x ≥ 1 . Dibujar su gráfica.
9) Idem para la función
3
W ( x ) = (1 − x ) 3 si x < 1 , mientras que W ( x ) = 0 si x ≥ 1 .

45
ECUACIONES DIFERENCIALES
Haciendo uso de las instrucciones ya examinadas es posible construir un archivo.m para resolver,
numéricamente, ecuaciones diferenciales por el método de aproximación de Euler. Tal como se muestra a
continuación, el siguiente programa que pudiera llamarse euler.m permite encontrar la solución
aproximada del problema de valores iniciales asociado a la ecuación logística
dy
= y (3 − y ) , y (0.1) = 0.2 .
dx

%MÉTODO DE EULER PARA RESOLVER ECUACIONES DIFERENCIALES

%y'=f(x,y);
%Cada nueva función f(x,y) debe introducirse en el FOR

%h=input('longitud del paso =');

%n=input('número de iteracciones =');
%x0=input('valor inicial para x =');
%y0=input('valor inicial para y =');
h=0.1;n=25;x0=0.1;y0=0.2;

x =[ ] ; y =[ ] ; x(1) = x0 ; y(1) = y0;

for i=1:n

%y(i+1)=y(i)+0.1*(x(i)^2-y(i)^2)*h;

y(i+1)=y(i)+(y(i)*(3-y(i)))*h;

x(i+1)=x(i)+h;

end

plot(x,y)

dy
Téngase en cuenta que para resolver numéricamente otra nueva ecuación como la = 0.1( x 2 − y 2 ) es
dx
preciso introducir la nueva forma funcional dentro del bucle for , haciendo activa la línea del programa

>> y(i+1) = y(i)+0.1*(x(i)^2 - y(i)^2)*h ;

El MATLAB emplea dos instrucciones, elaboradas específicamente, para resolver tanto ecuaciones
diferenciales como sistemas de ecuaciones diferenciales ordinarias.
La instrucción ode23 obtiene la solución aproximada, tanto de ecuaciones diferenciales como de sistemas,
empleando métodos de Runge-Kutta de segundo y tercer orden.
La instrucción ode45 obtiene la solución aproximada, tanto de ecuaciones diferenciales como de sistemas,
empleando métodos de Runge-Kutta de tercero y cuarto orden.
El funcionamiento de ambos procedimientos es similar.
dx
Supongamos que queremos resolver la ecuación diferencial logística = 2 x(5 − x) , con la
dt
condición inicial x(0) = 1.
El primer paso a dar es el de crear un archivo.m , que llamaremos logistic.m donde describiremos la
función que da lugar a la ecuación diferencial:

46
function dx = logistic(t,x)

%Esta función logistic.m evalúa la ecuación diferencial logística

dx = 2 * x .* ( 5 – x ) ;

Para resolver la ecuación diferencial es preciso teclear en el espacio de trabajo la siguiente instrucción:

>> [ t , x ] =ode23 ( ' logistic ' , ti , tf , X0 ) ,

o bien

>> [ t , x ] =ode45 ( ' logistic ' , ti , tf , X0 ) .

En ambos casos la instrucción requiere cuatro entradas que son:

El nombre del archivo.m que describe la ecuación escrito entre comillas simples.
El tiempo inicial ti donde comienza la integración.
El tiempo final tf donde finaliza la integración.
La condición inicial X0 que deseemos emplear.

La instrucción tiene además dos salidas que son los vectores columna t e x.
En t son almacenados los correspondientes instantes del tiempo donde se ha calculado la solución.
En x son almacenados los valores de la solución en dichos instantes del tiempo.
Por ejemplo podemos realizar la siguiente integración:

>> [ t , x ] =ode23 ( ' logistic ' , 0 , 10 , 1 )

o bien

>> [ t , x ] =ode45 ( ' logistic ' , 0 , 10 , 1 )

que corresponde a la elección ti=0, tf=10, X0=1.

La solución obtenida puede entonces ser examinada gráficamente si escribimos

>> plot(t,x)

o también

>> plot(x) .

Si queremos cambiar los valores de los parámetros que aparecen en la logistic.m esto puede hacerse sin
necesidad de acceder al archivo en el bloc de notas añadiendo la instrucción global , tal como se muestra
a continuación:

function dx = logistic(t,x)

%Esta función logistic.m evalúa la ecuación diferencial logística

global a b

dx = a * x. * ( b – x ) ;

En tal caso la introducción de los valores a=3 y b=4, por ejemplo, se realizará, simplemente, escribiendo
en la ventana de comandos :

47
>> global a b
>> a = 3 ; b = 4 ;

Ejercicios.
Resolver las siguientes ecuaciones diferenciales con las condiciones iniciales que se especifican.

dx
1) = 1.1x + 2 , con la condición inicial x(0) = 0.
dt
dx
2) = 1.1x 2 + sen(t ) , con la condición inicial x(0) = -1.
dt
dx
3) = 2t 2 + 4 , con la condición inicial x(0) = 1.
dt

La forma de resolver sistemas de ecuaciones diferenciales ordinarias es completamente análoga.

 .  −1 1   x 
x
Supongamos que deseamos resolver el sistema  .  =     junto con las condiciones
 y   − 1 − 1  y 
 
iniciales x(0)=1 , y(0)=3.
Primeramente crearemos una función para describir el sistema de ecuaciones diferenciales, que podríamos
llamar sis2x2.m .

function dx= sis2x2 ( t , x )

%Esta función evalúa un sistema lineal de ecuaciones diferenciales

dx = [ - x(1) + x(2) ; - x(1) – x(2) ] ;

A continuación formaremos un vector X0=[ -1 , 3 ] de condiciones iniciales, y eligiendo ti=0 y tf=50,

obtendríamos la solución de la forma:

>> [ t , x ] =ode23 ( ' sis2x2 ' , 0 , 50 , X0 )

En este caso la variable de salida t corresponde a los instantes de tiempo donde se ha calculado la
solución, mientras que y es una matriz de dos columnas, cada una de las cuales representa las funciones
x(t) e y(t) de la solución estimada por el programa.
En este caso es posible visualizar la solución a diversos niveles:
Dibujando la primera columna de la y con

>> plot( x ( : , 1 ) )

obtendremos la trayectoria de x(t).

Dibujando la segunda columna de la y con

>> plot( x ( : , 2 ) )

obtendremos la trayectoria de y(t).

El espacio de fases ( x(t) , y(t) ) puede ser obtenido mediante

>> plot( x ( : , 1 ) , x( : , 2 ) ) .

48
Podemos considerar los coeficientes de la matriz que define el sistema de ecuaciones diferenciales como
parámetros que pueden ser introducidos desde el exterior mediante el comando global :

function dx= sis2x2 ( t , x )

%Esta función evalúa un sistema lineal de ecuaciones diferenciales

global a11 a12 a21 a22

dx = [ a11* x(1) + a12* x(2) ; a21* x(1) + a22* x(2) ] ;

Ejercicio
Construyamos el sistema anterior con coeficientes aleatorios

function dx= sis2x2 ( t , x )

%Esta función evalúa un sistema lineal de ecuaciones diferenciales

a=randn(2,2)

dx = [ a(1,1)* x(1) + a(1,2)* x(2) ; a(2,1)* x(1) + a(2,2)* x(2) ] ;

ejecutar sucesivas veces las instrucciones

>> [ t , x ] =ode23 ( ' sis2x2 ' , 0 , 50 , X0 )

>> plot( x ( : , 1 ) , x( : , 2 ) )

y adivinar a partir de la gráfica del espacio de fases la naturaleza real o compleja de los autovalores así
como el signo de su parte real.

Ejercicios.
Resolver numéricamente y dibujar el espacio de fases de los siguientes sistemas de ecuaciones
diferenciales. Discutir, previamente, la configuración del espacio de fases de cada uno de estos sistemas
encontrando sus autovalores y autovectores mediante el uso del comando [ V,D]=eig(A) .

 .   0 4  x 
x
1)  .  =     , con la condiciones iniciales x(0)=-1 , y(0)=1 (Centro).
 y   − 9 0   y 
 
 .  − 2 −1  x 
x
2)  .  =     , con la condiciones iniciales x(0)=1 , y(0)=1 (Nodo estable).
 y   4 − 7   y 
 
 .   − 1 − 3  x 
x
3)  .  =     , con la condiciones iniciales x(0)=2 , y(0)=3 (Punto de silla).
 y   − 3 1   y 
 
 .  −1 1   x 
x
4)  .  =     , con la condiciones iniciales x(0)=-1 , y(0)=2 (Espiral estable).
 y   − 1 − 1  y 
 
.  0 2   x  1
x
5)  .  =     +   , con la condiciones iniciales x(0)=-1 , y(0)=1(Espiral estable).
 y   − 2 − 1  y   3 
 

49
  .  1 − 1  x   2 
x
6)  .  =     +   , con la condiciones iniciales x(0)=-1 , y(0)=0(Espiral estable).
 y   5 − 3   y   1 
 
.  2 1   x   1
x
7)  .  =     +   , con la condiciones iniciales x(0)=-1 , y(0)=1 (Centro).
 y   − 5 − 2   y  1
 
8) Comprobar que si en el sistema anterior sustituimos a22 = −2.1 el centro se convierte en una espiral
inestable. Si hacemos a22 = −1.9 el centro se convierte en una espiral estable.
−1
9) Teniendo en cuanta que las matrices MAM y A tienen los mismos autovalores, poner ejemplos
de sistemas lineales de ecuaciones diferenciales 2x2 que tengan como espacio de fase
• Un nodo estable.
• Un nodo inestable.
• Un punto de silla.
• Una estrella estable.
• Una estrella inestable.

El tratamiento de los sistemas 3x3 y, en general, nxn se realiza de forma completamente análoga a los
sistemas 2x2.
Si queremos tratar una ecuación diferencial de orden superior es preciso convertirla en un sistema, en
forma normal, mediante el consabido cambio de variables.
El MATLAB es especialmente eficiente en el tratamiento de sistemas y ecuaciones no lineales.

Ejemplo.
Consideremos la ecuación no lineal, de segundo orden, conocida como ecuación de Van der Pol:

.. .
x + ( x 2 − 1) x + x = 0

.
Este ecuación puede ser transformada en un sistema mediante el cambio de variables x = x1 , x = x 2 . El
sistema resultante será:
.
x1 = x 2
.
x 2 = x 2 (1 − x12 ) − x1

Para buscar sus soluciones aproximadas crearemos una función que llamaremos vdpol.m como la
siguiente:

function dx = vdpol(t,x)

dx=[ x(2) ; x(2).*(1 – x(1).^2) – x(1) ]

Para simular, finalmente, la solución de la ecuación diferencial de la ecuación de Van der Pol en el
intervalo 0 ≤ t ≤ 20 , con las condiciones iniciales x1 (0) = 1 / 3 , x 2 (0) = 0 , es suficiente escribir en
el espacio de trabajo las siguientes instrucciones:

>> ti=0 ; tf=20 ; x0=[1/3 0] ;

>> [ t , x ] = ode23( 'vpol ', t0 , tf , x0 );
>> plot(t,x)

Algunas ecuaciones no lineales tienen un comportamiento extraordinariamente complejo tal como puede
apreciarse resolviendo el siguiente sistema debido a Lorenz:

50
 .
x = 10( y − x)
.
y = 28 x − y − xz
. 8
z = − z + xy
3
Para la simulación de sus soluciones crearemos la función loren.m, que describimos a continuación:

function dx = loren ( t , x )

dx=[ 10*(x(2)-x(1)) ; 28*x(1)-x(2)-x(1)*x(3) ; (-8/3)*x(3)+x(1)*x(2) ];

Vamos a simular su solución aproximada en el intervalo de tiempo 0 ≤ t ≤ 20 , con las condiciones

iniciales x1 (0) = 0.2 , x 2 (0) = 0.4 , x3 (0) = 20 . Para ello es suficiente escribir en el espacio de
trabajo las siguientes instrucciones:

>> ti=0 ; tf=20 ; x0=[0.2 0.4 20] ;

>> [ t , x ] = ode23( 'loren' , t0 , tf , x0 );

El resultado que se almacena en la variable y es una matriz de tres columnas, cada una de las cuales
corresponde con la respectiva variable del sistema. La complejidad de la solución se contempla por medio
de los siguientes gráficos.

Trayectorias de las distintas variables frente al tiempo:

>> plot(t , x(:,1))

>> plot(t , x(:,2))
>> plot(t , x(:,3))

Representación bidimensional de cada una de las variables frente a otra:

>> plot(x(:,1) , x(:,2))

>> plot(x(:,1) , x(:,3))
>> plot(x(:,2) , x(:,3))

Espacio de fases tridimensional:

>> plot3(x(:,1) , x(:,2) , x(:,3))

51
5.2 ECUACIONES DIFERENCIALES

Muchos problemas de las ciencias aplicadas pueden formularse, matemáticamente, por medio
de la determinación de una función desconocida que satisface una relación determinada entre
ella y sus derivadas. Surgen, de este modo, ecuaciones cuyas incógnitas no representan algún
tipo de número sino una función. Tales ecuaciones se denominan ecuaciones diferenciales.

Consideremos, como ejemplo específico, la ecuación diferencial

dy
= 0.1( x 2 − y 2 )
dx
Para que esta ecuación tenga solución única es necesario imponer una condición inicial que
especifique el valor de la función incógnita y (x ) para un determinado valor x0 , es decir,
y ( x0 ) = y 0 . En este caso la ecuación diferencial se transforma en un problema de valores

iniciales como, por ejemplo, el siguiente:

dy
= 0.1( x 2 − y 2 ) , y (0) = 1 (5.1)
dx

EL MÉTODO DE EULER
MATLAB cuenta con potentes instrucciones propias para resolver numéricamente
ecuaciones diferenciales que emplearemos más adelante. No obstante, con el fin de
ejercitar las instrucciones ya examinadas, vamos a programar un sencillo algoritmo para
resolver, numéricamente, ecuaciones diferenciales conocido como método de
aproximación de Euler.
Consideremos un problema de valores iniciales de la forma
y ' = f ( x, y ) , y ( x 0 ) = y 0 (5.2)
Sea y = φ ( x ) la solución que deseamos buscar. Podemos considerar que una primera
aproximación a dicha solución en un punto x1 cercano a x0 es la recta tangente a y = φ ( x ) , cuya
pendiente tendrá el valor φ ' ( x0 ) = f ( x0 , y0 ) , por satisfacer la ecuación diferencial (5.2).
El valor aproximado y1 de φ ( x1 ) , puede obtenerse, según indica la Figura 5.5, de la forma
y1 = y0 + φ ' ( x0 )( x1 − x0 ) = y0 + f ( x0 , y 0 )( x1 − x0 ) .

Una vez calculado y1 es posible obtener el valor aproximado y 2 de la solución en otro

punto x 2 cercano a x1 , y así sucesivamente. En general la sucesión de puntos que se aproxima a

la solución puede obtenerse de forma recursiva mediante la expresión
y n +1 = y n + f ( xn , y n )( xn +1 − xn ) ,

52
 y2

según se muestra en la Figura 5.5.

y = φ (x )

y3

y1

y0

x0 x1 x2 x3

Si consideramos una separación uniforme entre los puntos xi en los que aproximamos la

solución, la expresión anterior toma la forma

y n +1 = y n + f ( x n , y n )h . (5.3)

El siguiente programa que mostramos a continuación y que hemos denominado euler.m

permite encontrar la solución aproximada del problema de valores iniciales (5.1) empleando el
algoritmo de Euler que se muestra en (5.3).
%METODO DE EULER PARA RESOLVER ECUACIONES DIFERENCIALES
%y'=f(x,y);
%Cada nueva función f(x,y) debe introducirse en el FOR
h=0.1;n=25;x0=0;y0=1;
x =[ ] ; y =[ ] ; x(1) = x0 ; y(1) = y0;
for i=1:n
y(i+1)=y(i)+0.1*(x(i)^2-y(i)^2)*h;
x(i+1)=x(i)+h;
end
plot(x,y)

La representación gráfica de la solución obtenida que se obtiene mediante la instrucción

53
plot(x,y) se muestra en la Figura 5.6.

Método de Euler
1.25

1.2

1.15

1.1

1.05

0.95

0.9
0 0.5 1 1.5 2 2.5 3

Figura 5.6

Téngase en cuenta que para resolver numéricamente otra nueva ecuación como, por ejemplo la
dy
= 3(2 x 2 − 4 y 2 ) , sería preciso introducir la nueva forma funcional dentro del bucle for del
dx
programa anterior.

5.2.2 LAS FUNCIONES ODE23 ODE45

MATLAB posee dos instrucciones específicas, para resolver tanto ecuaciones diferenciales
como sistemas de ecuaciones diferenciales ordinarias. En ambas instrucciones se hace uso de
algoritmos mucho más potentes que el método de Euler visto anteriormente.
Las funciones ode23 y ode45 obtienen la solución aproximada, tanto de ecuaciones
diferenciales como de sistemas, empleando métodos de Runge-Kutta de segundo y tercer orden
o cuarto y quinto orden respectivamente. El funcionamiento de ambas funciones es similar.

Ejemplo. Supongamos que queremos resolver el problema de valores iniciales (5.1). El primer
paso es crear un archivo.m, que llamaremos ecua.m donde describiremos la función que da
lugar a la ecuación diferencial:

function yprima=ecua(x,y)
%Esta función evalúa la ecuación diferencial
yprima=0.1*(x^2-y^2);

54
Para resolver la ecuación diferencial es preciso teclear en el espacio de trabajo la función:
» [x, y]=ode23('ecua',ti, tf, X0);
o bien
» [x, y]=de45('ecua', ti, tf, X0);

En ambos casos las funciones requieren cuatro entradas que son:

• El nombre del archivo.m que describe la ecuación escrito entre comillas simples.
• El valor inicial ti donde comienza la integración.
• El valor final tf donde finaliza la integración.
• La condición inicial X0 que deseemos emplear.
Además se dispone de dos salidas que son los vectores columna x e y.
• En x son almacenados los valores de la variable independiente donde se ha calculado la
solución.
• En y son almacenados los valores de la solución que se corresponde con los valores de
x.

Ejemplo. La integración que corresponde a la elección ti=0, tf=10, X0=1, es

» [ x , y ] =ode23 ( 'ecua' , 0 , 10, 1 )

La solución obtenida puede entonces ser examinada gráficamente si escribimos plot(x,y) o

también plot(y), tal como muestra la Figura 5.7.
Trayectoria de la ecuación diferencial (5.1)
10

1
0 2 4 6 8 10 12 14

Figura 5.7

Sea ahora una ecuación diferencial que depende de diversos parámetros

dy
= a (bx 2 − cy 2 ) , y (0) = 1 .
dx

55
Si queremos resolver esta ecuación para diferentes valores de los parámetros, MATLAB
permite los cambios de valores en los parámetros sin necesidad de acceder al archivo en el bloc
de notas. Para ello es suficiente añadir la instrucción global, tal como se muestra a
continuación.

function yprima=ecua(x,y)
%Esta función evalúa la ecuación diferencial
global a b c
yprima=a*(b*x^2–c*y^2);

En tal caso la introducción de los valores a=0.2, b=4 y c=3, se realizará, simplemente,
escribiendo en la pantalla del espacio de trabajo:
» global a b c
» a=0.2; b=4; c=3;

Ejercicio. Resolver las siguientes ecuaciones diferenciales con las condiciones iniciales que se
especifican.
dx
a) = 1.1x + 2 , con la condición inicial x(0)=0.
dt
dx
b) = 1.1x 2 + sen(t ) , con la condición inicial x(0)=-1.
dt
dx
c) = 2t 2 + 4 , con la condición inicial x(0)=1.
dt

A continuación consideraremos diversos ejemplos prácticos que ilustran el uso de las

ecuaciones diferenciales y su resolución por medio de MATLAB.

5.2.3 MODELOS DE CRECIMIENTO POBLACIONAL

Consideremos una población que tiene x(t ) individuos en el instante t . Llamemos r (t , x ) a tasa
de crecimiento de la población, es decir, la diferencia entre su tasa de nacimientos y mortalidad.
La tasa de cambio respecto al tiempo de dicha población obedecerá, entonces, a la siguiente
ecuación diferencial
dx
= r (t , x ) x .
dt
Existe multitud de formas de concretar la función r (t , x ) . La más sencilla consiste en
suponer que r (t , x ) = a , siendo a una constante. En este caso estaremos en presencia de la

56
ecuación
dx
= ax .
dt
Otra forma más sofisticada de especificar la tasa de crecimiento de la población es suponer
que
r (t , x ) = a (b − x ),

es decir, que dicha tasa de crecimiento es proporcional a la distancia que separa el tamaño de la
población de un determinado parámetro b . En tal caso la evolución de la población, a través del
tiempo, queda recogida por medio de la ecuación diferencial
dx
= a (b − x ) x .
dt
Como es lógico, el tamaño que adquiere una población en el transcurso del tiempo está muy
relacionado con el número de individuos del que inicialmente parte. Cualquiera de las
ecuaciones anteriores tiene, entonces, una infinidad de soluciones. Para garantizar la unicidad
de la solución es suficiente con imponer una determinada población inicial x0 en el instante t0 a
partir del que comience a contarse el tiempo. Las ecuaciones anteriores se transforman entonces
en problemas de valores iniciales de la forma:
dx
= ax , x (t0 ) = x0 , (5.4),
dt
o bien
dx
= ax (b − x ) , x (t0 ) = x0 . (5.5)
dt

Consideremos que los valores de los parámetros a y b son, respectivamente a = 0.05 y

b = 100 , y que la población parte de un tamaño inicial x0 = 3 . Haciendo uso de las instrucciones

de MATLAB para resolver ecuaciones diferenciales, es fácil obtener las gráficas de la Figura
5.8 donde se muestra la evolución de una población durante 10 unidades de tiempo, siguiendo
los modelos (5.4) y (5.5).

57
 Trayectoria ecuación diferencial (5.4)
5

4.8

4.6

4.4

4.2

3.8

3.6

3.4

3.2

3
1 2 3 4 5 6 7 8 9 10 11

Trayectoria ecuación diferencial (5.5)

120

100

80

60

40

20

0
0 5 10 15 20 25 30 35

5.2.5 Sistemas de ecuaciones diferenciales

La forma de resolver los sistemas de ecuaciones diferenciales es análoga a las ecuaciones
ordinarias.
Supongamos que deseamos resolver el sistema
 x   − 1 1   x 
  =     (5.6)
 y   − 1 − 1  y 
junto con las condiciones iniciales x(0)=1, y(0)=3.
Primeramente crearemos una función para describir el sistema de ecuaciones diferenciales,
que podríamos llamar sis2x2.m.

function y=sis2x2(t, x)
%Esta función evalúa un sistema lineal de ecuaciones diferenciales
y=[-x(1)+x(2); -x(1)–x(2)];

A continuación formaremos un vector X0=[-1, 3] de condiciones iniciales, y eligiendo ti=0 y

tf=50, obtendríamos la solución de la forma:

58
» X0=[-1, 3];
» [t, y ]=ode23('sis2x2', [0, 50], X0)

En este caso la variable de salida t corresponde a los instantes de tiempo donde se ha calculado
la solución, mientras que y es una matriz de dos columnas, cada una de las cuales representa las
funciones x(t) e y(t) de la solución estimada por el programa. En este caso es posible visualizar
la solución a diversos niveles. Dibujando la primera columna de la y obtendremos la trayectoria
de x(t) con la instrucción
» plot(y(:, 1))
Dibujando la segunda columna de la y obtendremos la trayectoria de y(t) con la instrucción
» plot(y(:, 2))
El espacio de fases (x(t), y(t)) puede ser obtenido mediante
» plot(x(:, 1), y(:, 1))
tal como se muestra en la Figura 5.12.

Sistema de ecuaciones diferenciales (5.6)

3

2.5

1.5

0.5

-0.5
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

Figura 5.12

Podemos considerar los coeficientes de la matriz que define el sistema de ecuaciones

diferenciales como parámetros que pueden ser introducidos desde el exterior mediante el
comando global

function y=sis2x2(t, x)
%Esta función evalúa un sistema lineal de ecuaciones diferenciales
global a11 a12 a21 a22
y=[a11*x(1)+a12*x(2);a21*x(1)+a22*x(2)];

Ejercicio. Resolver numéricamente y dibujar el espacio de fases de los siguientes sistemas de

59
ecuaciones diferenciales. Discutir, previamente, la configuración del espacio de fases de cada
uno de estos sistemas encontrando sus autovalores y autovectores mediante el uso del comando
[V,D]=eig(A).
 x   0 4  x 
a)   =     , con las condiciones iniciales x(0)=-1, y(0)=1.
 y   − 9 0   y 

 x  − 2 − 1  x 
b)   =     , con las condiciones iniciales x(0)=1, y(0)=1.
 y   4 − 7   y 

 x   − 1 − 3   x 
c)   =     , con las condiciones iniciales x(0)=2, y(0)=3.
 y   − 3 1  y 

 x  −1 1  x 
d)   =     , con las condiciones iniciales x(0)=-1, y(0)=2.
 y   − 1 − 1  y 

 x   0 2   x  1
e)   =     +   , con las condiciones iniciales x(0)=-1, y(0)=1.
 y   − 2 − 1  y   3 

 x   1 − 1  x   2 
f)   =     +   , con las condiciones iniciales x(0)=-1, y(0)=0.
 y   5 − 3   y   1 

 x   2 1  x  1
g)   =     +   , con las condiciones iniciales x(0)=-1, y(0)=1.
 y   − 5 − 2   y  1

Obsérvese que si queremos tratar una ecuación diferencial de orden superior es preciso
convertirla en un sistema, en forma normal, mediante un cambio de variables. MATLAB es
especialmente eficiente en el tratamiento de sistemas y ecuaciones no lineales.

Ejemplo. Consideremos la ecuación no lineal, de segundo orden, conocida como ecuación de

Van der Pol
2
x + ( x − 1) x + x = 0 .

Esta ecuación puede ser transformada en un sistema mediante el cambio de variables

x = x1 , x = x 2 . El sistema resultante será:

x1 = x 2
x 2 = x 2 (1 − x12 ) − x1

Buscaremos sus soluciones aproximadas creando una función que llamaremos vdpol.m como la
siguiente.

function y=vdpol(t,x)
y=[ x(2) ; x(2).*(1 – x(1).^2) – x(1) ]

Para simular la solución de la ecuación diferencial de la ecuación de Van der Pol en el intervalo

60
0 ≤ t ≤ 20 , con las condiciones iniciales x1 (0) = 1 / 3, x2 (0) = 0 , es suficiente escribir en el espacio

de trabajo, para obtener la representación gráfica de las trayectorias de las variables contenidas
en y que se muestra en la Figura 5.13, las siguientes instrucciones
» ti=0 ; tf=20 ; x0=[1/3 0] ;
» [t, y]=ode23('vdpol', [ti, tf], x0);
» plot(t,y)
Trayectorias de la ecuación diferencial de Van der Pol
3

-1

-2

-3
0 2 4 6 8 10 12 14 16 18 20

Figura 5.13

Como se observa en la Figura 5.14 las trayectorias del sistema son atraídas, con el transcurso
del tiempo, hacia una curva del espacio de fase de las variables almacenadas en y que se
denomina ciclo límite.
Espacio de fases tridimensional: ecuación de Van der Pol
3

-1

-2

-3
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Figura 5.14

El tratamiento de los sistemas 3x3 y, en general, nxn se realiza de forma completamente análoga
a los sistemas 2x2.

5.2.6 El atractor de Lorenz

Algunas ecuaciones no lineales tienen un comportamiento extraordinariamente complejo, tal

61
como puede apreciarse resolviendo el siguiente sistema debido a Lorenz
x = 10( y − x )
y = 28 x − y − xz
8
z = − z + xy
3
Para la simulación de sus soluciones crearemos la función lorenz.m, que describimos a
continuación.

function y=lorenz(t, x)
y=[ 10*(x(2)-x(1)) ; 28*x(1)-x(2)-x(1)*x(3) ; (-8/3)*x(3)+x(1)*x(2) ];

Vamos a simular su solución aproximada en el intervalo de tiempo 0 ≤ t ≤ 20 , con las

condiciones iniciales x1 (0) = 0.2 , x 2 (0) = 0.4 , x 3 (0) = 20 . Para ello es suficiente escribir en el
espacio de trabajo las instrucciones
» ti=0; tf=20; x0=[0.2 0.4 20];
» [t, y ]=ode23('lorenz', [t0, tf], x0);
El resultado que se almacena en la variable y es una matriz de tres columnas, cada una de las
cuales corresponde con la respectiva variable del sistema. La complejidad de la solución se
contempla por medio de los siguientes gráficos.
Las trayectorias de las distintas variables frente al tiempo se escriben como
» plot(t, y(:,1)), plot(t , y(:,2)), plot(t, y(:,3))

Trayectorias de las distintas variables frente al tiempo

20

-20
0 2 4 6 8 10 12 14 16 18 20
50

-50
0 2 4 6 8 10 12 14 16 18 20
60

40

20

0
0 2 4 6 8 10 12 14 16 18 20

Figura 5.15

Y la representación bidimensional de cada una de las variables frente a otra se expresa

62
» plot(y(:,1), y(:,2)), plot(y(:,1), y(:,3)), plot(y(:,2), y(:,3))

Trayectorias de las distintas variables enfrentadas

50

-50
-20 -15 -10 -5 0 5 10 15 20
60

40

20

0
-20 -15 -10 -5 0 5 10 15 20
60

40

20

0
-25 -20 -15 -10 -5 0 5 10 15 20 25

Figura 5.16

Espacio de fases tridimensional:

» plot3(y(:,1), y(:,2), y(:,3))

Espacio de fases tridimensional

50

40

30

20

10

0
40
20 20
0 10
0
-20
-10
-40 -20

Figura 5.17

Debido a la complejidad que presentan las trayectorias del espacio de fase del sistema de
Lorenz, a tal configuración se le denomina un atractor extraño. En este caso las trayectorias son
atraídas, no hacia un ciclo límite como ocurrían en el sistema de Van der Pol sino, hacia una
región del espacio de fase de carácter extraño e indefinido. Se dice, entonces, que la evolución
de las variables, respecto al tiempo, es de tipo caótico.

63