Belajar Math - Differential Equations & Dynamic Systems

import math
 
# Contoh 1: dy/dt = 2t → y = t² + C
def f1(t, y):
    return 2 * t
 
# Solusi analitis
def y1_analytical(t):
    return t**2  # C=0 untuk y(0)=0
 
# Contoh 2: dy/dt = -2y → y = y₀ × e^(-2t)
def f2(t, y):
    return -2 * y
 
def y2_analytical(t, y0=1):
    return y0 * math.exp(-2 * t)
 
print("Contoh ODE:")
print(f"  dy/dt = 2t → y(t) = t²")
print(f"  dy/dt = -2y → y(t) = y₀ × e^(-2t)")
print(f"\ny(3) untuk dy/dt=2t: analitis = {y1_analytical(3)}")
print(f"y(1) untuk dy/dt=-2y: analitis = {y2_analytical(1):.4f}")

import math
 
def euler_method(f, t0, y0, t_end, h):
    """Euler method untuk ODE: dy/dt = f(t, y)."""
    t_values = [t0]
    y_values = [y0]
 
    t, y = t0, y0
    while t < t_end:
        y = y + h * f(t, y)
        t = t + h
        t_values.append(t)
        y_values.append(y)
 
    return t_values, y_values
 
# Gerak projectile: dy/dt = v₀ - g×t (vertikal)
g = 9.81
v0 = 50  # m/s
 
def projectile_vertical(t, y):
    """dy/dt = v₀ - g×t (height,忽略 drag)."""
    return v0 - g * t
 
# Solusi analitis: y(t) = v₀t - 0.5gt²
def projectile_analytical(t):
    return v0 * t - 0.5 * g * t**2
 
# Numerik dengan Euler
t_num, y_num = euler_method(projectile_vertical, 0, 0, 10, 0.1)
 
# Bandingkan
print("Projectile motion: v₀=50 m/s, g=9.81 m/s²")
print(f"{'t':>5} {'Euler':>8} {'Exact':>8} {'Error':>8}")
for i in range(0, len(t_num), 10):
    t = t_num[i]
    y_euler = y_num[i]
    y_exact = projectile_analytical(t)
    error = abs(y_euler - y_exact)
    print(f"{t:>5.1f} {y_euler:>8.2f} {y_exact:>8.2f} {error:>8.4f}")

print("\nPengaruh step size terhadap akurasi (t=5s):")
exact = projectile_analytical(5)
 
for h in [1.0, 0.5, 0.1, 0.01, 0.001]:
    _, y_vals = euler_method(projectile_vertical, 0, 0, 5, h)
    y_approx = y_vals[-1]
    error = abs(y_approx - exact)
    print(f"  h={h}: y={y_approx:.4f}, error={error:.4f}")

def rk4_method(f, t0, y0, t_end, h):
    """Runge-Kutta 4th order — lebih akurat dari Euler."""
    t_values = [t0]
    y_values = [y0]
 
    t, y = t0, y0
    while t < t_end:
        k1 = f(t, y)
        k2 = f(t + h/2, y + h*k1/2)
        k3 = f(t + h/2, y + h*k2/2)
        k4 = f(t + h, y + h*k3)
 
        y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        t = t + h
        t_values.append(t)
        y_values.append(y)
 
    return t_values, y_values
 
# Bandingkan Euler vs RK4
print("Euler vs RK4 (h=0.1, t=5):")
exact = projectile_analytical(5)
 
_, y_euler = euler_method(projectile_vertical, 0, 0, 5, 0.1)
_, y_rk4 = rk4_method(projectile_vertical, 0, 0, 5, 0.1)
 
print(f"  Euler: {y_euler[-1]:.4f} (error: {abs(y_euler[-1]-exact):.4f})")
print(f"  RK4:   {y_rk4[-1]:.4f} (error: {abs(y_rk4[-1]-exact):.4f})")

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
 
def lotka_volterra(t, state, alpha=1.1, beta=0.4, delta=0.1, gamma=0.4):
    """Predator-prey: dx/dt = αx - βxy, dy/dt = δxy - γy."""
    x, y = state
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]
 
def rk4_system(f, t0, y0, t_end, h, **params):
    """RK4 untuk sistem ODE."""
    t_values = [t0]
    y_values = [list(y0)]
 
    t = t0
    y = list(y0)
    while t < t_end:
        k1 = f(t, y, **params)
        k2 = f(t+h/2, [y[i]+h*k1[i]/2 for i in range(len(y))], **params)
        k3 = f(t+h/2, [y[i]+h*k2[i]/2 for i in range(len(y))], **params)
        k4 = f(t+h, [y[i]+h*k3[i] for i in range(len(y))], **params)
 
        y = [y[i] + (h/6)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]) for i in range(len(y))]
        t = t + h
        t_values.append(t)
        y_values.append(list(y))
 
    return t_values, y_values
 
# Simulasi
t, states = rk4_system(lotka_volterra, 0, [10, 5], 50, 0.01)
prey = [s[0] for s in states]
predator = [s[1] for s in states]
 
print("Predator-Prey (Lotka-Volterra):")
print(f"  Initial: prey=10, predator=5")
print(f"  Final (t=50): prey={prey[-1]:.1f}, predator={predator[-1]:.1f}")
print(f"  Max prey: {max(prey):.1f}, Max predator: {max(predator):.1f}")

import math
 
def projectile_2d_euler(v0, angle_deg, g=9.81, h=0.01, t_end=10):
    """Simulasi projectile 2D dengan Euler method."""
    angle = math.radians(angle_deg)
    vx = v0 * math.cos(angle)
    vy = v0 * math.sin(angle)
 
    x, y = 0, 0
    trajectory = [(x, y)]
 
    t = 0
    while t < t_end and y >= 0:
        # dy/dt = vy, dx/dt = vx
        # dvy/dt = -g, dvx/dt = 0
        x = x + h * vx
        y = y + h * vy
        vy = vy + h * (-g)
        t += h
        if y >= 0:
            trajectory.append((x, y))
 
    return trajectory
 
traj = projectile_2d_euler(50, 45)
print(f"Projectile (v₀=50, 45°):")
print(f"  Landing: x={traj[-1][0]:.1f}m")
print(f"  Max height: {max(p[1] for p in traj):.1f}m")
print(f"  Data points: {len(traj)}")