#### Introduction
Numerical methods for solving ordinary differential equations (ODEs) provide practical approaches to approximate solutions when analytical expressions are unavailable or difficult to obtain. These methods discretize the time domain and iteratively compute approximate values of the solution to the initial value problem
$
\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0,
$
by stepping forward in time. The simplest among these is the [[Explicit Euler Method]], which uses the slope at the current point to estimate the solution at the next time step. While straightforward and computationally inexpensive, this method is only first-order accurate and can suffer from stability issues for stiff problems.
To improve accuracy, the [[Runge-Kutta Midpoint Method]] evaluates the slope at the midpoint of the interval, effectively incorporating more information about the solution's behavior within each step. This second-order Runge-Kutta method achieves better accuracy and stability compared to Euler's method by computing an intermediate slope and using it to update the solution. For even higher accuracy, the [[Fourth-Order Runge-Kutta]] method takes four slope evaluations per step and combines them in a weighted average, resulting in a fourth-order accurate method widely used in practice due to its balance of precision and computational cost.
The accompanying code examples and comparative analyses in these pages demonstrate the implementation and performance of these methods on test problems such as the exponential decay ODE $dy/dt = -2y$ with initial condition $y(0) = 1$. The [[Numerical ODE Code Comparisons]] page provides visualizations and error analyses that highlight the convergence properties and relative accuracies of these methods, illustrating the trade-offs between simplicity, computational effort, and solution quality.
#### Code Snippets
``` MATLAB
function [t,y] = myEuler(f,range,iv,h)
t = range(1):h:range(2);
y = zeros(length(t),1);
y(1) = iv;
for i = 1:length(t) - 1
y(i+1) = y(i) + h*f(t(i),y(i));
end
```
``` Python
def euler(f, range_vals, iv, h):
"""
Euler's method for solving ODEs
Parameters:
f: function that takes (t, y) and returns dy/dt
range_vals: tuple (t_start, t_end) defining integration range
iv: initial value y(t_start)
h: step size
Returns:
t: array of time points
y: array of solution values
"""
t = np.arange(range_vals[0], range_vals[1] + h, h)
y = np.zeros(len(t))
y[0] = iv
for i in range(len(t) - 1):
y[i+1] = y[i] + h * f(t[i], y[i])
return t, y
```
#### Example Driver
Note, these drivers assume that you've also created methods for [[Explicit Euler Method]], [[Runge-Kutta Midpoint Method]], and [[Fourth-Order Runge-Kutta]].
``` MATLAB
% Example usage of numerical integration methods
% Solve dy/dt = -2*y with y(0) = 1
% Analytical solution: y(t) = exp(-2*t)
clear; clc;
% Define the ODE function
f = @(t, y) -2 * y;
% Parameters
range_vals = [0, 2]; % integrate from t=0 to t=2
initial_value = 1.0;
step_size = 0.1;
% Solve using all three methods
[t_euler, y_euler] = myEuler(f, range_vals, initial_value, step_size);
[t_rk4, y_rk4] = rk4(f, range_vals, initial_value, step_size);
[t_mid, y_mid] = rkMid2(f, range_vals, initial_value, step_size);
% Analytical solution for comparison
t_analytical = t_euler; % same time points
y_analytical = exp(-2 * t_analytical);
% Display results
fprintf('At t = %.1f:\n', t_euler(end));
fprintf('Euler: %.6f\n', y_euler(end));
fprintf('RK4: %.6f\n', y_rk4(end));
fprintf('Midpoint: %.6f\n', y_mid(end));
fprintf('Analytical: %.6f\n', y_analytical(end));
% Plot comparison
figure;
plot(t_analytical, y_analytical, 'k-', 'LineWidth', 2, 'DisplayName', 'Analytical');
hold on;
plot(t_euler, y_euler, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Euler');
plot(t_rk4, y_rk4, 'b:', 'LineWidth', 1.5, 'DisplayName', 'RK4');
plot(t_mid, y_mid, 'g-.', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
hold off;
xlabel('Time (t)');
ylabel('y(t)');
title('Comparison of Numerical Integration Methods');
legend('Location', 'best');
grid on;
% Calculate and display errors
error_euler = abs(y_euler - y_analytical);
error_rk4 = abs(y_rk4 - y_analytical);
error_mid = abs(y_mid - y_analytical);
fprintf('\nMaximum absolute errors:\n');
fprintf('Euler: %.2e\n', max(error_euler));
fprintf('RK4: %.2e\n', max(error_rk4));
fprintf('Midpoint: %.2e\n', max(error_mid));
% Plot errors
figure;
semilogy(t_euler, error_euler, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Euler');
hold on;
semilogy(t_rk4, error_rk4, 'b:', 'LineWidth', 1.5, 'DisplayName', 'RK4');
semilogy(t_mid, error_mid, 'g-.', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
hold off;
xlabel('Time (t)');
ylabel('Absolute Error');
title('Error Comparison of Numerical Integration Methods');
% legend('Location', 'best');
grid on;
% Test with different step sizes to show convergence
step_sizes = [0.2, 0.1, 0.05, 0.025, 0.0125];
final_errors = zeros(3, length(step_sizes));
for i = 1:length(step_sizes)
h = step_sizes(i);
[t_e, y_e] = myEuler(f, range_vals, initial_value, h);
[t_r, y_r] = rk4(f, range_vals, initial_value, h);
[t_m, y_m] = rkMid2(f, range_vals, initial_value, h);
% Calculate analytical solution at final time
y_true = exp(-2 * range_vals(2));
final_errors(1, i) = abs(y_e(end) - y_true); % Euler
final_errors(2, i) = abs(y_r(end) - y_true); % RK4
final_errors(3, i) = abs(y_m(end) - y_true); % Midpoint
end
% Plot convergence
figure;
loglog(step_sizes, final_errors(1, :), 'r--o', 'LineWidth', 1.5, 'DisplayName', 'Euler');
hold on;
loglog(step_sizes, final_errors(2, :), 'b:s', 'LineWidth', 1.5, 'DisplayName', 'RK4');
loglog(step_sizes, final_errors(3, :), 'g-.^', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
% Add theoretical convergence lines
h1 = loglog(step_sizes, step_sizes, 'k--', 'DisplayName', 'O(h)');
h2 = loglog(step_sizes, step_sizes.^2, 'k:', 'DisplayName', 'O(h²)');
h3 = loglog(step_sizes, step_sizes.^4, 'k-.', 'DisplayName', 'O(h⁴)');
% Make theoretical lines semi-transparent by adjusting color
set(h1, 'Color', [0.5 0.5 0.5]);
set(h2, 'Color', [0.5 0.5 0.5]);
set(h3, 'Color', [0.5 0.5 0.5]);
hold off;
xlabel('Step Size (h)');
ylabel('Absolute Error at t = 2');
title('Convergence Analysis');
legend('Location', 'best');
grid on;
fprintf('\nConvergence analysis complete. Check the plots!\n');
```
``` MATLAB
% Example usage of numerical integration methods
% Solve dy/dt = -2*y with y(0) = 1
% Analytical solution: y(t) = exp(-2*t)
clear; clc;
% Define the ODE function
f = @(t, y) -2 * y;
% Parameters
range_vals = [0, 2]; % integrate from t=0 to t=2
initial_value = 1.0;
step_size = 0.1;
% Solve using all three methods
[t_euler, y_euler] = myEuler(f, range_vals, initial_value, step_size);
[t_rk4, y_rk4] = rk4(f, range_vals, initial_value, step_size);
[t_mid, y_mid] = rkMid2(f, range_vals, initial_value, step_size);
% Analytical solution for comparison
t_analytical = t_euler; % same time points
y_analytical = exp(-2 * t_analytical);
% Display results
fprintf('At t = %.1f:\n', t_euler(end));
fprintf('Euler: %.6f\n', y_euler(end));
fprintf('RK4: %.6f\n', y_rk4(end));
fprintf('Midpoint: %.6f\n', y_mid(end));
fprintf('Analytical: %.6f\n', y_analytical(end));
% Plot comparison
figure;
plot(t_analytical, y_analytical, 'k-', 'LineWidth', 2, 'DisplayName', 'Analytical');
hold on;
plot(t_euler, y_euler, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Euler');
plot(t_rk4, y_rk4, 'b:', 'LineWidth', 1.5, 'DisplayName', 'RK4');
plot(t_mid, y_mid, 'g-.', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
hold off;
xlabel('Time (t)');
ylabel('y(t)');
title('Comparison of Numerical Integration Methods');
legend('show', 'Location', 'best');
grid on;
% Calculate and display errors
error_euler = abs(y_euler - y_analytical);
error_rk4 = abs(y_rk4 - y_analytical);
error_mid = abs(y_mid - y_analytical);
fprintf('\nMaximum absolute errors:\n');
fprintf('Euler: %.2e\n', max(error_euler));
fprintf('RK4: %.2e\n', max(error_rk4));
fprintf('Midpoint: %.2e\n', max(error_mid));
% Plot errors
figure;
semilogy(t_euler, error_euler, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Euler');
hold on;
semilogy(t_rk4, error_rk4, 'b:', 'LineWidth', 1.5, 'DisplayName', 'RK4');
semilogy(t_mid, error_mid, 'g-.', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
hold off;
xlabel('Time (t)');
ylabel('Absolute Error');
title('Error Comparison of Numerical Integration Methods');
legend('show', 'Location', 'best');
grid on;
% Test with different step sizes to show convergence
step_sizes = [0.2, 0.1, 0.05, 0.025, 0.0125];
final_errors = zeros(3, length(step_sizes));
for i = 1:length(step_sizes)
h = step_sizes(i);
[t_e, y_e] = myEuler(f, range_vals, initial_value, h);
[t_r, y_r] = rk4(f, range_vals, initial_value, h);
[t_m, y_m] = rkMid2(f, range_vals, initial_value, h);
% Calculate analytical solution at final time
y_true = exp(-2 * range_vals(2));
final_errors(1, i) = abs(y_e(end) - y_true); % Euler
final_errors(2, i) = abs(y_r(end) - y_true); % RK4
final_errors(3, i) = abs(y_m(end) - y_true); % Midpoint
end
% Plot convergence
figure;
loglog(step_sizes, final_errors(1, :), 'r--o', 'LineWidth', 1.5, 'DisplayName', 'Euler');
hold on;
loglog(step_sizes, final_errors(2, :), 'b:s', 'LineWidth', 1.5, 'DisplayName', 'RK4');
loglog(step_sizes, final_errors(3, :), 'g-.^', 'LineWidth', 1.5, 'DisplayName', 'Midpoint');
% Add theoretical convergence lines
loglog(step_sizes, step_sizes, 'k--', 'Alpha', 0.5, 'DisplayName', 'O(h)');
loglog(step_sizes, step_sizes.^2, 'k:', 'Alpha', 0.5, 'DisplayName', 'O(h²)');
loglog(step_sizes, step_sizes.^4, 'k-.', 'Alpha', 0.5, 'DisplayName', 'O(h⁴)');
hold off;
xlabel('Step Size (h)');
ylabel('Absolute Error at t = 2');
title('Convergence Analysis');
legend('show', 'Location', 'best');
grid on;
fprintf('\nConvergence analysis complete. Check the plots!\n');
```
``` Python - Example Usage
# Example usage:
if __name__ == "__main__":
# Example: solve dy/dt = -2*y with y(0) = 1
# Analytical solution: y(t) = exp(-2*t)
def example_ode(t, y):
return -2 * y
# Parameters
range_vals = (0, 2) # integrate from t=0 to t=2
initial_value = 1.0
step_size = 0.1
# Solve using all three methods
t_euler, y_euler = euler(example_ode, range_vals, initial_value, step_size)
t_rk4, y_rk4 = rk4(example_ode, range_vals, initial_value, step_size)
t_mid, y_mid = rk_mid2(example_ode, range_vals, initial_value, step_size)
# Analytical solution for comparison
t_analytical = t_euler # same time points
y_analytical = np.exp(-2 * t_analytical)
print(f"At t = {t_euler[-1]:.1f}:")
print(f"Euler: {y_euler[-1]:.6f}")
print(f"RK4: {y_rk4[-1]:.6f}")
print(f"Midpoint: {y_mid[-1]:.6f}")
print(f"Analytical: {y_analytical[-1]:.6f}")
```