Ecuaciones Diferenciales EDO MATLAB

UNIVERSIDAD DE CONCEPCION
FACULTAD DE CIENCIAS
FISICAS Y MATEMATICAS

DEPARTAMENTO DE INGENIERIA MATEMATICA

C´alculo Num´erico (521230)
Laboratorio 9
Ecuaciones Diferenciales Ordinarias
El objetivo de este laboratorio es aprender t´ecnicas para la resoluci´on num´erica de problemas de
valores iniciales (P.V.I.) para ecuaciones diferenciales ordinarias (E.D.O.) y sistemas de E.D.O.
Matlab tienevarios comandos para la resoluci´on num´erica de P.V.I. para E.D.O.:
Ordinary differential equation solvers.
(If unsure about stiffness, try ODE45 first, then ODE15S.)
ode45
- Solve non-stiff differential equations, medium order method.
ode23
- Solve non-stiff differential equations, low order method.
ode113
- Solve non-stiff differential equations, variable order method.
ode23t
- Solve moderatelystiff differential equations, trapezoidal rule.
ode15s
- Solve stiff differential equations, variable order method.
ode23s
- Solve stiff differential equations, low order method.
ode23tb - Solve stiff differential equations, low order method.
Como se ve en esta lista, hay m´etodos para resolver E.D.O. stiff y no stiff. Adem´as hay m´etodos de orden
bajo, medio, alto y variable.
Todos ellos tienenuna sintaxis semejante. Por ejemplo, para resolver el P.V.I.
y ′ = f (t, y),
y(to ) = yo ,
en el intervalo [to , tf ] mediante el comando ode45 en su opci´on m´as sencilla, debe ejecutarse:

[t,y]=ode45(’f’,[to tf],yo);
donde:
• f es el nombre de la funci´
on f (t, y) (t´ıpicamente definida mediante un programa function en un
archivo f.m);
• to y tf son los extremos del intervalo donde se deseaconocer la soluci´on;
• yo es el valor de la soluci´on en to (es decir el valor de la condici´on inicial y(to ) = yo );
• t devuelve los valores de la variable independiente t donde el m´etodo calcula el valor de la soluci´on;
• y devuelve los valores de la soluci´on en cada uno de los puntos t.
Estos comandos no requieren como dato un paso de integraci´
on h pues todos ellos determinan de
maneraautom´
atica en cada paso k, el tama˜
no del paso de integraci´
on hk necesario para mantener los
errores por debajo de una tolerancia determinada. Los valores de t que entrega corresponden a los puntos
tk = tk−1 + hk , k = 1, 2, . . ., en los que el comando necesit´o calcular el valor de y(tk ).
1

Si se desea conocer la soluci´on para ciertos valores de t, puede alternativamente ejecutarse:[t,y]=ode45(’f’,tspan,yo);
donde tspan es el vector de valores donde se desea conocer la soluci´on. Por ejemplo, tspan=0:0.1:1.
En ese caso, la salida t coincide con tspan e y contiene los valores de la soluci´on en esos puntos.
La tolerancia predeterminada de estos m´etodos es 10−3 , para el error relativo, y 10−6 , para el error
absoluto. Si se desea calcular la soluci´on con otras tolerancias,deben prefijarse las opciones elegidas
mediante el comando odeset. Adem´as, en la ejecuci´
on del comando para resolver la E.D.O., debe
agregarse el par´
ametro adicional de opciones. La sintaxis para realizar esto es, por ejemplo:

options=odeset(’RelTol’,1e-6,’AbsTol’,1.e-8);
[t,y]=ode45(’f’,[to tf],yo,options);
Si se ejecuta options=odeset(’RelTol’,1e-6,’AbsTol’,1.e-8) sin el “;” puede verse quehay otras
opciones que pueden prefijarse, adem´as de las tolerancias de los errores.
Por ejemplo, si se desea resolver el P.V.I.
y ′ = y,
y(0) = 1,
en el intervalo [0, 1.5], mediante el comando ode45 y visualizar la soluci´on obtenida, debe crearse un
fichero f.m como sigue:

function z=f(t,y)
z=y;
y ejecutarse:

[t,y]=ode45(’f’,[0 1.5],1);
plot(t,y)
As´ı se obtiene la siguiente gr´
afica:
4.5

43.5

3

2.5

2

1.5

1

0

0.5

1

1.5

El siguiente ejemplo resuelve la misma ecuaci´
on en los puntos t=0:0.1:1.5, con error absoluto menor
a 10−6 y calcula los errores cometidos restando los valores calulados a los de la soluci´on verdadera, que
en este caso es y(t) = et :
2

options=odeset(’AbsTol’,1.e-6);
tspan=0:.1:1.5;
[t,y]=ode45(’f’,tspan,1,options);
error=exp(t)-y
error =
1.0e-06 *
0...