% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root = ../Metodos_Numericos.tex  

\chapter[Ecuaciones diferenciales ordinarias]{Ecuaciones diferenciales ordinarias}

En este capítulo se lleva a cabo la solución de ecuaciones diferenciales de la forma

$$
\frac{dy}{dx} = f(x,y)
$$

En general los métodos que veremos hacen una formulación lineal de la forma

$$
y^{(k+1)} = y^{(k)} + \phi h
$$

donde $\phi$ es una estimación de la pendiente de la línea recta y utilizaremos esta linea para hacer una predicción de los valores de la integral.

\section{Integración por el método de Euler}

La primera derivada proporciona una estimación directa de la pendiente en $x^{(k)}$ dada por

$$
\phi = f(x^{(k)}, y^{(k)})
$$

donde $f(x^{(k)}, y^{(k)})$ es la ecuación diferencial evaluada en $x^{(k)}$ y $y^{(k)}$. Utilizando esta estimación tenemos

$$ 
y^{(k+1)} = y^{(k)} + f(x^{(k)}, y^{(k)}) h
$$

La ecuación anterior es conocida como el método de Euler (o de Euler-Cauchy o de punto medio). Su implementación en Matlab es:

\begin{verbatim}
function [x,y] = Euler(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N
        y(k+1,:) = y(k,:) + f(x(k), y(k,:))*h;
        x(k+1) = x(k) + h;
    end;
end
\end{verbatim}

donde x\_ini el el vector de valores iniciales, h es el incremento, N es el conjunto de valores donde se integra y f son el conjunto de ecuaciones diferenciales.

\section{Integración por el método de Heun con solo uno y con varios predictores}

Un método para mejorar la estimación de la pendiente involucra la determinación de dos derivadas para cada uno de los  intervalos (una en el punto inicial y otra en el punto final). Las dos derivadas se promedian con el fin de obtener una estimación mejorada de las pendientes para todo el intervalo. 

Recordemos que el método de Euler, la pendiente al inicio de un intervalo es

$$
y'^{(k)} = f(x^{(k)}, y^{(k)})
$$
y se usa para extraporlar linealmente a $y^{(k+1)}$

$$
y_p^{(k+1)} = y^{(k)} + f(x^{(k)}, y^{(k)}) h
$$

La pendiente en el punto final del intervalo la podemos calcular como:
$$
y'^{(k+1)} = f(x^{(k+1)}, y_p^{(k+1)})
$$

y la pendiente promedio la podemos calcular como

$$
\bar{y}' = \frac{y'^{(k+1)}+ y'^{(k)}}{2} 
$$

De lo anterior podemos decir que el método tiene dos pasos, uno de predicción que es una iteración de Euler y un Corrector donde se utiliza la estimación de la pendiente al final del bloque. Estos pasos son:

Predicción
$$
y_p^{(k+1)} = y^{(k)} + f(x^{(k)}, y^{(k)}) h
$$

Corrección
$$
y^{(k+1)} = y^{(k)} + \frac{f(x^{(k)}, y^{(k)}) + f(x^{(k+1)}, y_p^{(k+1)})}{2} h
$$

La implementación en Matlab para el Método de Heun es:

\begin{verbatim}
function [x, y] = Heun(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N
        %Predictor
        k1 = f(x(k),   y(k,:)       );
        k2 = f(x(k)+h, y(k,:) + k1*h);
        %corrector
        %disp([k1, k2, y + k1]);
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1 + k2 )*h/2.0;        
    end;
end\end{verbatim}

\section{Integración por el método del punto medio}

Esta técnica usa el método de Euler para predecir un valor de $y$ en el punto medio del intervalo

$$
y^{(k + 1/2)} = y^{(k)} + f(x^{(k)}, y^{(k)}) h/2
$$

Después este valor es utilizado para calcular la pendiente en el punto medio y hacer la solución en todo el intervalo como

$$
y^{(k+1)} = y^{(k)} + f(x^{(k+1/2)}, y^{(k+1/2)}) h
$$

La implementación en Matlab es

\begin{verbatim}
function [x,y] = Punto_Medio(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N

        xm = x(k) + h/2;
        ym = y(k,:) + f(x(k), y(k,:))*h/2;

        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + f(xm, ym)*h;        
    end;
end
\end{verbatim}

\section[Runge-Kutta 2do orden]{Fórmulas de integración por algunos métodos de Runge-Kutta de segundo orden }

Los algoritmos de Runge-Kutta son los métodos de integración con mejor desempeño y exactitud. En el caso del método de Runge-Kutta tenemos calcular dos valores $k_1$ y $k_2$ utilizando

$$
k_1 = f(x^{(k)}, y^{(k)})
$$

$$
k_2 = f(x^{(k)} + \frac{3}{4} h, y^{(k)}+\frac{3}{4}k_1 h)
$$

La actualización se lleva a cabo haciendo

$$
y^{(k+1)} = y^{(k)} + \left ( \frac{1}{3} k_1 + \frac{2}{3} k_2 \right) h
$$


\begin{verbatim}
function [x,y] = RK(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N     
        k1 = f(x(k), y(k,:));
        k2 = f(x(k) + 3*h/4, y(k,:)+3*h*k1/4);
        
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1/3 + 2*k2/3)*h;        

    end;
end
\end{verbatim}

\section[Runge-Kutta 3er orden]{Fórmula de integración por el método de Runge-Kutta clásico de tercer orden}

En el caso del método de Runge-Kutta de tercer orden tenemos que calcular tres valores $k_1$, $k_2$ y $k_3$ utilizando las siguientes formulas

$$
k_1 = f(x^{(k)}, y^{(k)})
$$

$$
k_2 = f(x^{(k)} + \frac{1}{2} h, y^{(k)}+\frac{1}{2}k_1 h)
$$

$$
k_3 = f(x^{(k)} + h, y^{(k)} - k_1 h + 2 k_2 h)
$$

La actualización se lleva a cabo haciendo

$$
y^{(k+1)} = y^{(k)} + \frac{1}{6} \left ( k_1 + 4 k_2  + k_3 \right) h
$$

\begin{verbatim}
function [x, y] = RK3(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N
       
        k1 = f(x(k),       y(k,:));
        k2 = f(x(k) + h/2, y(k,:) + k1*h/2);
        k3 = f(x(k) + h,   y(k,:) - k1*h + 2*k2*h);
 
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1 + 4*k2 + k3)*h/6;        
    end;
end
\end{verbatim}

\section[Runge-Kutta 4to orden]{Formula de integración por el método de Runge-Kutta clásico de cuarto orden}

El más popular de los métodos RK es el de cuarto orden. En el caso del método de Runge-Kutta de cuarto orden tenemos que calcular cuatro valores $k_1$, $k_2$, $k_3$ y $k_4$ utilizando las siguientes formulas

$$
k_1 = f(x^{(k)}, y^{(k)})
$$

$$
k_2 = f(x^{(k)} + \frac{1}{2} h, y^{(k)}+\frac{1}{2}k_1 h)
$$

$$
k_3 = f(x^{(k)} + \frac{1}{2} h, y^{(k)}+\frac{1}{2}k_2 h)
$$

$$
k_4 = f(x^{(k)} + h, y^{(k)}+ k_3 h)
$$
La actualización se lleva a cabo haciendo

$$
y^{(k+1)} = y^{(k)} + \frac{1}{6} \left ( k_1 + 2 k_2  + 2 k_3 + k_4\right) h
$$

La implementación en Matlab para el método de Runge-Kutta de cuarto orden es:

\begin{verbatim}
function [x, y] = RK4(x0, y0, h, N, f)
    D = length(y0);
    
    x = zeros(N+1, 1);
    y = zeros(N+1, D);
    
    x(1) = x0;
    y(1,:) = y0;
    
    for k = 1:N
       
        k1 = f(x(k),       y(k,:));
        k2 = f(x(k) + h/2, y(k,:) + k1*h/2);
        k3 = f(x(k) + h/2, y(k,:) + k2*h/2);
        k4 = f(x(k) + h,   y(k,:) + k3*h);
        
        x(k+1) = x(k) + h;
        y(k+1,:) = y(k,:) + (k1 + 2*k2 + 2*k3 + k4)*h/6;        
    end;
end
\end{verbatim}

\section{Ejemplo}

\subsection {Una ecuación diferencial sencilla}

Consideremos la ecuación diferencial

$$
f(x,y) = -2x^3 + 12x^2 -20x +8.5
$$

La integral la podemos calcular como

$$
\int \frac{dy}{dx} dx = \int f(x,y) dx = -0.5x^4 + 4x^3 -10x^2 +8.5 x 
$$

La solución para $y$ para un valor inicial $y(0) = 1$ lo cual da la siguiente expresión

$$
y(x) = -0.5x^4 + 4x^3 -10x^2 +8.5 x +1
$$

Encontrar la solución numérica para esta ecuación diferencial, suponiendo un valor inicial $x^{(0)} = 0$, $y^{(0)} = 1$ y $h= 0.5$, utilizando los métodos analizados en este capítulo.

\textbf{Solución}

La implementación de la función en Matlab es

\begin{verbatim}
    function dy = ecuacion(x, y)
        dy = -2*x^3 + 12*x^2- 20*x + 8.5;
    end
\end{verbatim}

En la siguiente tabla se muestran los resultados al aplicar la integración numérica.

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
    \hline
           &               & Real & Euler & Heun & P. Medio & RK  & RK3 & RK4\\
    $k$ & $x^{(k)}$ & $y^{(k)}$ & $y_{E}^{(k)}$ & $y_{H}^{(k)}$ &$y_{PM}^{(k)}$& $y_{RK}^{(k)}$& $y_{RK3}^{(k)}$& $y_{RK4}^{(k)}$\\
    \hline
    0 & 0.0  &  1.0000  &  1.0000  &  1.0000  &  1.0000  &  1.0000  &  1.0000  &  1.0000\\
    1 & 0.5  &  3.2188  &  5.2500  &  3.4375  &  3.1094  &  3.2773  &  3.2188  &  3.2188\\
    2 & 1.0  &  3.0000  &  5.8750  &  3.3750  &  2.8125  &  3.1016  &  3.0000  &  3.0000\\
    3 & 1.5  &  2.2188  &  5.1250  &  2.6875  &  1.9844  &  2.3477  &  2.2188  &  2.2188\\
    4 & 2.0  &  2.0000  &  4.5000  &  2.5000  &  1.7500  &  2.1406  &  2.0000  &  2.0000\\
    5 & 2.5  &  2.7188  &  4.7500  &  3.1875  &  2.4844  &  2.8555  &  2.7188  &  2.7188\\
    6 & 3.0  &  4.0000  &  5.8750  &  4.3750  &  3.8125  &  4.1172  &  4.0000  &  4.0000\\
    7 & 3.5  &  4.7188  &  7.1250  &  4.9375  &  4.6094  &  4.8008  &  4.7188  &  4.7188\\
    8 & 4.0  &  3.0000  &  7.0000  &  3.0000  &  3.0000  &  3.0312  &  3.0000  &  3.0000\\
    \hline
\end{tabular}
\end{center}

Para evaluar la exactitud de los métodos calculamos el valor absoluto de la diferencia entre el valor real y el calculado.

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
    \hline
    $k$ & $|y^{(k)}-y_{E}^{(k)}|$ & $|y^{(k)}-y_{H}^{(k)}|$ &$|y^{(k)}-y_{PM}^{(k)}|$& $|y^{(k)}-y_{RK}^{(k)}|$& $|y^{(k)}-y_{RK3}^{(k)}|$& $|y^{(k)}-y_{RK4}^{(k)}|$\\
        \hline
    0 &       0     &    0        &  0    &     0     &    0      &   0\\
    1 &  2.0312  &  0.2188  &  0.1094  &  0.0586  &       0    &     0\\
    2 &  2.8750  &  0.3750  &  0.1875  &  0.1016  &       0    &     0\\
    3 &  2.9062  &  0.4688  &  0.2344  &  0.1289  &       0    &     0\\
    4 &  2.5000  &  0.5000  &  0.2500  &  0.1406   &      0    &     0\\
    5 &  2.0312  &  0.4688  &  0.2344  &  0.1367   &      0    &     0\\
    6 &  1.8750  &  0.3750  &  0.1875  &  0.1172    &     0    &     0\\
    7 &  2.4062  &  0.2188  &  0.1094  &  0.0820    &     0    &     0\\
    8 &  4.0000  &       0      &   0    & 0.0312      &   0     &    0\\
    \hline
\end{tabular}
\end{center}

Para reproducir estos valores utilizar el código 

\begin{verbatim}
function ecuacion_simple() 
    % Numero de puntos
    N = 8; 
    
    % Valores Iniciales
    xini = 0;
    xfin = 4;
    y_ini = 1;
    h = (xfin-xini)/N;
    x = [xini:h:xfin]';
    

    [x1,y1] = Euler(xini, y_ini, h, N, @ecuacion);
    [x2,y2] = Heun(xini, y_ini, h, N, @ecuacion);
    [x3,y3] = Punto_Medio(xini, y_ini, h, N, @ecuacion);
    [x4,y4] = RK(xini, y_ini, h, N, @ecuacion);
    [x5,y5] = RK3(xini, y_ini, h, N, @ecuacion);
    [x6,y6] = RK4(xini, y_ini, h, N, @ecuacion);
    y2
  
    [x, solucion_real(x), y1, y2, y3, y4, y5, y6]

  
      [x, abs(solucion_real(x)-y1),...
          abs(solucion_real(x)-y2),...
          abs(solucion_real(x)-y3),...
          abs(solucion_real(x)-y4),... 
          abs(solucion_real(x)-y5),...
          abs(solucion_real(x)-y6)]
    
    function dy = ecuacion(x, y)
        dy = -2*x^3 + 12*x^2- 20*x + 8.5;
    end

    function y = solucion_real(x)
        y = -0.5*x.^4 + 4*x.^3 - 10*x.^2 + 8.5*x +1;
    end
end
\end{verbatim}

\subsection{Circuito RL}

Para el circuito RL que se muestra en la figura \ref{fig:Cto_RL} calcular el valor de la corriente $i(t)$ en función del tiempo $t$, asumiendo que el interruptor se cierra en el tiempo $t=0$

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/circuitoRL.eps, width=8.0cm}}
   \caption{Circuito RL}
   \label{fig:Cto_RL}
