Grados de la EPSE

Prácticas de Matlab
Práctica 2

1. Resolución analítica de ecuaciones diferenciales de primer or-

den
Consideremos una ecuación de primer orden en forma normal:
dy
= f (t; y) (1)
dt
con la condición inicial
y (t0 ) = y0 :
Para resolver de forma exacta ecuaciones diferenciales utilizaremos el comando dsolve.
Si queremos hallar la solución general de una ecuación de primer orden la sintaxis es la
siguiente:
dsolve(eqn);
donde eqn es la ecuación que queremos resolver. Ésta se de…ne de la siguiente manera:

syms y(t)

eqn = dif f (y; t) == f (t; y)

Por ejemplo, para resolver la ecuación
dy
= 3y
dt
ponemos
syms y(t)
eqn = dif f (y; t) == 3 y
S = dsolve(eqn)
o directamente
S = dsolve(dif f (y; t) == 3 y)
y se obtiene
S = C1 exp(3 t):

1
 Si queremos resolver el problema de Cauchy pondremos

dsolve(eqn; condicion inicial):

Si resolvemos el problema de Cauchy

dy
= 3y; y (0) = 4;
dt
pondremos
syms y(t)
S = dsolve(dif f (y; t) == 3 y; y(0) == 4);
obteniendo
S = 4 exp(3 t):
El resultado S es una variable de caracteres que se puede asignar a una función:

f = inline(S):

A partir de la función f se pueden calcular los valores de la solución o dibujarla con el

comando plot:
Si la solución es un punto …jo, debemos hacer la grá…ca de otra forma, ya que en
ese caso la función inline no funciona correctamente. En el ejemplo anterior, tenemos el
punto …jo y1 (t) = 0. Para dibujar la grá…ca en el intervalo t 2 [0; 2] podemos hacer lo
siguiente:

t = 0 : 0;01 : 2;
y = 0: t:=t;
plot(t; y)

dy
Ejercicio 1 Dada la ecuación = (y 1)t :
dt
1. Hallar la solución general.

2. Hallar la única solución del problema de Cauchy con y (0) = 2 y asignársela a la

variable y:

3. Crear una función a partir de la solución y calcular y (0.75).

4. Dibujar la solución en verde en el intervalo [ 1; 1]:

5. Para la condición inicial y (0) = 1, dibujar la solución en negro en el intervalo

[ 1; 1]:

2
 Las ecuaciones de orden superior se resuelven de manera similar. Por ejemplo, la
sintaxis en una ecuación de segundo orden sería:

dsolve(eqn)

dsolve(eqn; [condicion inicial_1; condicion inicial_2])

La sintaxis de las derivadas es la siguiente:

du d2 u d3 u
= dif f (y; t); = dif f (y; t; 2); = dif f (y; t; 3); etc:
dt dt2 dt3
Por ejemplo, si queremos resolver el problema

y 00 + y = 0; y(0) = 1; y 0 (0) = 2;

ponemos
syms y(t)
z = dif f (y; t)
S = dsolve(dif f (y; t; 2) == y; [y(0) == 1; z(0) == 2]);
obteniendo
S = cos(t) + 2 sin(t):
Si la ecuación es de orden n la única diferencia es que se ponen n condiciones iniciales.

Ejercicio 2 Hallar la única solución del problema

y 00 2y 0 + y = et ; y (0) = 1; y 0 (0) = 0:

2. Resolución numérica de ecuaciones diferenciales

Consideremos una ecuación de primer orden en forma normal:
dy
= f (t; y) (2)
dt
con la condición inicial
y (t0 ) = y0 :
La función
[t; y] = ode45(g; [t0 ; t1 ]; inicial);
devuelve un conjunto de coordenadas t e y que representan la solución de la ecuación en
el intervalo t 2 [t0 ; t1 ]; y que se calculan usando métodos de Runge-Kutta de cuarto y
quinto orden.
Aquí
g = @(t; y) f uncion;
donde f uncion es la función f en (2), t es la variable independiente e y la variable
dependiente..

3
 Los valores de t0 y t1 especi…can los extremos del intervalo en el cual deseamos calcular
la solución.
El valor inicial especi…ca el valor inicial de la función en el extremo izquierdo del
intervalo, es decir, inicial=y (t0 ).
En este caso podemos utilizar indistintamente : , .^ y := o , ^ y =: para el producto,
división y elevación a una potencia.
Por ejemplo, si queremos resolver el problema y 0 = y t, y (1) = 2; en t 2 [1; 3];
pondremos
[t; y] = ode45(@(t; y) y t; [1; 3]; 2):
Finalmente, notemos que si escribimos

