Solución análítica de ecuaciones diferenciales ordinarias

En esta entrada usaremos MatLab para encontrar la solución analítica de una ecuación diferencial ordinaria. Para ello utilizaremos la función dsolve.
--------------------------------------------------------------------------------------------------------------------------------
dsolve
Soluciona ecuaciones diferenciales ordinarias y sistemas de ecuaciones diferenciales ordinarias
Sintaxis
S = dsolve(eqn)
S = dsolve(eqn,cond)
S = dsolve(eqn,cond,Name,Value)
Y = dsolve(eqns)
Y = dsolve(eqns,conds)
Y = dsolve(eqns,conds,Name,Value)
[y1,...,yN] = dsolve(eqns)
[y1,...,yN] = dsolve(eqns,conds)
[y1,...,yN] = dsolve(eqns,conds,Name,Value)
Descripción 
En esta entrada solo describiremos las dos primeras.
1). S = dsolve(eqn): Soluciona la ecuación diferencial ordiaria eqn. Aquí eqn es una ecuación simbólica que contiene diff para indicar derivadas.
 2). S = dsolve(eqn,cond): Soluciona la ecuación diferencial ordinaria eqn con la condición inicial o condición de contorno cond.

Las variables simbólicas se pueden definir con la función syms. Alternativamente, se puede usar un string con la letra D para indicar derivadas.
----------------------------------------------------------------------------------------------------------------------------- 
Ejemplo 1: Obtengamos la solución particular de la ecuación diferencial ordinaria de primer orden:

                            y'= sin(2*x)-y*tan(x);   y' = dy/dx,
con la condición inicial
                            y(0)=1.
----------------INICIO------------------------------------------------------------------------------------------------------

%---creamos las variables simbólicas con syms:
syms y(x)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve:
y(x) = dsolve(diff(y,x) == sin(2*x)-y*tan(x),y(0)==1);
%----Evaluamos y para valores de x entre 0 y 4 a pasos de 0.1 y los valores lo asignamos a la variable y1: 
x = 0:0.1:4;
y1=eval(y);
%----graficamos y1 usando la función plot:
figure
plot(x,y1)
grid on
%---Evaluamos y para el valor x = 1 usando eval: 
 x=1;
yo=eval(y);
%---graficamos el punto  (x,yo) sobre la gráfica de y(x), de color rojo 'r' y marca 'o' :
hold on
plot(x,yo,'ro')
xlabel('x')
ylabel('y')
grid on

-------------------FIN-----------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
y(x) = 3*cos(x)-cos(2x)-1
Evaluada para x=0 resulta:
yo = 1.0371

La gráfica resultante es:


------------------------------------------------------------------------------------------------------------------------------
Ejemplo 2: Obtengamos la solución particular de la ecuación diferencial ordinaria de segundo orden:
                          h'' + 2h'+2h = 0;  h' = dh/dt,
con la condición inicial
                           h(0)=0  ;  h''(0)=15 .
----------INICIO-------------------------------------------------------------------------------------------------------
%---creamos las variables simbólicas con syms:
syms h(t)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve. Usamos D en lugar de diff para definir la ODE colocándola como string (entre comillas simples):
h(t) = dsolve('D2h == -2*Dh-2*h', h(0)==0, 'Dh(0)==15')
%----Evaluamos h para 100 valores de t entre 0 y 10 y lo asignamos a la variable h1:
N=100;
t = linspace(0,10,N);
h1=eval(h);
%----graficamos h1 usando la función plot:
figure
plot(t,h1)
grid on
%---Evaluamos h para el valor t = 2:
t=2;
ho=eval(h);
%---graficamos el punto (0,ho) sobre la gráfica de h(t), de color rojo 'r' y marca '*' :
hold on
plot(t,ho,'r*')
xlabel('t')
ylabel('h')
grid on
-------------FIN-------------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
h(t) = 15*exp(-t)*sin(t)
Evaluada para t=2 resulta:
ho = 1.8459
La grafica resultante es:

Gráficas de campos direccionales con MatLab

Esta rutina en Matlab realiza la gráfica de un campo direccional a partir de la ecuación diferencial :

y'=f(x,y)=-5*x.^4.*y.^2.

Sobre el campo de direcciones se realiza la gráfica de la solución particular:

y = 1./(x.^5+1).

%=====INICIO========
%-1)---Región--------
x1=0;
x2=2;
y1=0;
y2=2;
N = 20;
x = linspace(x1,x2,N);
y = linspace(y1,y2,N);
[X,Y]= meshgrid(x,y);

%-2)---calculo del angulo de inclinación de cada vector
m = -5*X.^4.*Y.^2;
T = atan(m);

%-3)--componentes de cada vector
U = cos(T);
V = sin(T);

%-4)--Dibujo de los vectores tangentes
h = quiver(X,Y,U,V);

%-5)--superposición en el campo de direcciones de la solución particular
hold on
yp = 1./(x.^5+1);
plot(x,yp)
grid on
%===FIN========

El resultado es el siguiente:

Fig.1. Campo de direcciones (azul). Solución particular (rojo)