\end{figure}

\textbf{Solución}

La ecuación diferencial para este problema es

$$
\frac{di}{dt} = \frac{V - R i}{L}
$$

La ecuación diferencial en Matlab es:

\begin{verbatim}
    function di = Circuito_RL(t, i)
        %Parametros

        R = 2.0;  % Ohms
        L = 0.3;  % Henries
        V = 10.0; % Volts
        di = (V - R *i)/L;
    end
\end{verbatim}

La solución real de la ecuación diferencial, en función del tiempo,  es:

$$
i(t) =  \frac{V}{R}\left(1 - e^{(-Rt/L)}\right)
$$

La solución utilizando el Método de Runge-Kutta 4 para un valor inicial $t=0$, $i(0)=0$ y un incremento $h=0.2$ es: 

\begin{center}
\begin{tabular}{|c|c|}
    \hline
    $t$ & $i(t)$\\
    \hline
    0.0  &  0\\
    0.2  &  3.5391\\
    0.4  &  4.5732\\
    0.6  &  4.8753\\
    0.8  &  4.9636\\
    1.0  &  4.9894\\
    \hline
 \end{tabular}
\end{center}

En la figura \ref{fig:Sol_Cto_RL} se muestra la solución gráfica de los métodos estudiados en este capítulo. En negro se da la solución real de la ecuación diferencial y en amarillo la solución del Runge-Kutta de cuarto orden que es el que mejor se desempeña para el incremento dado. La peor solución es la calculada con el método de Euler la cual se muestra en azul en esta figura.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/sol_circuitoRL.eps, width=10.0cm}}
   \caption{Circuito RL}
   \label{fig:Sol_Cto_RL}