ode45(g; [t0 ; t1 ]; inicial);

entonces Matlab dibuja la grá…ca de la solución en lugar de dar datos numéricos.

Ejercicio 3 Resolvemos la la ecuación y 0 = 2ty en el intervalo [ 2; 2]: Teclear:

1. [t; y] = ode45(@(t; y) 2 y t; [ 2; 2]; 2)

2. [t; y] = ode45(@(t; y) 2 y t; [ 2; 2]; 2);

3. plot(t; y;0 r0 )

4. ode45(@(t; y) 2 y t; [ 2; 2]; 2)

¿Qué obtenemos en cada caso?

MATLAB tiene integrado su propio editor. MATLAB sólo puede ejecutar funciones
(archivos .m) que estén en sus librerías o en el directorio actual. Por ello es necesario
cambiar al directorio donde salvamos nuestro archivo antes de poder ejecutarlo.
Usando estos archivos podemos ejecutar una sucesión de comandos de golpe ejecutando
desde la terminal el archivo *.m en el que hemos escrito los mismos.

Ejercicio 4 Vamos a dibujar varias soluciones de la ecuación y 0 = 2ty en el intervalo

[ 2; 2] y en distintos colores.

1. Creamos un archivo nombre.m y lo guardamos.

2. Dibujamos la solución para la condición inicial 10:

hold on
axis([ 2; 2; 20; 20])
f = @(t; y) 2 y t
[t1; y1] = ode45(f; [ 2; 2]; 10);
plot(t1; y1;0 r0 )

4
 3. Dibujamos en el mismo grá…co las soluciones para las condiciones iniciales 15 y
20

4. Ejecutamos en la terminal de Matlab el archivo escribiendo el nombre del mismo.

En muchas ocasiones hemos de resolver una ecuación diferencial que depende de uno
o varios parámetros.
Por ejemplo, si queremos resolver el problema

y 0 = at2 y; y (0) = 2;

en el intervalo [0; T ] = [0; 10] con a = 3, pondremos

a=3
[t; y] = ode45(@(t; y) a t^2 y; [0; 10]; 2)

Ejercicio 5 Dada la ecuación y 0 = at2 y en el intervalo [0; 1.5] con la condición inicial
y (0) = 2:

1. Resolver la ecuación para a = 4:

2. Dibujar el grá…co de la solución.

Vamos a considerar la ecuación logística

y0 = y y2;
y (t0 ) = y0 :

Según el análisis cualitativo de las soluciones sabemos que, según sea la condición inicial,
las soluciones se comportan del siguiente modo:

1. y0 = 0 e y0 = 1 son puntos …jos.

2. Si y0 > 1, entonces y (t) ! 1 cuando t ! +1:

3. Si 0 < y0 < 1, entonces y (t) ! 1 cuando t ! +1:

4. Si y0 < 0, entonces y (t) ! 1 cuando t ! +1:

Los distintos tipos de soluciones se pueden ver en la Grá…ca 1.

Ejercicio 6 Obtener un grá…co con distintos tipos de soluciones de la ecuación logística

en el intervalo [ 2; 10]. Realizar los siguientes pasos:

1. Fijar los ejes [ 2; 10] [ 2; 2]:

2. Hallar las soluciones para y ( 2) = 1, y ( 2) = 0 y dibujarlas en color negro.

5
 2

1.5

0.5

-0.5

-1

-1.5

-2
-2 0 2 4 6 8 10

Figura 1: Soluciones de la ecuación logística

3. Hallar las soluciones para y ( 2) = 2 y dibujarla en color rojo.

4. Hallar las soluciones para y ( 2) = 0.01 y dibujarla en color azul.
5. Hallar las soluciones para y ( 2) = 0.01 y dibujarla en color verde.
6. Ponerle a la grá…ca el título "Soluciones de la ecuación logística".

Ejercicios para casa:

1. Dibujar en un grá…co distintas soluciones de la ecuación y 0 = y +e t en el intervalo
[0; 10]. Elegid un tamaño adecuado de los ejes. ¿A qué convergen las soluciones
cuando t ! +1?
2. Notemos que al usar el comando ode45 Matlab escoge el valor del paso h. Podemos
saber cuál este paso imprimiendo el archivo de tiempos t y viendo la diferencia
entre dos valores consecutivos. Sin embargo, hemos de tener en cuenta que Matlab
no siempre coge un tamaño de paso …jo, sino que este valor puede variar según el
momento del tiempo.

