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-8El 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)
recibeX = [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 enY
(cada columna deY
corresponde a un instantet(k)
). - Se grafica
Y(1,:)
frente at
para mostrar la posición \(x(t)\) amortiguada, y en una figura aparteY(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 matrizY
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
enodeset
. - 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 comandoode45
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)
tomaY = [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 resultadox
es un vector de puntos de malla, yY
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\), dondeY(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)\,]\) dadoyn = [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.