\end{figure}

El código implementado para este ejemplo es:

\begin{verbatim}
function Circuito_RL()
    
    % Numero de puntos
    N = 100; 
    
    % Valores Iniciales
    tini = 0;
    tfin = 1;
    i_ini = 0;
    h = (tfin-tini)/N;
    

    [t0, i0] = Solucion_Real(tini, i_ini, h, N);
    [t1, i1] = Euler(tini, i_ini, h, N, @Circuito_RL);
    
    [t2, i2] = RK(tini, i_ini, h, N, @Circuito_RL);
    [t3, i3] = RK3(tini, i_ini, h, N, @Circuito_RL);
    [t4, i4] = RK4(tini, i_ini, h, N, @Circuito_RL);
   
    [t5, i5] = Heun(tini, i_ini, h, N, @Circuito_RL);
 

    
    plot(t0, i0, 'ko-', t1, i1, 'b .-', ...
         t2, i2, 'r*-', t3, i3, 'g--',  ...
         t4, i4, 'y+-', t5, i5, 'c');

     xlabel('tiempo (seg.)');
     ylabel('Corriente (Amp.)');
     title('Solucion Circuito RL');

    % Solucion exacta de la ecuacion Diferencial
    
    function [t, ic] =Solucion_Real(t0, ir, h, N) 
       %Parametros

        R = 2.0;  % Ohms
        L = 0.3;  % Henries
        V = 10.0; % Volts
        
        t = [t0:h:N*h]';
        Tau = L/R;
        ic = 5 *(1 - exp(-t/Tau));   
    end

    % Ecuación diferencial del circuito RL

    function di = Circuito_RL(t, i)
        %Parametros

        R = 2.0;  % Ohms
        L = 0.3;  % Henries
        V = 10.0; % Volts
        di = (V - R *i)/L;
    end
