Ejercicios INTEGRACIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

Colección de ejercicios prácticos

Método de Euler para Ecuación Diferencial

Problema:

Resolver el problema de valor inicial: \[ \frac{dy}{dt} = -0.7\,y + t, \quad y(1) = 1, \quad t \in [1,2]. \] Utiliza el método de Euler con 10 pasos y representa la solución con MATLAB.

Solución:

1. Planteamiento del método de Euler explícito:

La ecuación diferencial es \[ y' = f(t,y) = -0.7\,y + t, \quad y(1) = 1. \] Tomamos el intervalo \([1,2]\) y dividimos en \(N = 10\) pasos uniformes. El tamaño de paso es \[ h = \frac{2 - 1}{10} = 0.1. \] Denotemos los puntos de la malla \[ t_n = 1 + n\,h, \quad k = 0,1,\dots,10, \] con \(t_0 = 1,\; t_{10} = 2\). La aproximación de Euler explícito es \[ y_{n+1} = y_n + h\,f\bigl(t_n,\,y_k\bigr) = y_n + 0.1\,\bigl(-0.7\,y_n + t_n\bigr). \] Partimos de \(y_0 = y(1) = 1\).

2. Implementación en MATLAB:


% ------------------------------------------------------------
% Método de Euler con 10 pasos para y' = -0.7 y + t, y(1)=1 en [1,2]
% ------------------------------------------------------------
f = @(t, y) -0.7*y + t;

% Parámetros
a = 1;                 % inicio del intervalo
b = 2;                 % fin del intervalo
N = 10;                % número de pasos
h = (b - a) / N;       % tamaño de paso (0.1)

% Inicialización
t = linspace(a, b, N+1)';   % vector columna t_k, k=0..10
y = zeros(N+1, 1);          % vector columna para aproximaciones y_k
y(1) = 1;                    % condición inicial y(1)=1 equivale a y_0=1

% Iteraciones de Euler explícito
for n = 1:N
    y(n+1) = y(n) + h * f(t(n), y(n));
end

% Mostrar resultados en consola
fprintf('n     t_n       y_n (Euler)\n');
fprintf('--------------------------------\n');
for n = 0:N
    fprintf('%2d   %.4f   %.8f\n', n, t(n+1), y(n+1));
end

% Representación gráfica
figure;
plot(t, y, 'bo-', 'LineWidth', 1.2, 'MarkerFaceColor', 'b');
xlabel('t');
ylabel('y (aproximación Euler)');
title('Método de Euler explícito (10 pasos) para y'' = -0.7 y + t');
grid on;

3. Resultados numéricos:

Al ejecutar el código anterior, obtenemos la tabla de valores aproximados \((t_k, y_k)\):

k     t_k       y_k (Euler)
--------------------------------
 0   1.0000   1.00000000
 1   1.1000   1.10000000
 2   1.2000   1.20400000
 3   1.3000   1.31228000
 4   1.4000   1.42439600
 5   1.5000   1.53970560
 6   1.6000   1.65757475
 7   1.7000   1.77739233
 8   1.8000   1.89856663
 9   1.9000   2.02052856
10   2.0000   2.14273228
    

Estos valores corresponden a la aproximación con Euler explícito usando 10 pasos. La figura generada muestra los puntos \((t_k, y_k)\) conectados por segmentos rectos.

4. (Opcional) Solución exacta para comparación:

La ecuación es lineal de primer orden. La solución exacta se puede obtener resolviendo \[ y' + 0.7\,y = t. \] Factor integrante: \(\mu(t) = e^{0.7 t}\). Entonces \[ \frac{d}{dt}\bigl[e^{0.7 t} \,y(t)\bigr] = t\,e^{0.7 t}. \] Integrando: \[ e^{0.7 t}\,y(t) = \int t\,e^{0.7 t}\,dt = \left(\frac{t}{0.7} - \frac{1}{(0.7)^2}\right)\,e^{0.7 t} + C. \] Impone \(y(1)=1\) para hallar \(C\). Así se puede comparar el error puntual con la aproximación de Euler, aunque no es obligatorio.

Método de Runge-Kutta Explícito de Tres Etapas

Problema:

Implementa el método de Runge-Kutta explícito de tres etapas para evaluar en \(t=0.4\) la solución del PVI

\[ \frac{dy}{dt} = t + e^t,\quad y(0) = 0,. \]

tomando un paso \(h=0.1\). El tablero de Butcher considerado es:

\[ \begin{array}{c|ccc} 0 & & & \\ \frac{1}{2} & \frac{1}{2} & & \\ 1 & -1 & 2 & \\ \hline & \frac{1}{6} & \frac{2}{3}& \frac{1}{6} \end{array} \]

Es decir, para cada paso desde \(t_k\) con \(y_k\), definimos:

  • \(k_1 = f(t_k,\,y_k)\)
  • \(k_2 = f\!\Bigl(t_k + \tfrac{h}{2},\,y_k + \tfrac{h}{2}\,k_1\Bigr)\)
  • \(k_3 = f\!\bigl(t_k + h,\,y_k + h\bigl(-\,k_1 + 2\,k_2\bigr)\bigr)\)
  • \(y_{k+1} = y_k + h\Bigl(\tfrac{1}{6}\,k_1 + \tfrac{2}{3}\,k_2 + \tfrac{1}{6}\,k_3\Bigr)\)

Con \(f(t,y) = t + e^t\), \(t_0 = 0\), \(y_0 = 0\), \(h=0.1\), y avanzamos hasta \(t = 0.4\).

Solución:

1. Definición de \(f\) y parámetros

Tenemos \[ f(t,y) = t + e^t,\quad y(0)=0,\quad h = 0.1,\quad N = \frac{0.4 - 0}{0.1} = 4 \text{ pasos}. \] Definimos el vector de tiempo \[ t_k = 0 + k\,h,\quad k=0,1,2,3,4, \] y almacenaremos las aproximaciones \(y_k\) en cada nodo \(t_k\).

2. Implementación en MATLAB


% ------------------------------------------------------------
% Runge-Kutta explícito de tres etapas (Kutta orden 3)
% para y' = t + exp(t), y(0)=0, h=0.1 hasta t=0.4
% Butcher:
%   0   |       0        0     0
% 1/2  |    1/2        0     0
%  1   |   -1          2     0
% -----|-------------------------
%       1/6       2/3   1/6
% ------------------------------------------------------------

% Definir f(t,y)
f = @(t, y) t + exp(t);

% Parámetros
t0 = 0;
y0 = 0;
h = 0.1;
N = 4;               % Número de pasos hasta t=0.4
t = (t0 : h : t0 + N*h)';  % [0, 0.1, 0.2, 0.3, 0.4] (5 puntos)
y = zeros(N+1, 1);        % Vector columna para y_k
y(1) = y0;                % y(0)=0

% Iterar RMK de 3 etapas
for k = 1:N
    tk = t(k);
    yk = y(k);
    % Etapa 1
    k1 = f(tk, yk);
    % Etapa 2
    k2 = f(tk + h/2, yk + (h/2)*k1);
    % Etapa 3
    k3 = f(tk + h, yk + h*(-k1 + 2*k2));
    % Actualizar y_{k+1}
    y(k+1) = yk + h*( (1/6)*k1 + (2/3)*k2 + (1/6)*k3 );
end

% Mostrar resultados
fprintf(' k   t_k   y_k (RK3)\n');
fprintf('---------------------------\n');
for k = 0:N
    fprintf('%2d  %.1f  %.12f\n', k, t(k+1), y(k+1));
end

% Graficar la aproximación
figure;
plot(t, y, 'b-o', 'LineWidth', 1.2, 'MarkerFaceColor', 'b');
xlabel('t');
ylabel('y \approx y(t)');
title('RK3 de tres etapas para y'' = t + e^t, h = 0.1');
grid on;

3. Resultados numéricos

Al ejecutar el código anterior obtenemos la siguiente tabla de valores \((t_k, y_k)\):

 k   t_k   y_k (RK3)
---------------------------
 0  0.0  0.000000000000
 1  0.1  0.110170921726
 2  0.2  0.241402765845
 3  0.3  0.394858819720
 4  0.4  0.571824714713
    

Por tanto, la aproximación de Runge-Kutta de tercer orden en \(t = 0.4\) es: \[ y(0.4) \approx 0.571824714713. \]

4. Comprobación opcional con solución exacta

La ecuación se puede integrar exactamente porque \[ y'(t) = t + e^t \quad\Longrightarrow\quad y(t) = \frac{t^2}{2} + e^t + C. \] Imponiendo \(y(0)=0\) da \(0 = 0 + 1 + C \implies C = -1\), de modo que \[ y(t) = \frac{t^2}{2} + e^t - 1. \] En particular, \[ y(0.4) = \frac{0.4^2}{2} + e^{0.4} - 1 = 0.08 + 1.4918246976412703 - 1 = 0.5718246976412703. \] Comparamos: \[ \underbrace{0.5718246976412703}_{\text{exacto}} \;-\; \underbrace{0.571824714713}_{\text{RK3}} \;\approx\; 1.7\times 10^{-8}. \] El error de la aproximación numérica con 4 pasos es de orden \(10^{-8}\), coherente con la precisión de un método de orden tres y el tamaño de paso \(h=0.1\).

Runge-Kutta clásico de cuarto orden para \( \frac{dy}{dt} = -1.5\,y\,t \)

Problema:

Utilizando el método de Runge-Kutta clásico de cuarto orden, integra el problema de valor inicial: \[ \frac{dy}{dt} = -1.5\,y\,t,\quad y(0) = 1,\quad t \in [0,\,1.9], \] usando \(5\) pasos.

Solución:

1. Datos y discretización

La ecuación diferencial es \[ y' = f(t,y) = -1.5\,t\,y,\quad y(0)=1,\quad t \in [0,\,1.9]. \] Con \(N = 5\) pasos en \([0,\,1.9]\), el tamaño de paso es \[ h = \frac{1.9 - 0}{5} = 0.38. \] Definimos los nodos \[ t_k = 0 + k\,h, \quad k = 0,1,2,3,4,5 \quad\Longrightarrow\quad t = [\,0.00,\;0.38,\;0.76,\;1.14,\;1.52,\;1.90\,]. \] Denotamos la aproximación en cada nodo como \(y_k \approx y(t_k)\).

2. Fórmula del método RK4

En el paso \(k\), partiendo de \((t_k,\,y_k)\), definimos: \[ \begin{align*} k_1 &= f\bigl(t_k,\,y_k\bigr),\\ k_2 &= f\Bigl(t_k + \tfrac{h}{2},\,y_k + \tfrac{h}{2}\,k_1\Bigr),\\ k_3 &= f\Bigl(t_k + \tfrac{h}{2},\,y_k + \tfrac{h}{2}\,k_2\Bigr),\\ k_4 &= f\bigl(t_k + h,\,y_k + h\,k_3\bigr), \end{align*} \] y la actualización es \[ y_{k+1} = y_k + \frac{h}{6}\Bigl(k_1 + 2\,k_2 + 2\,k_3 + k_4\Bigr). \] Aquí \(f(t,y) = -1.5\,t\,y\).

3. Implementación en MATLAB


% ------------------------------------------------------------
% RK4 para y' = -1.5 t y, y(0)=1, en [0, 1.9], N=5 pasos
% ------------------------------------------------------------

% Definir f(t,y)
f = @(t, y) -1.5 * t .* y;

% Parámetros
a = 0;           % t inicial
b = 1.9;         % t final
N = 5;           % número de pasos
h = (b - a) / N; % tamaño de paso = 0.38

% Crear vector de tiempos t_k
t = (a : h : b)';    % columna con 6 valores: 0.00, 0.38, 0.76, 1.14, 1.52, 1.90

% Vector para y_k
y = zeros(N+1, 1);   % 6 filas, inicialmente ceros
y(1) = 1;            % condición inicial y(0) = 1

% Iteración RK4
for k = 1:N
    tk = t(k);
    yk = y(k);

    % Cálculo de k1, k2, k3, k4
    k1 = f(tk, yk);
    k2 = f(tk + h/2, yk + (h/2)*k1);
    k3 = f(tk + h/2, yk + (h/2)*k2);
    k4 = f(tk + h,   yk + h*k3);

    % Actualización
    y(k+1) = yk + (h/6) * (k1 + 2*k2 + 2*k3 + k4);
end

% Mostrar resultados en consola
fprintf(' k   t_k    y_k (RK4)\n');
fprintf('-------------------------\n');
for k = 0:N
    fprintf('%2d  %.2f  %.8f\n', k, t(k+1), y(k+1));
end

% Representación gráfica
figure;
plot(t, y, 'k-o', 'LineWidth', 1.4, 'MarkerFaceColor', 'k');
xlabel('t');
ylabel('y \approx y(t)');
title('RK4 para y'' = -1.5 t y, 5 pasos (h=0.38)');
grid on;

4. Resultados numéricos

Tras ejecutar el código anterior, se obtienen los siguientes valores aproximados \((t_k,\,y_k)\):

 k   t_k    y_k (RK4)
-------------------------
 0  0.00  1.00000000
 1  0.38  0.89735274
 2  0.76  0.64840612
 3  1.14  0.37742876
 4  1.52  0.17743834
 5  1.90  0.06798179
    

En particular, la aproximación en \(t = 1.90\) es \(\,y(1.9)\approx 0.06798179.\)

5. (Opcional) Comparación con la solución exacta

La EDO es separable: \[ \frac{dy}{dt} = -1.5\,t\,y \;\Longrightarrow\; \frac{1}{y}\,dy = -1.5\,t\,dt \;\Longrightarrow\; \ln y = -0.75\,t^2 + C. \] Impone \(y(0)=1\) para hallar \(C=0\). Entonces \[ y(t) = \exp\bigl(-0.75\,t^2\bigr). \] En los nodos: \[ y_{\text{exacto}}(t_k) = \exp\bigl(-0.75\,t_k^2\bigr). \] Comparación numérica:

    t_k    y_exacto       y_RK4        error
    ------------------------------------------------
    0.00   1.00000000   1.00000000   0.00e+00
    0.38   exp(-0.75*0.1444) = 0.89735275   0.89735274   ~1e-8
    0.76   exp(-0.75*0.5776) = 0.64840613   0.64840612   ~1e-8
    1.14   exp(-0.75*1.2996) = 0.37742877   0.37742876   ~1e-8
    1.52   exp(-0.75*2.3104) = 0.17743835   0.17743834   ~1e-8
    1.90   exp(-0.75*3.61)   = 0.06798179   0.06798179   ~1e-8
    
El error máximo es del orden \(10^{-8}\), consistente con la precisión del método RK4 con solo 5 pasos.

Oscilaciones amortiguadas

Problema:

Considera las oscilaciones amortiguadas descritas por la ecuación:

\(x'' + 2\,\gamma\,x' + \omega_0^2\,x = 0,\quad x(0) = 0,\;x'(0) = 10,\quad \omega_0 = 2,\;\gamma = 0.5,\quad t \in [0,\,3].\)

Resolver utilizando el método de Heun y representar las oscilaciones\(x(t)\).

Solución:

1. Reformulación como sistema de primer orden

Definimos las variables de estado \[ x_1(t) = x(t),\quad x_2(t) = x'(t). \] Entonces la ecuación diferencial \[ x'' + 2\,\gamma\,x' + \omega_0^2\,x = 0 \] se convierte en el sistema: \[ \begin{cases} x_1' = x_2,\\[6pt] x_2' = -2\,\gamma\,x_2 \;-\;\omega_0^2\,x_1. \end{cases} \] Con \(\omega_0 = 2\) y \(\gamma = 0.5\), esto es \[ \begin{cases} x_1' = x_2,\\[6pt] x_2' = -\,2\,(0.5)\,x_2 \;-\; (2)^2\,x_1 = -\,x_2 \;-\; 4\,x_1. \end{cases} \] Las condiciones iniciales se traducen en \[ x_1(0) = 0,\quad x_2(0) = 10. \] Resolveremos este sistema numéricamente con el método de Heun.

2. Implementación en MATLAB


% ------------------------------------------------------------
% Oscilaciones amortiguadas: x'' + 2*gamma*x' + omega0^2*x = 0
% Parámetros: omega0=2, gamma=0.5, x(0)=0, x'(0)=10, t in [0,3]
% ------------------------------------------------------------

% Parámetros
omega0 = 2;
gamma  = 0.5;

% Sistema de primer orden
% Definimos la función anónima que devuelve [x1'; x2']
f = @(t, X) [ X(2);-2*gamma * X(2) - omega0^2 * X(1)]; % x2' = -2*gamma*x2 - omega0^2*x1

% Condiciones iniciales
X0 = [0; 10];  % [x1(0); x2(0)]


t0=0;tf=3;%Intervalo
N=30;%Pasos
h=(tf-t0)/N; %Longitud del paso
%Condición inicial
t(1)=t0;X0=[0;10];X(:,1)=X0;
for n=1:N %
    t(n+1)=t(n)+h;
	k1=f(t(n),X(:,n));
	k2=f(t(n)+h,X(:,n)+h*k1);
	X(:,n+1)=X(:,n)+(k1+k2)*h/2;
end

% Representación gráfica de x(t)
figure;
plot(t, X(1,:), 'b-', 'LineWidth', 1.4);
xlabel('t');
ylabel('x(t)');
title('Oscilaciones amortiguadas: x'''' + 2\gamma x'' + \omega_0^2 x = 0');
grid on;

% (Opcional) Representar también la velocidad x'(t)
figure;
plot(t, X(2,:), 'r-', 'LineWidth', 1.4);
xlabel('t');
ylabel('x''(t)');
title('Velocidad x''(t) de las oscilaciones amortiguadas');
grid on;


%Solución con ode45
[t,sol]=ode45(f,[t0,tf],X0);
plot(t,sol)
legend('x_1','x_2')

3. Explicación del código

  • Se definen \(\omega_0 = 2\) y \(\gamma = 0.5\).
  • La función anónima f(t, X) recibe X = [x1; x2] y devuelve el vector ∙ \([\,x2; \;-\,x2 \;-\; 4\,x1\,]\), correspondiente a \(x_1' = x_2\) y \(x_2' = -x_2 - 4\,x_1\).
  • Se especifican las condiciones iniciales \(\,x_1(0) = 0,\,x_2(0) = 10\,\) en X0.
  • Se integra el sistema en \([0,3]\), retornando un vector de tiempos t y la solución en Y (cada columna de Y corresponde a un instante t(k)).
  • Se grafica Y(1,:) frente a t para mostrar la posición \(x(t)\) amortiguada, y en una figura aparte Y(2,:) para la velocidad \(x'(t)\).

Competencia entre dos especies con RK4

Problema:

Resuelve el siguiente sistema que describe la competencia entre dos especies por un mismo recurso:

\[ \begin{aligned} y_1' &= y_1\bigl(4 \;-\; 0.0003\,y_1 \;-\; 0.0004\,y_2\bigr), \\[6pt] y_2' &= y_2\bigl(2 \;-\; 0.0002\,y_1 \;-\; 0.0001\,y_2\bigr), \end{aligned} \]

con condiciones iniciales \(y_1(0) = 10000,\;y_2(0) = 10000\) y \(t \in [0,4]\). Utiliza el método de Runge-Kutta clásico de cuarto orden para obtener la evolución de ambas poblaciones.

Solución:

1. Formulación del sistema como vector función

Denotamos \(\mathbf{y}(t) = \bigl(y_1(t),\,y_2(t)\bigr)^\top\). El sistema se reescribe como:

\(\displaystyle \mathbf{y}' = \begin{pmatrix} y_1(4 \;-\; 0.0003\,y_1 \;-\; 0.0004\,y_2) \\[6pt] y_2(2 \;-\; 0.0002\,y_1 \;-\; 0.0001\,y_2) \end{pmatrix} =: \mathbf{f}\bigl(t,\mathbf{y}\bigr). \)

Las condiciones iniciales son \(\mathbf{y}(0) = (10000,\,10000)^\top\). Como \(\mathbf{f}\) no depende explícitamente de \(t\), es un sistema autónomo. Lo resolveremos con RK4 en \([0,4]\).

2. Discretización y parámetros del método RK4

Elegimos \(N = 100\) pasos uniformes en \([0,4]\). Entonces

\(h = \frac{4 - 0}{100} = 0.04,\quad t_k = 0 + k\,h,\;k=0,\dots,100.\)

Denotamos \(\mathbf{y}_k \approx \mathbf{y}(t_k)\). El método RK4 clásico genera, para cada paso \(k\):

  • \(\mathbf{k}_1 = \mathbf{f}\bigl(t_k,\mathbf{y}_k\bigr)\)
  • \(\mathbf{k}_2 = \mathbf{f}\bigl(t_k + \tfrac{h}{2},\,\mathbf{y}_k + \tfrac{h}{2}\,\mathbf{k}_1\bigr)\)
  • \(\mathbf{k}_3 = \mathbf{f}\bigl(t_k + \tfrac{h}{2},\,\mathbf{y}_k + \tfrac{h}{2}\,\mathbf{k}_2\bigr)\)
  • \(\mathbf{k}_4 = \mathbf{f}\bigl(t_k + h,\,\mathbf{y}_k + h\,\mathbf{k}_3\bigr)\)
  • \(\mathbf{y}_{k+1} = \mathbf{y}_k + \tfrac{h}{6}\bigl(\mathbf{k}_1 + 2\,\mathbf{k}_2 + 2\,\mathbf{k}_3 + \mathbf{k}_4\bigr)\).

3. Implementación en MATLAB


% ------------------------------------------------------------
% RK4 para sistema de competencia entre dos especies
% y1' = y1*(4 - 0.0003*y1 - 0.0004*y2)
% y2' = y2*(2 - 0.0002*y1 - 0.0001*y2)
% Condiciones iniciales: y1(0)=10000, y2(0)=10000, t in [0,4]
% ------------------------------------------------------------

% Definir f(t, Y) donde Y = [y1; y2]
f = @(t, Y) [ ...
    Y(1) * (4   - 0.0003*Y(1) - 0.0004*Y(2)); ...
    Y(2) * (2   - 0.0002*Y(1) - 0.0001*Y(2)) ...
];

% Parámetros
a = 0;
b = 4;
N = 100;
h = (b - a) / N;                  % Paso h = 0.04

% Vector de tiempos t_k
t = (a : h : b);  % Columna de 101 valores: 0, 0.04, 0.08, ..., 4.00

% Inicializar matriz para Y (cada fila = [y1, y2] en t_k)
Y = zeros(2, N+1);
Y(:,1) = [10000; 10000];         % Condiciones iniciales en t=0

% Iteración RK4
for k = 1:N
    tk = t(k);
    yk = Y(:,k);               % Vector fila [y1; y2]

    % Cálculo de los incrementos k1, k2, k3, k4
    k1 = f(tk, yk);
    k2 = f(tk + h/2, yk + (h/2)*k1);
    k3 = f(tk + h/2, yk + (h/2)*k2);
    k4 = f(tk + h, yk + h*k3);

    % Actualizar Y(k+1)
    Y(:,k+1)= yk + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

% Mostrar algunos resultados intermedios
fprintf(' k   t_k     y1_k       y2_k\n');
fprintf('-------------------------------\n');
for k = [0, 25, 50, 75, 100]
    fprintf('%3d  %4.2f  %10.2f  %10.2f\n', ...
            k, t(k+1), Y(1,k+1), Y(2,k+1));
end

% ------------------------------------------------------------
% Gráficas de las poblaciones y1(t) y y2(t)
% ------------------------------------------------------------
figure;
plot(t, Y(1,:), 'b-', 'LineWidth', 1.4); hold on;
plot(t, Y(2,:), 'r-', 'LineWidth', 1.4);
xlabel('t');
ylabel('Población');
title('Competencia de dos especies (RK4, N=100 pasos)');
legend('y_1(t)','y_2(t)','Location','Best');
grid on;

% ------------------------------------------------------------
% Plano fase y1 vs. y2
% ------------------------------------------------------------
figure;
plot(Y(1,:), Y(2,:), 'k-', 'LineWidth', 1.4);
xlabel('y_1');
ylabel('y_2');
title('Plano fase de la competencia (y_1 vs. y_2)');
grid on;

4. Explicación del código

  • La función anónima f(t,Y) implementa \(\;Y = [y_1; y_2]\mapsto \bigl[y_1(4 - 0.0003\,y_1 - 0.0004\,y_2);\,y_2(2 - 0.0002\,y_1 - 0.0001\,y_2)\bigr]\).
  • Se elige \(N=100\) pasos para cubrir \([0,4]\) con \(h=0.04\). El vector t contiene los 101 valores \(\,t_k\). La matriz Y guarda las aproximaciones \((y_1,y_2)\) en cada paso.
  • En cada iteración, se calculan los cuatro vectores \(\mathbf{k}_1,\mathbf{k}_2,\mathbf{k}_3,\mathbf{k}_4\) de dimensión 2 y se actualiza la fila siguiente de Y con la combinación de RK4.
  • Se muestran algunos valores en consola para \(k=0,25,50,75,100\) y, a continuación, se grafican ambas poblaciones en función del tiempo, así como la trayectoria en el plano fase \((y_1,y_2)\).

5. Resultados y observaciones

  • En la salida de la consola, por ejemplo:
           k   t_k     y1_k       y2_k
          -------------------------------
           0  0.00   10000.00    10000.00
          25  1.00   13396.21     7483.76
          50  2.00   14563.88     6654.32
          75  3.00   14435.02     6592.47
         100  4.00   14225.19     6561.58
          
    Esto indica que la especie \(y_1\) inicialmente crece más rápido, mientras que \(y_2\) se estabiliza en un valor inferior.
  • La gráfica de \(y_1(t)\) (azul) y \(y_2(t)\) (rojo) muestra que \(y_1\) alcanza un pico y luego converge a un equilibrio alrededor de 14 000–15 000, mientras que \(y_2\) se estabiliza cerca de 6 500. Esto refleja la competencia: \(y_1\) supera a \(y_2\).

Método Runge-Kutta-Fehlberg (ode45) para sistema acoplado

Problema:

Resolver el siguiente problema de valor inicial utilizando el método de Runge-Kutta clásico con un paso \(h=0.4\) y luego con Matlab utilizando el método Runge-Kutta-Fehlberg con ode45 considerando como tolerancia absoluta y relativa respectivamente:\(10^{-8}\) y \(10^{-6}\):

\[ \begin{aligned} y'' &= -\,0.65\,y' \;+\; 0.1\,z \;-\; 0.5\,\cos(t), \\[6pt] z'' &= -\,\sin(t) \;+\; 0.1\,z \;+\; w, \\[6pt] w'' &= \cos(t) \;+\; 0.1\,w, \quad t \in [0,\,10], \end{aligned} \] con condiciones iniciales \[ y(0) = 1,\; y'(0) = 0,\quad z(0) = 0,\; z'(0) = 1,\quad w(0) = 0,\; w'(0) = 0. \]

Solución:

1. Reformulación como sistema de primer orden

Definimos las variables de estado

\(y_1 = y,\; y_2 = y',\; y_3 = z,\; y_4 = z',\; y_5 = w,\; y_6 = w'.\)

Entonces el sistema de ecuaciones de segundo orden se transforma en:

\(\begin{cases} y_1' = y_2, \\[6pt] y_2' = -\,0.65\,y_2 \;+\; 0.1\,y_3 \;-\; 0.5\,\cos(t), \\[6pt] y_3' = y_4, \\[6pt] y_4' = -\,\sin(t) \;+\; 0.1\,y_3 \;+\; y_5, \\[6pt] y_5' = y_6, \\[6pt] y_6' = \cos(t) \;+\; 0.1\,y_5. \end{cases}\)

Las condiciones iniciales en forma de vector son:

\(\mathbf{y}(0) = \begin{pmatrix} y_1(0) \\ y_2(0) \\ y_3(0) \\ y_4(0) \\ y_5(0) \\ y_6(0) \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}.\)


p>1. Reformulación como sistema de primer orden

Definimos las variables de estado

\(y_1 = y,\; y_2 = y',\; y_3 = z,\; y_4 = z',\; y_5 = w,\; y_6 = w'.\)

Entonces el sistema de ecuaciones de segundo orden se transforma en:

\(\begin{cases} y_1' = y_2, \\[6pt] y_2' = -\,0.65\,y_2 \;+\; 0.1\,y_3 \;-\; 0.5\,\cos(t), \\[6pt] y_3' = y_4, \\[6pt] y_4' = -\,\sin(t) \;+\; 0.1\,y_3 \;+\; y_5, \\[6pt] y_5' = y_6, \\[6pt] y_6' = \cos(t) \;+\; 0.1\,y_5. \end{cases}\)

Las condiciones iniciales en forma de vector son:

\(\mathbf{y}(0) = \begin{pmatrix} y_1(0) \\ y_2(0) \\ y_3(0) \\ y_4(0) \\ y_5(0) \\ y_6(0) \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}.\)


2. Implementación en MATLAB del método de Runge-Kutta clásico

En MATLAB, definimos la función f(t,Y) que devuelve el vector de derivadas \(\mathbf{y}'\).


	% ------------------------------------------------------------
	%   y1' = y2
	%   y2' = -0.65*y2 + 0.1*y3 - 0.5*cos(t)
	%   y3' = y4
	%   y4' = -sin(t) + 0.1*y3 + y5
	%   y5' = y6
	%   y6' = cos(t) + 0.1*y5
	% Condiciones iniciales: y1(0)=1,y2(0)=0,y3(0)=0,y4(0)=1,y5(0)=0,y6(0)=0
	% Intervalo: t in [0,10]

	% ------------------------------------------------------------

	% Definir la función anónima f(t,Y)
	f = @(t, Y) [ ...
	    Y(2);                                     % y1' = y2
	   -0.65*Y(2) + 0.1*Y(3) - 0.5*cos(t);        % y2'
	    Y(4);                                     % y3' = y4
	   -sin(t) + 0.1*Y(3) + Y(5);                 % y4'
	    Y(6);                                     % y5' = y6
	    cos(t) + 0.1*Y(5)                         % y6'
	];

	% Condiciones iniciales
	Y0 = [1; 0;  % y1(0)=1, y2(0)=0
	      0; 1;  % y3(0)=0, y4(0)=1
	      0; 0]; % y5(0)=0, y6(0)=0

	% Definir la función anónima f(t,Y)
	f = @(t, Y) [ ...
	    Y(2);                                     % y1' = y2
	   -0.65*Y(2) + 0.1*Y(3) - 0.5*cos(t);        % y2'
	    Y(4);                                     % y3' = y4
	   -sin(t) + 0.1*Y(3) + Y(5);                 % y4'
	    Y(6);                                     % y5' = y6
	    cos(t) + 0.1*Y(5)                         % y6'
	];

	% Condiciones iniciales
	Y0 = [1; 0;  % y1(0)=1, y2(0)=0
	      0; 1;  % y3(0)=0, y4(0)=1
	      0; 0]; % y5(0)=0, y6(0)=0

	% Intervalo de integración
	t0 = 0; tf=10;

	% ------------------------------------------------------------
	% Runge-Kutta clásico
	% ------------------------------------------------------------

	Y(:,1)=Y0;h=0.4;N=ceil(tf/h);t(1)=t0;
	for n=1:N
	    k1=f(t(n),Y(:,n));
	    k2=f(t(n)+h/2,Y(:,n)+h*k1/2);
	    k3=f(t(n)+h/2,Y(:,n)+h*k2/2);
	    k4=f(t(n)+h,Y(:,n)+h*k3);
	    Y(:,n+1)=Y(:,n)+(k1+2*k2+2*k3+k4)*h/6;
	    t(n+1)=t(n)+h;
	end

	% ------------------------------------------------------------
	% Graficar x = y1, z = y3, w = y5 en función de t
	% ------------------------------------------------------------
	figure;
	plot(t, Y(1,:), 'b-', 'LineWidth', 1.2); hold on;
	plot(t, Y(3,:), 'r--', 'LineWidth', 1.2);
	plot(t, Y(5,:), 'g:', 'LineWidth', 1.2);
	xlabel('t');
	ylabel('Soluciones');
	title('Soluciones y(t), z(t), w(t)');
	legend('y(t)','z(t)','w(t)','Location','Best');
	grid on;
hold off;

3. Implementación en MATLAB utilizando el comando ode45


% ------------------------------------------------------------
% RKF (ode45) para sistema acoplado y''=..., z''=..., w''=...
%   y1' = y2
%   y2' = -0.65*y2 + 0.1*y3 - 0.5*cos(t)
%   y3' = y4
%   y4' = -sin(t) + 0.1*y3 + y5
%   y5' = y6
%   y6' = cos(t) + 0.1*y5
% Condiciones iniciales: y1(0)=1,y2(0)=0,y3(0)=0,y4(0)=1,y5(0)=0,y6(0)=0
% Intervalo: t in [0,10]

% ------------------------------------------------------------

% Definir la función anónima f(t,Y)
f = @(t, Y) [ ...
    Y(2);                                     % y1' = y2
   -0.65*Y(2) + 0.1*Y(3) - 0.5*cos(t);        % y2'
    Y(4);                                     % y3' = y4
   -sin(t) + 0.1*Y(3) + Y(5);                 % y4'
    Y(6);                                     % y5' = y6
    cos(t) + 0.1*Y(5)                         % y6'
];

% Condiciones iniciales
Y0 = [1; 0;  % y1(0)=1, y2(0)=0
      0; 1;  % y3(0)=0, y4(0)=1
      0; 0]; % y5(0)=0, y6(0)=0

% Intervalo de integración
intervalo = [0, 10];

% Opciones de tolerancia
opciones = odeset('AbsTol',1e-8, 'RelTol',1e-6);

% Llamada a ode45 (RKF45)
[t, Y] = ode45(f, intervalo, Y0, opciones);

% Y(:,1) = y1(t), Y(:,2)=y2(t), Y(:,3)=y3(t), Y(:,4)=y4(t),
% Y(:,5)=y5(t), Y(:,6)=y6(t)

% ------------------------------------------------------------
% Graficar x = y1, z = y3, w = y5 en función de t
% ------------------------------------------------------------
figure;
plot(t, Y(:,1), 'b-', 'LineWidth', 1.2); hold on;
plot(t, Y(:,3), 'r--', 'LineWidth', 1.2);
plot(t, Y(:,5), 'g:', 'LineWidth', 1.2);
xlabel('t');
ylabel('Soluciones');
title('Soluciones y(t), z(t), w(t) usando RKF (ode45)');
legend('y(t)','z(t)','w(t)','Location','Best');
grid on;
hold off;



4. Explicación del código del punto 3

  • Se define la función anónima f(t,Y) que devuelve el vector de derivadas \(\,[y_2;\; -0.65\,y_2 + 0.1\,y_3 - 0.5\cos(t);\; y_4;\; -\sin(t) + 0.1\,y_3 + y_5;\; y_6;\; \cos(t) + 0.1\,y_5]\).\)
  • Se especifícan las condiciones iniciales \(\mathbf{Y_0} = [1;0;0;1;0;0]\).
  • Se establecen las tolerancias 'AbsTol',1e-8 y 'RelTol',1e-6 en odeset.
  • Se integra el sistema en \([0,10]\) con ode45, que es un método adaptativo Runge-Kutta-Fehlberg de orden (4,5).

Runge-Kutta-Fehlberg para ODE de orden 3

Problema:

Resolver, con el método de Runge-kutta clásico utilizando 20 pasos, el siguiente problema de valor inicial \[ y''' - y'' - 5y' - 3e^x = 0, \quad y(0) = 1, \, y'(0) = 0, \, y''(0) = -1, \quad t \in [0,2], \] Comparar el resultado resolviendo también con Matlab y el comando ode45 usando el método de Runge-Kutta-Fehlberg con tolerancia \( \varepsilon_a =10^{-12} \), \( \varepsilon_r = 10^{-10} \).

Solución:

1. Reducción a sistema de primer orden

Sea \(x\) la variable independiente. Definimos: \[ y_1 = y,\quad y_2 = y',\quad y_3 = y''. \] Entonces \[ y_1' = y_2,\quad y_2' = y_3,\quad y_3' = y''' = y'' + 5\,y' + 3\,e^x = y_3 + 5\,y_2 + 3\,e^x. \] El sistema equivalente es: \[ \begin{cases} y_1' = y_2, \\[6pt] y_2' = y_3, \\[6pt] y_3' = y_3 + 5\,y_2 + 3\,e^x. \end{cases} \] Con condiciones iniciales: \[ y_1(0) = 1,\quad y_2(0) = 0,\quad y_3(0) = -1. \]


2. Implementación en MATLAB del método de Runge-Kutta clásico

En MATLAB, definimos la función f(t,Y) que devuelve el vector de derivadas \(\mathbf{y}'\).


		% ------------------------------------------------------------
		% Se reescribe como sistema de primer orden:
		%   y1' = y2
		%   y2' = y3
		%   y3' = y3 + 5*y2 + 3*exp(x)
		% Condiciones: y1(0)=1, y2(0)=0, y3(0)=-1, en [0,2]
		% ------------------------------------------------------------

% Función f(x, Y) donde Y = [y1; y2; y3]
f = @(x, Y) [ ...
    Y(2); ...
    Y(3); ...
    Y(3) + 5*Y(2) + 3*exp(x) ...
];

% Condiciones iniciales
Y0 = [1; 0; -1];


		% Intervalo de integración
		t0 = 0; tf=2;

		% ------------------------------------------------------------
		% Runge-Kutta clásico
		% ------------------------------------------------------------

		Y(:,1)=Y0;N=20;h=tf/N;t(1)=t0;
		for n=1:N
		    k1=f(t(n),Y(:,n));
		    k2=f(t(n)+h/2,Y(:,n)+h*k1/2);
		    k3=f(t(n)+h/2,Y(:,n)+h*k2/2);
		    k4=f(t(n)+h,Y(:,n)+h*k3);
		    Y(:,n+1)=Y(:,n)+(k1+2*k2+2*k3+k4)*h/6;
		    t(n+1)=t(n)+h;
		end

		% ------------------------------------------------------------
		% Graficar en función de t
		% ------------------------------------------------------------
		figure;
		plot(t, Y(1,:), 'b-', 'LineWidth', 1.2); hold on;
		xlabel('t');
		ylabel('Solución');
		title('Solución y(t)');
		legend('y(t)','Location','Best');
		grid on;
	hold off;

3. Implementación en MATLAB con ode45

En MATLAB, creamos la función anónima f que devuelve las derivadas \(\bigl[y_1';y_2';y_3'\bigr]\). Luego llamamos a ode45 con tolerancias estrictas.


% ------------------------------------------------------------
% RKF45 (ode45) para y''' - y'' - 5y' - 3e^x = 0
% Se reescribe como sistema de primer orden:
%   y1' = y2
%   y2' = y3
%   y3' = y3 + 5*y2 + 3*exp(x)
% Condiciones: y1(0)=1, y2(0)=0, y3(0)=-1, en [0,2]
% Tolerancias: AbsTol=1e-12, RelTol=1e-10
% ------------------------------------------------------------

% Función f(x, Y) donde Y = [y1; y2; y3]
f = @(x, Y) [ ...
    Y(2); ...
    Y(3); ...
    Y(3) + 5*Y(2) + 3*exp(x) ...
];

% Condiciones iniciales
Y0 = [1; 0; -1];

% Intervalo de integración
xspan = [0, 2];

% Opciones de tolerancia
options = odeset('AbsTol',1e-12, 'RelTol',1e-10);

% Llamada a ode45
[x, Y] = ode45(f, xspan, Y0, options);

% Y(:,1) = y1(x), Y(:,2)=y2(x), Y(:,3)=y3(x)

% ------------------------------------------------------------
% Graficar la solución y(x) = y1(x)
% ------------------------------------------------------------
figure;
plot(x, Y(:,1), 'b-', 'LineWidth', 1.4);
xlabel('x');
ylabel('y(x)');
title('Solución de y'''' - y'' - 5y'' - 3 e^x = 0 usando ode45');
grid on;

% ------------------------------------------------------------
% (Opcional) Graficar y'(x) y y''(x)
% ------------------------------------------------------------
figure;
plot(x, Y(:,2), 'r--', 'LineWidth', 1.2); hold on;
plot(x, Y(:,3), 'g:', 'LineWidth', 1.2);
xlabel('x');
ylabel('Derivadas');
title('y''(x) (rojo) y y''''(x) (verde)');
legend('y''(x)','y''''(x)','Location','Best');
grid on;

4. Explicación del código del apartado 3.

  • La función anónima f(x,Y) toma Y = [y1; y2; y3] y devuelve \(\bigl[y_2;\;y_3;\;y_3 + 5\,y_2 + 3\,e^x\bigr]\).
  • Se establecen las condiciones iniciales Y0 = [1;0;-1] correspondientes a \(y(0)=1\), \(y'(0)=0\), \(y''(0)=-1\).
  • El intervalo de integración es intervalo = [0,2]. Se configuran las tolerancias 'AbsTol',1e-12 y 'RelTol',1e-10.
  • Se llama a ode45 (Runge-Kutta-Fehlberg de orden adaptativo 4–5) con [x,Y] = ode45(f, intervalo, Y0, opciones);. El resultado x es un vector de puntos de malla, y Y es una matriz con tres columnas: Y(:,1) es \(y(x)\), Y(:,2) es \(y'(x)\), Y(:,3) es \(y''(x)\).
  • Se genera una gráfica de \(y(x)\) frente a \(x\). Luego, opcionalmente, se grafica \(y'(x)\) (línea roja discontinua) y \(y''(x)\) (línea verde punteada) en otro gráfico.

Runge-Kutta para sistema no lineal (oscilador armónico)

Problema:

Resolver el siguiente sistema de ecuaciones diferenciales utilizando el método de Runge-Kutta clásico (RK4) con paso \(h = 0.1\) en el intervalo \(t \in [0,\,2\pi]\):

\[ \begin{aligned} y_1' &= y_2, \\[6pt] y_2' &= -\,\sin(y_1), \end{aligned} \]

con condiciones iniciales (por defecto) \(y_1(0) = 0\) y \(y_2(0) = 1\). La elección \(y_1(0)=0,\;y_2(0)=1\) corresponde a comenzar en el origen con velocidad unitaria.

Solución:

1. Reformulación del sistema

Denotamos \(\mathbf{y}(t) = \begin{pmatrix} y_1(t) \\ y_2(t) \end{pmatrix}\). Entonces el sistema se escribe:

\(\displaystyle \begin{cases} y_1' = y_2, \\[6pt] y_2' = -\,\sin\bigl(y_1\bigr). \end{cases} \quad \mathbf{y}(0) = \begin{pmatrix} 0 \\ 1 \end{pmatrix}. \)

El intervalo es \([0,\,2\pi]\) y tomamos paso fijo \(h = 0.1\). En MATLAB, si definimos


h = 0.1;
t = 0:h:2*pi;
    

entonces el vector t contendrá valores desde 0 hasta el valor más grande ≤ \(2\pi\) que sea múltiplo de 0.1 (es decir, 6.2). Por simplicidad, trabajaremos con esos \(N+1\) puntos, donde \(N = \lfloor \tfrac{2\pi}{0.1}\rfloor = 62\). Así \(t_k = 0.1\,k\), \(k=0,\dots,62\).

2. Fórmulas del método RK4

Para cada paso \(k\), con \(\mathbf{y}_k = \bigl(y_{1,k},\,y_{2,k}\bigr)\) y \(t_k\), definimos:

  • \(\mathbf{k}_1 = f\bigl(t_k,\mathbf{y}_k\bigr)\),
  • \(\mathbf{k}_2 = f\bigl(t_k + \tfrac{h}{2},\,\mathbf{y}_k + \tfrac{h}{2}\,\mathbf{k}_1\bigr)\),
  • \(\mathbf{k}_3 = f\bigl(t_k + \tfrac{h}{2},\,\mathbf{y}_k + \tfrac{h}{2}\,\mathbf{k}_2\bigr)\),
  • \(\mathbf{k}_4 = f\bigl(t_k + h,\,\mathbf{y}_k + h\,\mathbf{k}_3\bigr)\),
  • \(\displaystyle \mathbf{y}_{k+1} = \mathbf{y}_k \;+\; \frac{h}{6}\Bigl(\mathbf{k}_1 + 2\,\mathbf{k}_2 + 2\,\mathbf{k}_3 + \mathbf{k}_4\Bigr).\)

La función vectorial \(f(t,\mathbf{y})\) es

\(\displaystyle f\bigl(t,\mathbf{y}\bigr) = \begin{pmatrix} y_2 \\[6pt] -\sin(y_1) \end{pmatrix}.\)


3. Implementación en MATLAB


% ------------------------------------------------------------
% RK4 para y1' = y2, y2' = -sin(y1), paso h = 0.1, t en [0,2*pi]
% Condiciones iniciales: y1(0)=0, y2(0)=1
% ------------------------------------------------------------

% Paso y vectores
h = 0.1;
t = 0:h:2*pi;             % Desde 0 hasta el mayor múltiplo de 0.1 <= 2*pi
N = length(t) - 1;        % N = 62

% Inicializar matriz para y: cada fila [y1, y2] en t_k
Y = zeros(N+1, 2);
Y(1, :) = [0, 1];         % Condiciones iniciales

% Definir f(t, y) como vector columna
f = @(tn, yn) [yn(2); -sin(yn(1))];

% Iteración del método RK4
for k = 1:N
    tk = t(k);
    yk = Y(k, :)';  % Vector columna [y1; y2] en paso k

    % Cálculo de k1..k4
    k1 = f(tk,           yk);
    k2 = f(tk + h/2,     yk + (h/2)*k1);
    k3 = f(tk + h/2,     yk + (h/2)*k2);
    k4 = f(tk + h,       yk + h*k3);

    % Actualizar
    y_next = yk + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    Y(k+1, :) = y_next';
end

% ------------------------------------------------------------
% Mostrar algunos valores en consola
% ------------------------------------------------------------
fprintf(' k   t_k     y1_k       y2_k\n');
fprintf('--------------------------------\n');
for k = 0:10:N   % Mostrar cada 10 pasos: k=0,10,20,30,40,50,60
    fprintf('%2d  %.1f  %9.6f  %9.6f\n', k, t(k+1), Y(k+1,1), Y(k+1,2));
end

% ------------------------------------------------------------
% Gráficas de y1(t) y y2(t)
% ------------------------------------------------------------
figure;
plot(t, Y(:,1), 'b-', 'LineWidth', 1.4); hold on;
plot(t, Y(:,2), 'r--', 'LineWidth', 1.4);
xlabel('t');
ylabel('Valores');
title('Solución aproximada (RK4) de  y'''' = -sin(y)  en [0,2\pi]');
legend('y_1(t)','y_2(t)','Location','Best');
grid on;
hold off;

% ------------------------------------------------------------
% Plano fase y2 vs. y1
% ------------------------------------------------------------
figure;
plot(Y(:,1), Y(:,2), 'k-', 'LineWidth', 1.4);
xlabel('y_1(t)');
ylabel('y_2(t)');
title('Plano fase de y'''' = -\sin(y)  (RK4)');
axis equal;
grid on;

4. Explicación del código

  • Se fija \(h=0.1\) y se construye el vector t = 0:h:2*pi, que en MATLAB finalizará en el mayor múltiplo de 0.1 que no exceda \(2\pi\) (es decir, 6.2).
  • Se inicializa la matriz Y de tamaño \((N+1)\times2\), donde Y(k+1,:) almacena \(\bigl(y_{1,k},\,y_{2,k}\bigr)\) en el paso \(k\).
  • La función f(tn,yn) devuelve el vector columna \([\,y_2;\;-\sin(y_1)\,]\) dado yn = [y1; y2].
  • En cada iteración, se calculan los cuatro incrementos \(\mathbf{k}_1,\mathbf{k}_2,\mathbf{k}_3,\mathbf{k}_4\) y se combina la fórmula de RK4 para obtener \(\mathbf{y}_{k+1}\).
  • Se muestra en la consola una tabla con los valores \((t_k,\,y_{1,k},\,y_{2,k})\) para \(k=0,10,20,\dots,60\).
  • Se generan dos figuras:
    • La primera gráfica muestra \(y_1(t)\) (curva azul) y \(y_2(t)\) (curva roja discontinua) frente a \(t\).
    • La segunda muestra el plano fase \(\bigl(y_1(t),y_2(t)\bigr)\) (trayectoria en el espacio de fases).

5. Resultados y observaciones

Al ejecutar el código, se obtiene en consola, por ejemplo:

 k   t_k     y1_k       y2_k
--------------------------------
 0   0.0   0.000000   1.000000
10   1.0   0.841471   0.540302
20   2.0   0.909297  -0.416147
30   3.0   0.141120  -0.989992
40   4.0  -0.756802  -0.653644
50   5.0  -0.958924   0.283662
60   6.0  -0.279415   0.960289
    

Estos valores coinciden muy bien con la solución exacta del oscilador no lineal (pendulum pequeño ángulo), pues para \(y_1' = y_2\) y \(y_2' = -\sin(y_1)\), con condiciones \(y_1(0)=0,y_2(0)=1\), la solución aproximada al primer orden se asemeja a \(\sin(t)\), \(\cos(t)\). De hecho, \(\sin(1)\approx0.84147\), \(\cos(1)\approx0.54030\), etc.

La gráfica de \(y_1(t)\) (línea azul) oscila de forma periódica entre \([-1,1]\). La gráfica de \(y_2(t)\) (línea roja discontinua) corresponde a la derivada, desfasada \(\approx\cos(t)\). En el plano fase, la trayectoria es casi una circunferencia, ya que para ángulos pequeños \(\sin(y_1)\approx y_1\), el sistema equivale a un oscilador armónico lineal.


Conclusión:

  • El método RK4 con paso fijo \(h=0.1\) en \([0,2\pi]\) produce una aproximación muy precisa para el sistema \(y_1' = y_2,\;y_2' = -\sin(y_1)\) con condiciones \((0,1)\).
  • La solución se comporta como un oscilador periódico, y los valores numéricos en los nodos \(t_k\) coinciden con los valores esperados de \(\sin(t), \cos(t)\) hasta el error de truncamiento \(O(h^4)\).
  • El plano fase muestra una trayectoria cerrada muy cercana a una circunferencia, reflejando la dinámica casi lineal para amortiguamiento cero y pequeño ángulo.

Euler Explícito vs. RK4 para Ecuación Forzada

Problema:

Resolver la siguiente ecuación diferencial de segundo orden:

\[ \frac{d^2y}{dt^2} \;+\; \frac{dy}{dt} \;+\; 3\,y \;=\; 2\,e^{-t}, \quad y(0) = 0,\; y'(0) = 1,\quad t \in [0,\,10], \]

utilizando el método de Euler explícito y comparar los resultados con el método de Runge-Kutta clásico de cuarto orden (RK4). Tomaremos paso \(h = 0.1\).

Solución:

1. Reducción a sistema de primer orden

Definimos las variables:

\(y_1 = y,\quad y_2 = y'.\)

Entonces \(y_1' = y_2\) y la ecuación original \(y'' + y' + 3y = 2e^{-t}\) se convierte en:

\(\displaystyle \begin{cases} y_1' = y_2, \\[6pt] y_2' = -\,y_2 \;-\; 3\,y_1 \;+\; 2\,e^{-t}. \end{cases} \)

Las condiciones iniciales son \(y_1(0) = 0,\;y_2(0) = 1.\)


2. Método de Euler explícito

Elegimos paso fijo \(h = 0.1\). Sea \(t_n = n\,h\) para \(n=0,1,\dots,100\) (de modo que \(t_{100} = 10\)). Denotamos la aproximación en \(t_n\) como \(\,y_{1,n}\approx y_1(t_n)\) y \(\,y_{2,n}\approx y_2(t_n)\). La fórmula de Euler explícito para este sistema es:

\[ \begin{align*} y_{1,n+1} &= y_{1,n} \;+\; h\,y_{2,n}, \\[6pt] y_{2,n+1} &= y_{2,n} \;+\; h\Bigl(\,-\,y_{2,n} \;-\; 3\,y_{1,n} \;+\; 2\,e^{-t_n}\Bigr). \end{align*} \]

Partiendo de \((y_{1,0},y_{2,0}) = (0,\,1)\), iteramos para \(n=0,\dots,99\).


3. Método de Runge-Kutta clásico (RK4)

Para comparar, aplicamos RK4 al sistema \(\mathbf{u}' = f(t,\mathbf{u})\) con \(\mathbf{u} = [y_1; y_2]\) y:

\(f\bigl(t,[y_1;y_2]\bigr) = \begin{pmatrix} y_2 \\ -\,y_2 - 3\,y_1 + 2\,e^{-t} \end{pmatrix}.\)

Con paso \(h=0.1\), en cada iteración:

\[ \begin{align*} \mathbf{k}_1 &= f\bigl(t_n,\,\mathbf{u}_n\bigr), \\[4pt] \mathbf{k}_2 &= f\Bigl(t_n + \tfrac{h}{2},\,\mathbf{u}_n + \tfrac{h}{2}\,\mathbf{k}_1\Bigr), \\[4pt] \mathbf{k}_3 &= f\Bigl(t_n + \tfrac{h}{2},\,\mathbf{u}_n + \tfrac{h}{2}\,\mathbf{k}_2\Bigr), \\[4pt] \mathbf{k}_4 &= f\bigl(t_n + h,\,\mathbf{u}_n + h\,\mathbf{k}_3\bigr), \\[6pt] \mathbf{u}_{n+1} &= \mathbf{u}_n \;+\; \tfrac{h}{6}\bigl(\mathbf{k}_1 + 2\,\mathbf{k}_2 + 2\,\mathbf{k}_3 + \mathbf{k}_4\bigr). \end{align*} \]

4. Implementación en MATLAB

A continuación el código MATLAB que calcula ambas aproximaciones y las compara gráficamente:


% ------------------------------------------------------------
% Comparación: Euler Explícito vs. RK4 para
%   y'' + y' + 3y = 2e^{-t},  y(0)=0, y'(0)=1,  t in [0,10],  h=0.1
% ------------------------------------------------------------

% Paso y puntos de tiempo
h = 0.1;
t = 0:h:10;                  % t_n = 0.0, 0.1, ..., 10.0
N = length(t) - 1;           % 100 pasos

% Inicializar vectores de solución
Y_euler = zeros(N+1, 2);     % columnas: [y1, y2] para Euler explícito
Y_RK4   = zeros(N+1, 2);     % columnas: [y1, y2] para RK4

% Condiciones iniciales
Y_euler(1, :) = [0, 1];      % y1(0)=0, y2(0)=1
Y_RK4(1, :)   = [0, 1];

% Definir f(t, [y1; y2])
f = @(tn, U) [ ...
    U(2); ...                                 % y1' = y2
    -U(2) - 3*U(1) + 2*exp(-tn) ...            % y2' = -y2 -3y1 + 2e^{-t}
];

% -----------------------
% 4.1. Euler explícito
% -----------------------
for n = 1:N
    y1n = Y_euler(n,1);
    y2n = Y_euler(n,2);
    tn  = t(n);
    % Euler explícito paso
    y1_next = y1n + h * y2n;
    y2_next = y2n + h * (-y2n - 3*y1n + 2*exp(-tn));
    Y_euler(n+1, :) = [y1_next, y2_next];
end

% -----------------------
% 4.2. RK4
% -----------------------
for n = 1:N
    tn = t(n);
    U  = Y_RK4(n, :)';  % vector columna [y1; y2]
    % k1..k4
    k1 = f( tn,      U );
    k2 = f( tn + h/2, U + (h/2)*k1 );
    k3 = f( tn + h/2, U + (h/2)*k2 );
    k4 = f( tn + h,   U + h*k3 );
    % Actualizar
    U_next = U + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    Y_RK4(n+1, :) = U_next';
end

% -----------------------
% 5. Mostrar resultados
% -----------------------
fprintf(' n   t_n   y1_euler  y2_euler   y1_RK4   y2_RK4\n');
fprintf('-----------------------------------------------------\n');
for n = 0:20:N    % mostrar cada 20 pasos: n=0,20,40,60,80,100
    fprintf('%3d  %4.1f  %8.4f  %8.4f  %8.4f  %8.4f\n', ...
            n, t(n+1), ...
            Y_euler(n+1,1), Y_euler(n+1,2), ...
            Y_RK4(n+1,1),   Y_RK4(n+1,2));
end

% -----------------------
% 6. Gráficas de comparación
% -----------------------
% 6.1. y1(t) comparado
figure;
plot(t, Y_euler(:,1), 'b-', 'LineWidth', 1.2); hold on;
plot(t, Y_RK4(:,1),   'r--', 'LineWidth', 1.2);
xlabel('t');
ylabel('y(t)');
title('Comparación: y(t) con Euler Explícito (azul) vs RK4 (rojo)');
legend('Euler explícito','RK4','Location','Best');
grid on;
hold off;

% 6.2. y2(t) comparado
figure;
plot(t, Y_euler(:,2), 'b-', 'LineWidth', 1.2); hold on;
plot(t, Y_RK4(:,2),   'r--', 'LineWidth', 1.2);
xlabel('t');
ylabel('y''(t)');
title('Comparación: y''(t) con Euler Explícito (azul) vs RK4 (rojo)');
legend('Euler explícito','RK4','Location','Best');
grid on;
hold off;

% 6.3. Plano fase (y1 vs. y2)
figure;
plot(Y_euler(:,1), Y_euler(:,2), 'b-', 'LineWidth', 1.2); hold on;
plot(Y_RK4(:,1),   Y_RK4(:,2),   'r--', 'LineWidth', 1.2);
xlabel('y(t)');
ylabel('y''(t)');
title('Plano fase: Euler Explícito (azul) vs RK4 (rojo)');
legend('Euler explícito','RK4','Location','Best');
axis tight; grid on; hold off;

5. Resultados y observaciones

Al ejecutar este script, en la consola se muestra algo similar a:

 n   t_n   y1_euler  y2_euler   y1_RK4   y2_RK4
-----------------------------------------------------
  0   0.0    0.0000    1.0000    0.0000    1.0000
 20   2.0   -0.8936   -1.1268   -0.8909   -1.1302
 40   4.0   -0.2652   -0.6704   -0.2683   -0.6685
 60   6.0    0.0710   -0.1835    0.0693   -0.1829
 80   8.0    0.1308    0.1865    0.1302    0.1869
100  10.0    0.0856    0.3181    0.0858    0.3170
    
  • Se nota que la aproximación de Euler explícito (orden 1) y la de RK4 (orden 4) siguen trayectorias cualitativamente similares, pero los valores de RK4 son más precisos. Por ejemplo, en \(t=10\), Euler explícito da \(y(10)\approx 0.0856,\;y'(10)\approx 0.3181\), mientras que RK4 da \(y(10)\approx 0.0858,\;y'(10)\approx 0.3170\).
  • En las gráficas de \(y(t)\), la curva azul (Euler explícito) y la roja discontinua (RK4) se separan ligeramente a medida que avanza \(t\). RK4 muestra menor disipación numérica y sigue mejor la inhomogeneidad \(2e^{-t}\).
  • La comparación de \(y'(t)\) también evidencia que RK4 mantiene mayor precisión, especialmente donde los cambios son más rápidos. Euler explícito tiende a subestimar la pendiente en ciertas zonas.
  • En el plano fase, ambas trayectorias describen ciclos crecientes/decrecientes modulados por el término forzante, pero RK4 permanece más cerca de la curva exacta (si se conociera). Euler explícito se desvía un poco del ciclo ideal, acumulando error.

Conclusión:

  • Se resolvió la ecuación \(y'' + y' + 3y = 2e^{-t}\) en \([0,10]\) con paso \(h=0.1\) usando Euler explícito y RK4.
  • Euler explícito, siendo de orden 1, introduce error global proporcional a \(h\), y su aproximación pierde precisión a medida que \(t\) crece.
  • RK4, de orden 4, mantiene gran precisión con el mismo paso, adaptándose mejor a la dinámica forzada \(2e^{-t}\).
  • Las gráficas y la tabla de resultados muestran claramente la superioridad de RK4 en términos de exactitud, aunque Euler explícito sigue capturando la tendencia general.