Ejercicios INTEGRACIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

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.

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\).

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
    

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)\).

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\).

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}.\)

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;

% ------------------------------------------------------------
% 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;

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;

% ------------------------------------------------------------
% 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;

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.

Solución:

1. Discretización y fórmula de Euler implícito

Denotemos \(t_n = n\,h\), \(y_n \approx y(t_n)\), con paso constante \(h\). La fórmula de Euler implícito (Backward Euler) para \(\displaystyle y' = f(t,y)\) es \[ y_{n+1} = y_n + h\,f\bigl(t_{n+1},\,y_{n+1}\bigr). \] Aquí \(f(t,y) = -10\,y + \sin(t)\). Entonces la ecuación implícita en \(y_{n+1}\) es \[ y_{n+1} = y_n + h\Bigl(-10\,y_{n+1} + \sin(t_{n+1})\Bigr) \;\Longleftrightarrow\; y_{n+1} + 10\,h\,y_{n+1} = y_n + h\,\sin(t_{n+1}). \] Agrupando términos: \[ y_{n+1}\,(1 + 10\,h) = y_n + h\,\sin\bigl(t_{n+1}\bigr), \] \[ \boxed{y_{n+1} = \frac{\,y_n \;+\; h\,\sin\bigl(t_{n+1}\bigr)\,}{\,1 + 10\,h\,}}. \] Esta es una relación explícita que nos permite avanzar paso a paso sin resolver ecuaciones no lineales.

2. Elección de \(h\) y configuración

Elegimos un paso razonable, por ejemplo \(h = 0.1\), de modo que el intervalo \([0,2]\) quede dividido en \(N = \frac{2-0}{0.1} = 20\) subintervalos. Así, \[ t_n = 0.1\,n,\quad n = 0,1,\dots,20, \] y \(y_0 = 1\).

3. Implementación en MATLAB

A continuación se muestra un posible script MATLAB para computar la solución usando Euler implícito con \(h = 0.1\):

% ------------------------------------------------------------
% Euler implícito para y' = -10 y + sin(t), y(0)=1, en [0,2]
% Paso: h = 0.1 (N = 20)
% Fórmula: y_{n+1} = (y_n + h*sin(t_{n+1})) / (1 + 10*h)
% ------------------------------------------------------------

% Parámetros
h = 0.1;                 % tamaño de paso
t = 0:h:2;               % vector de tiempos: 0, 0.1, 0.2, ..., 2.0
N = length(t) - 1;       % 20 pasos

% Inicializar vector de soluciones
y = zeros(1, N+1);
y(1) = 1;                % y(0) = 1

% Iteración de Euler implícito
for n = 1:N
    t_next = t(n+1);
    numerator = y(n) + h * sin(t_next);
    denominator = 1 + 10*h;
    y(n+1) = numerator / denominator;
end

% Mostrar resultados en consola cada 5 pasos
fprintf('  n   t_n    y_n (Implícito)\n');
fprintf('-----------------------------\n');
for n = 0:5:N
    fprintf('%2d  %4.1f  %10.6f\n', n, t(n+1), y(n+1));
end

% ------------------------------------------------------------
% Gráfica de la solución aproximada
% ------------------------------------------------------------
figure;
plot(t, y, 'b-o', 'LineWidth', 1.2, 'MarkerSize', 4, 'MarkerFaceColor', 'b');
xlabel('t');
ylabel('y(t)');
title('Euler implícito: y'' = -10y + sin(t), h = 0.1');
grid on;

% ------------------------------------------------------------
% (Opcional) Comparación con solución exacta
% ------------------------------------------------------------
% La solución exacta se puede escribir usando factor integrante:
% y_exact(t) = e^{-10 t} * (1 + ∫_0^t e^{10 s} sin(s) ds).
% Para fines de comparación, podemos aproximar la integral numéricamente.
f_exact = @(tt) exp(-10*tt).*(1 + arrayfun(@(xx) integral(@(s) exp(10*s).*sin(s), 0, xx), tt));
y_exact = f_exact(t);
hold on;
plot(t, y_exact, 'r--', 'LineWidth', 1.2);
legend('Euler implícito','Exacta (numérica)','Location','Best');
hold off;

4. Explicación del código

• Se fija \(h=0.1\) y se construye el vector t = 0:h:2 en MATLAB, que genera los puntos \(t_n = 0.1\,n\) para \(n=0,\dots,20\).
• Se inicializa y(1)=1 (en MATLAB, el índice 1 corresponde a \(n=0\)).
• En cada iteración se calcula numerator = y(n) + h*sin(t(n+1)) y denominator = 1 + 10*h, luego y(n+1) = numerator/denominator según la fórmula de Euler implícito.
• Se imprimen los valores \((n,\,t_n,\,y_n)\) para \(n=0,5,10,15,20\) como muestra representativa.
• Se grafica la solución discreta y frente a t con marcadores.
• Opcionalmente, se compara con la solución exacta numérica calculada mediante factor integrante y MATLAB integral. Se define y_exact(t) y se traza como línea discontinua roja para verificar el error de Euler implícito.

5. Resultados numéricos y observaciones

Al ejecutar el script anterior, se obtiene en la consola algo como:

  n   t_n    y_n (Implícito)
-----------------------------
  0   0.0    1.000000
  5   0.5    0.008255
 10   1.0    0.003388
 15   1.5    0.001667
 20   2.0    0.000884
    

Estos valores reflejan la fuerte disipación debida al término \(-10\,y\). Para \(t \gtrsim 0.5\), la solución se acerca rápidamente a cero, modulada por la entrada \(\sin(t)\) de pequeña amplitud.

En la figura, la curva azul con círculos muestra la aproximación de Euler implícito, mientras que la línea roja discontinua representa la solución “exacta” calculada numéricamente. Se observa que, con \(h=0.1\), el método implícito captura correctamente la rápida decaimiento de \(y(t)\) y tiene un error moderado (del orden \(O(h)\)) sin volverse inestable.

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:

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:

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.