end
\end{verbatim}

\subsection{Movimiento parabólico}

El movimiento parabólico lo podemos considerar como un movimiento en dos dimensiones; en una de las direcciones el objeto se desplaza con velocidad contante y en la otra con una aceleración constante

Las ecuaciones diferenciales que rigen este movimiento están dadas como

\begin{equation*}
\begin{array}{c}
\frac{d^2 x}{d t^2} = 0\\
\\
\frac{d^2 y}{d t^2} = a
\end{array}
\end{equation*}


Las ecuaciones diferenciales en el caso del tiro parabólico pueden representarse por un sistema de ecuaciones con cuatro variable el desplazamiento horizontal $x$, el desplazamiento vertical $y$, la velocidad horizontal $v_x$ y la velocidad vertical $v_y$ tal como se muestra a continuación:

\begin{equation*}
\begin{array}{c}
\frac{d^2 x}{d t^2} = 0  \equiv \left \{
	\begin{array}{c}
	     \frac{dx}{dt} = v_x \\
		\\
		\frac{dv_x}{dt} = 0		
	\end{array} \right . \\
\\
\frac{d^2 x}{d t^2} = -g \equiv \left \{
	\begin{array}{c}
		\frac{dy}{dt} = v_y\\
		\\
		\frac{dv_y}{dt} = -g
	\end{array} \right . \\
\end{array}
\end{equation*}

La implementación en Matlab para este conjunto de ecuaciones es:

\begin{verbatim}
function dv = parabolico(t, v)
   g = 9.8; % m/s2
    
   dv(1) = v(2);
   dv(2) = 0;
   dv(3) = v(4);
   dv(4) = -g;
end
\end{verbatim}

Para el caso de el tiro parabólico encontrar la solución del sistema de ecuaciones diferenciales, para un tiempo $t_f = 2.0$ seg,  asumiendo los siguientes inicialmente $x(0) = 0$ mts, $y_0 = 5$ mts, $v_x(0) = 10$ m/s, $v_y(0) = 10$ m/s y un incremento $h= 0.001$.

La solución en matlab es:

\begin{verbatim}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Datos
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


x0 = 0;    % desplazamiento horizontal inicial en mts
y0 = 5;    % desplazamiento vertical inicial en mts
vx0 = 10;  % velocidad horizontal inical en m/seg
vy0 = 10;  % velocidad vertical inicial en m/seg
t0 = 0;    % tiempo inicial seg.
tf = 2;    % tiempo final seg
h = 0.001; % intervalo de integración

N = floor((tf - t0)/h);

[t, d] = RK4(t0, [x0, vx0, y0, vy0], h, N, @parabolico);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Graficas de la solución
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


plot(d(:,1), d(:,3));
grid on
title('Tiro Parabólico');
xlabel('x (mts.)');
ylabel('y (mts.)');

\end{verbatim}

La solución es
 
\begin{verbatim}
>> tiro_parabolico


    1.9970   19.9700   10.0000    5.4288   -9.5706
    1.9980   19.9800   10.0000    5.4192   -9.5804
    1.9990   19.9900   10.0000    5.4096   -9.5902
    2.0000   20.0000   10.0000    5.4000   -9.6000
\end{verbatim}


En el caso particular de un tiro parabólico, donde un objeto es lanzado con una velocidad inicial $v_x(0)$ horizontalmente y $v_y(0)$ verticalmente y esta sujeta a la acción de la gravedad la solución de analítica de las ecuaciones diferenciales es el sistema mostrada a continuación
 
\begin{equation*}
\begin{array}{c}
x(t) = x(0) + v_x(0) t \\
\\
y(t) = y(0) + v_y(0) t - \frac{1}{2} g t^2
\end{array}
\end{equation*} 

La solución de acuerdo con la ecuaciones de movimiento es:

\begin{equation*}
\begin{array}{c}
x(t) = x(0) + v_x(0) t   \\
\\
x(2) = 0 + 10 * 2 = 20\\
\\
y(t) = y(0) + v_y(0) t - \frac{1}{2} g t^2\\
\\
y(2) = 5 + 10*2 - \frac{1}{2} (9.8) 2^2 = 5.40
\end{array}
\end{equation*} 

En la figura \ref{fig:parabolico} se muestra la gráfica del desplazamiento $x(t)$ vs $y(t)$ para la solución del tiro parabólico.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/parabolico.eps, width=12.0cm}}
   \caption{Gráfica de solución del tiro Parabólico}
   \label{fig:parabolico}
