import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def euler_method(f, y0, t, *args):
"""
Forward Euler method for solving ODEs.
dy/dt = f(y, t, *args)
y(t0) = y0
"""
y = np.zeros((len(t), len(y0) if hasattr(y0, '__len__') else 1))
y[0] = y0
for i in range(len(t) - 1):
dt = t[i+1] - t[i]
y[i+1] = y[i] + dt * np.array(f(y[i], t[i], *args))
return y
def rk4_method(f, y0, t, *args):
"""
4th order Runge-Kutta method for solving ODEs.
dy/dt = f(y, t, *args)
y(t0) = y0
"""
y = np.zeros((len(t), len(y0) if hasattr(y0, '__len__') else 1))
y[0] = y0
for i in range(len(t) - 1):
dt = t[i+1] - t[i]
k1 = np.array(f(y[i], t[i], *args))
k2 = np.array(f(y[i] + dt*k1/2, t[i] + dt/2, *args))
k3 = np.array(f(y[i] + dt*k2/2, t[i] + dt/2, *args))
k4 = np.array(f(y[i] + dt*k3, t[i] + dt, *args))
y[i+1] = y[i] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)
return y
# Test with exponential decay
def decay_ode(y, t, k):
return -k * y
k = 0.5
y0 = [10.0]
t_coarse = np.linspace(0, 10, 21) # 20 steps
t_fine = np.linspace(0, 10, 201) # 200 steps
# Analytical solution
y_exact = y0[0] * np.exp(-k * t_fine)
# Solve with different methods
y_euler_coarse = euler_method(decay_ode, y0, t_coarse, k)
y_euler_fine = euler_method(decay_ode, y0, t_fine, k)
y_rk4_coarse = rk4_method(decay_ode, y0, t_coarse, k)
y_rk4_fine = rk4_method(decay_ode, y0, t_fine, k)
y_odeint = odeint(decay_ode, y0, t_fine, args=(k,))
# Plot comparison
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Left plot: Visual comparison
ax1.plot(t_fine, y_exact, 'k-', linewidth=2, label='Exact', alpha=0.7)
ax1.plot(t_coarse, y_euler_coarse, 'r.-', linewidth=1.5, markersize=8, label='Euler (20 steps)')
ax1.plot(t_coarse, y_rk4_coarse, 'b.-', linewidth=1.5, markersize=8, label='RK4 (20 steps)')
ax1.plot(t_fine, y_odeint, 'g--', linewidth=2, label='LSODA (adaptive)', alpha=0.7)
ax1.set_xlabel('Time')
ax1.set_ylabel('y(t)')
ax1.set_title('Method Comparison: Exponential Decay')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Right plot: Error analysis
error_euler_coarse = np.abs(y_euler_coarse.flatten() - y0[0] * np.exp(-k * t_coarse))
error_euler_fine = np.abs(y_euler_fine.flatten() - y_exact)
error_rk4_coarse = np.abs(y_rk4_coarse.flatten() - y0[0] * np.exp(-k * t_coarse))
error_rk4_fine = np.abs(y_rk4_fine.flatten() - y_exact)
error_odeint = np.abs(y_odeint.flatten() - y_exact)
ax2.semilogy(t_coarse, error_euler_coarse, 'r.-', linewidth=1.5, markersize=8, label='Euler (20 steps)')
ax2.semilogy(t_fine, error_euler_fine, 'r--', linewidth=1, alpha=0.5, label='Euler (200 steps)')
ax2.semilogy(t_coarse, error_rk4_coarse, 'b.-', linewidth=1.5, markersize=8, label='RK4 (20 steps)')
ax2.semilogy(t_fine, error_rk4_fine, 'b--', linewidth=1, alpha=0.5, label='RK4 (200 steps)')
ax2.semilogy(t_fine, error_odeint, 'g-', linewidth=2, label='LSODA', alpha=0.7)
ax2.set_xlabel('Time')
ax2.set_ylabel('Absolute Error')
ax2.set_title('Error Analysis (Log Scale)')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Print final errors
print("Final Errors at t=10:")
print(f" Euler (20 steps): {error_euler_coarse[-1]:.6e}")
print(f" Euler (200 steps): {error_euler_fine[-1]:.6e}")
print(f" RK4 (20 steps): {error_rk4_coarse[-1]:.6e}")
print(f" RK4 (200 steps): {error_rk4_fine[-1]:.6e}")
print(f" LSODA (adaptive): {error_odeint[-1]:.6e}")