6
 a) Ejecutar [t; y] = ode45(@(t; y) 2 y t; [ 2; 2]; 10) y ver que el paso h no es …jo.
Para comprobarlo, podemos hacer lo siguiente:
b) Calculamos el número de …las de t: [m; n] = size(t), donde m es el número de
…las.
c) Calcular algunos valores de h tecleando, por ejemplo, t (2) t (1), t (3) t (2) ;
t(m) t(m 1), y cualquier t (j) t (j 1) con j m:

3. Ejecutar [t; y] = ode45(@(t; y) 1=y; [ 2; 2]; 1) y ver que el tamaño de h no es …jo.

4. Dada la ecuación y 0 = aty 2 en el intervalo [0; 0.9] con la condición inicial y (0) = 1,
hallar la solución para a = 2 y a = 2, y dibujar las dos soluciones en el mismo
grá…co con distintos colores.

3. Programar en Matlab (parte I)

Al igual que en los lenguajes de alto nivel, MATLAB permite crear programas utilizan-
do programación estructurada. Para ello cuenta con condicionales, bucles y funciones.

3.1. Funciones
Una función se de…ne mediante un …chero .m, cuyo nombre coincide con el de la
función. La primera línea ejecutable la palabra function. Su sintaxis es:

f unction argumentos_salida = nombre_f uncion(argumentos_entrada)

seguida de las instrucciones necesarias.

Cuando hay más de un argumento de salida, éstos deben ir entre corchetes y separados
por comas.
Es conveniente utilizar las primeras líneas del …chero para insertar un comentario
(iniciándolas con ’%’), explicando cómo debe usarse la función y sus argumentos (tanto
de entrada como de salida). Así, dicha de…nición será visible mediante la instrucción help
nombre-función.
La función puede …nalizarse en cualquier punto utilizando la instrucción return.
Una función utiliza las siguientes instrucciones para veri…car el número de argumentos:

nargin: número de argumentos de entrada que el usuario ha pasado a la función.

nargout: número de argumentos de salida que el usuario desea recibir de la función.

Las funciones se guardan en un archivo con extensión ".m". El nombredel …chero es

el que utilizamos para ejecutar el programa desde la línea de comandos de Matlab. Por
tanto, si tenemos el …chero Nombre.m, entonces escribiremos:

N ombre(Argumentos de entrada):

Ejemplo.

7
 Como ejemplo de lo anterior vamos a crear una función que calcule el valor de la
hipotenusa de un triángulo rectángulo a partir de sus dos catetos, así como la suma de
los dos catetos. A continuación se indican las instrucciones a programar. Las variable de
entrada serán los valores de los dos catetos, mientras que las variable de salida serán el
valor de la hipotenusa y la suma de los catetos.
En el archivo de nombre hipotenusa.m introducimos las siguientes líneas:

%La función hipotenusa calcula el valor de la hipotenusa de un triángulo dado el

valor de sus dos catetos, así como la suma de los dos catetos
function [hip,suma] = hipotenusa (cateto1, cateto2)
hip = sqrt(cateto1*cateto1+cateto2*cateto2);
suma=cateto1+cateto2;

Si ejecutamos el comando

[h; s] = hipotenusa(3; 4);

obtenemos
h = 5; s = 7:
La variable h contiene el valor de la hipotenusa y la variable s el valor de la suma de los
dos catetos x e y: Si ejecutamos

h = hipotenusa(3; 4)

sólo obtendremos el valor de la primera variable, es decir, el de la hipotenusa.

Es muy importante darle valor a las variables de entrada. Si escribimos [h; s] =
hipotenusa, Matlab imprimirá un mensaje de errorm ya que los valores de los catetos x; y
no están de…nidos.

Es importante resaltar que no es necesario (pero si conveniente) que el nombre de la

función coincida con el nombre del …chero. En cualquier caso, hay que recordar que para
ejecutar la función hemos de utilizar el nombre del …chero.

Ejercicio 7 Construir una función que calcule, dado el valor del radio r, el área del
círculo y la longitud de la circunferencia correspondientes.
La función tendrá una variable de entrada y dos de salida.
Escribir unos comentarios explicativos al principio del archivo.

3.2. Estructuras de control condicionadas

Permiten seleccionar entre dos conjuntos alternativos de instrucciones dependiendo de
que se veri…que una condición lógica (cuyo resultado es cierto o falso). Su sintaxis es de
la forma:

if condición
Instrucciones que deben ejecutarse si la condición es cierta.

8
 else