\end{figure}

\subsection{Sistema Masa Resorte}

El sistema masa resorte como el que se muestra en la figura \ref{fig:MR}, esta constituido por una masa $M$, sujeta a un resorte de constante $K$, el cual esta sujeto a fuerzas de amortiguamiento $C$ y se aplica una fuerza $F$. El modelo dinámico para este sistema es:

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/MasaResorte.eps, width=12.0cm}}
   \caption{Sistema Masa Resorte}
   \label{fig:MR}
\end{figure}

$$
M \frac{dx^2}{dt^2} + C\frac{dx}{dt} + K x = F
$$ 

Note que tenemos una ecuación diferencial de segundo orden la cual plantearemos como un par de ecuaciones de primer orden como:

$$
\frac{dx}{dt} = v
$$
$$
\frac{dv}{dt} = \frac{ F -   C v - K x}{M}
$$

La implementación de estas funciones en Matlab es:

\begin{verbatim}
    function D = Sistema_MR(t, r)
        M = 10.0;  % Masa  en Kg 
        C = 1;   % Coeficiente de Amortiguamiento Dinnamica
        K = 5;   % Constante del Resorte
        F = 0;    % Fuerza Aplicada
        d = r(1);
        v = r(2);
        
        D(1) = v;   % derivada del desplazamiento
        D(2) = (F - C*v - K*d)/M; % derivada de la velocidad
    end