Instrucciones a ejecutar si no se veri…ca la condición anterior
end

Cuando no hay instrucciones que ejecutar si la condición no se cumple, la sintaxis

anterior se reduce a

if condición
Instrucciones que deben ejecutarse
end

Al contrario, cuando se encadenan varios bloques alternativos, la sintaxis queda como:

if condición_1
Instrucciones a ejecutar cuando se veri…ca la condición 1
elseif condición_2
Instrucciones a ejecutar cuando no se veri…ca la condición 1 y sí la condición_2
...
else
Instrucciones a ejecutar cuando no se veri…can las condiciones anteriores
end

Podemos imponer más de una condición, o condiciones complejas, utilizando los op-
eradores relacionales (condiciones cuyo resultado es cierto o falso) combinados con oper-
adores lógicos (sirven como nexo entre varios relacionales).

Símbolo
Símbolo relacional
lógico
Negación lógica
< Menor ~A
o complementario
<= Menor o igual A&B Operador de intersección \
> Mayor AjB Operador de unión [
Unión exclusiva
>= Mayor o igual xor(A; B)
(sólo uno ha de ser cierto)
Unión de varios elementos
x == y Igualdad any([A1 ; A2 ; :::; An ])
(Algún elemento ha de ser cierto)
Intersección de varios elementos
x~=y Desigualdad all([A1 ; A2 ; :::; An ])
(Todos han de ser ciertos)

Las operaciones <; ; >; aplicadas a números complejos comparan las partes reales
de los mismos.
Si aplicamos los símbolos lógicos a un vector la operación se aplica en cada elemento
del vector, obteniendo un vector de ceros y unos.

Ejemplo.

9
 Vamos a escribir un programa que devuelva el valor de la hipotenusa dados dos catetos
si el valor absoluto de los dos es inferior a uno, el valor del perímetro del triángulo si alguno
de los dos catetos es igual a uno, o cero si no se cumple ninguna de las dos condiciones.

%Devuelve el valor de la hipotenusa si los dos catetos son menores

%que uno en valor absoluto,
%el valor del perímetro del triángulo si alguno de los dos catetos
%es igual a uno,
%o cero si no se cumple ninguna de las dos condiciones
function H=HipCond(x,y)
if abs(x)<1 & abs(y)<1
H=sqrt(x^2+y^2);
elseif abs(x)==1 j abs(y)==1
Hip=sqrt(x^2+y^2);
H=Hip+x+y;
else
H=0;
end

Ejercicio 8 Escribir un programa que dibuje la solución del problema de Cauchy

(
dy
= y 2 1;
dt
y (t0 ) = y0;

para una condición inicial dada en el intervalo [t0 ; T ]. Seguir las siguientes indicaciones:

1. Las variables de entrada son la condición inicial y0 , el tiempo inicial t0 y el tiempo

…nal T:

2. Las variables de salida son los vectores t e y con los valores de la solución.

3. Hallar la solución del problema usando el comando ode45:

4. Si la condición inicial es 1 o 1, dibujar la solución en color negro.

5. Si la condición inicial es mayor que 1 o menor que 1, dibujar la solución en color

rojo.

6. Si la condición inicial está entre 1 y 1, dibujar la solución en color verde.

7. Usar el comando hold on, de forma que al usar el programa repetidas veces no se
borren las grá…cas.

8. Establecer los ejes [t0 ; T ] [ 4; 5]:

9. Ejecutar el programa en el intervalo [0; 5] con las condiciones iniciales 1; 1; 0.9; 0.5; 3; 1.1:

10
Ejercicios para casa:

1. Escribir un programa que devuelva el valor de la hipotenusa dados dos catetos si

el valor absoluto de alguno de los dos es inferior a uno. En caso de no cumplirse la
condición el valor devuelto será cero.

2. Escribir un programa en el que, dados una función f (x) y dos puntos a; b, se calcule
el punto medio del intervalo [a; b] si f (a) f (b) < 0. El programa devolverá el valor
cero en caso de que no se cumpla la condición. El programa tendrá tres argumentos
de entrada: la función f (x) y los valores a; b.
Nota: Para introducir la función f se escribirá entre comillas la función que de-
seemos estudiar. Así, si f (x) = x3 y nuestra función se llama PuntoMedio, entonces
pondremos
P untoM edio(0 x^30 ; a; b)
Por otro lado, dentro del programa utilizaremos la función inline para crear la
función. Quedaría de la siguiente manera:
function h=PuntoMedio(f,x,y)
g=inline(f);
(Resto de comandos)

11