\end{verbatim}

En la figura \ref{fig:Sol_MR} se muestra la solución de la ecuaciones diferenciales utilizando el método de Runge-Kutta de cuarto orden, para un desplazamiento inicial $x(0) = 2$, una velocidad inicial $v(0) =0$, de una Masa $M= 10$, con un coeficiente de amortiguamiento $C=1$, una constante del resorte $K=5$ y una fuerza externa $F=0$. En la parte superior izquierda tenemos el desplazamiento en función del tiempo a la derecha la velocidad en función del tiempo y en la parte inferior el diagrama de fase que muestra la $v(t)$ vs $x(t)$.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/Sol_MasaResorte.eps, width=12.0cm}}
   \caption{Sistema Masa Resorte}
   \label{fig:Sol_MR}
\end{figure}

La implementación de este ejemplo en Matlab es:

\begin{verbatim}
function Sistema_Masa_Resorte(Metodo)
    
    % Numero de puntos
    N = 1000; 
    
    % Valores Iniciales
    tini = 0;
    tfin = 40;
    d_ini = 2;
    v_ini = 0;
    
    h = (tfin-tini)/N;
    
    [t,x] = Metodo(tini, [d_ini, v_ini], h, N, @Sistema_MR);
        
    title('Sistema Masa Resorte');
  
    subplot(2,2,1)
    plot(t, x(:,1));
    xlabel('tiempo (seg.)');
    ylabel('Desplazamiento (m.)');
    
    subplot(2,2,2)
    plot(t, x(:,2));
    xlabel('tiempo (seg.)');
    ylabel('Velocidad (m/s.)');

    subplot(2,2,3)
    plot(x(:,1), x(:,2));
    xlabel('Desplazamient (m.)');
    ylabel('Velocidad (m/s.)');


    % Ecuación diferencial del circuito RL

    function D = Sistema_MR(t, r)
        M = 10.0;  % Masa  en Kg 
        C = 1;   % Coeficiente de Amortiguamiento Dinnamica
        K = 5;   % Constante del Resorte
        F = 0;    % Fuerza Aplicada
        d = r(1);
        v = r(2);
        
        D(1) = v;   % derivada del desplazamiento
        D(2) = (F - C*v - K*d)/M; % derivada de la velocidad
    end
end
\end{verbatim}

\subsection{Solución Circuito RC}

Considere el circuito serie alimentado por una fuente $V_S$, con una resistencia $R$ y un capacitor $C$ tal como se muestra en la Figura \ref{fig:cot_RC}. Para este circuito la ecuación que lo rige esta dado por:

\[
V_s - R i - \frac{1}{C} \int i dt = 0 
\]

Para resolver el circuito, vamos a sustituir el valor de la corriente como el cambio de la carga con respecto al tiempo $i=\frac{dq} {dt}$. 

\[
V_s - R \frac{dq}{dt} - \frac{1}{C} q = 0 
\]


\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/CircuitoRC.jpg, width=12.0cm}}
   \caption{Circuito RC}
   \label{fig:cto_RC}
\end{figure}

Para llevar a cabo la integración numérica vamos a considerar que $V_s = 10$, $R=2$ y $C=0.01$. La función para llevar a cabo la integración numérica esta dada por:

\begin{verbatim}
    function dq = CtoRC(t, q)
        %Parametros
        R = 2.0;    % Ohms
        C = 0.01;   % Farads
        V = 10.0;   % Volts
        dq = (V*C - q)/(R*C);
    end
\end{verbatim}

Para calcular la corriente que circula a través del circuito es necesario calcular la derivada, la cual podemos  calcular mediante la siguiente expresión:

\[
i = \frac{d q}{dt } \approx \frac{q(h) - q(h-1)}{h}
\]

Para hacer la implementación en Matlab haremos

\begin{verbatim}
      ic = (q(2:N)-q(1:N-1))/h;
      tc = [t0:h:tf];
      tc = tc(1:N-1);
\end{verbatim}

El código completo para hacer la simulación se muestra a continuación y en la Figura \ref{fig:sol_RC} los resultados de la solución de la ecuación diferencial. Para esta Figura en la parte de arriba se muestra la solución de la carga con respecto al tiempo y en la parte de abajo la corriente en función del tiempo.

\begin{verbatim}
function Circuito_RC()
    
    % Numero de puntos
    N = 10000; 
    
    % Valores Iniciales
    t0 = 0;
    tf = 0.13;
    q0 = 0;
    h = (tf-t0)/N;

    [t, q] = RK4(t0, q0, h, N, @CtoRC);
    
    subplot(2,1,1);
    plot(t, q(:,1));
    xlabel('tiempo (seg.)');
    ylabel('Carga (Coulombs.)');
    title('Solucion Circuito RC');

    % Calculo de la derivada 
    ic = (q(2:N)-q(1:N-1))/h;
    tc = [t0:h:tf];
    tc = tc(1:N-1);
    

    subplot(2,1,2);
    plot(tc, ic);
    xlabel('tiempo (seg.)');
    ylabel('Corriente (Amp.)');
    title('Solucion Circuito RC');

    % Ecuación diferencial del circuito RL

    function dq = CtoRC(t, q)
        %Parametros
        R = 2.0;    % Ohms
        C = 0.01;   % Farads
        V = 10.0;   % Volts
        dq = (V*C - q)/(R*C);
    end
end

\end{verbatim}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_8/figuras/sol_RC.jpg, width=12.0cm}}
   \caption{Solución numérica del circuito RC. En la parte de arriba la carga con respecto al tiemo y abajo la corriente en función del tiempo.}
   \label{fig:sol_RC}
\end{figure}