|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +""" |
| 3 | +Model: Lotka-Volterra (Predator-Prey) |
| 4 | +Generated from: Voltera.ipynb |
| 5 | +Author: Charles Pimenta |
| 6 | + |
| 7 | +Exported at: 2025-10-31T12:40:21 |
| 8 | +
|
| 9 | +Note: |
| 10 | +- IPython magics ('%...') and shell lines ('!...') were commented out for .py validity. |
| 11 | +- Each cell is separated by a banner indicating its original index and nearest preceding markdown heading. |
| 12 | +""" |
| 13 | + |
| 14 | + |
| 15 | +#============================================================================== |
| 16 | +# Cell 3 | Section: Implementação em Python |
| 17 | +#============================================================================== |
| 18 | +from dataclasses import dataclass |
| 19 | +import numpy as np |
| 20 | +from scipy.integrate import solve_ivp |
| 21 | +from scipy.integrate import solve_ivp |
| 22 | +import matplotlib.pyplot as plt |
| 23 | + |
| 24 | + |
| 25 | +@dataclass |
| 26 | +class LVParams: |
| 27 | + """Parâmetros do modelo de Lotka–Volterra.""" |
| 28 | + alpha: float = 1.0 # taxa de crescimento das presas (sem predadores) [%] |
| 29 | + beta: float = 0.1 # taxa de predação (acoplamento X–Y) [%] |
| 30 | + gamma: float = 1.5 # taxa de mortalidade natural dos predadores (sem presas) [%] |
| 31 | + delta: float = 0.075 # eficiência de conversão de presas em predadores ou taxa de reprodução dos predadores [%] |
| 32 | + |
| 33 | +# Função diferencial |
| 34 | +def lv_ode(t: float, z: np.ndarray, p: LVParams) -> list: |
| 35 | + """ |
| 36 | + Campo vetorial do sistema de Lotka–Volterra. |
| 37 | + z = [X, Y], sendo X presas e Y predadores. |
| 38 | + """ |
| 39 | + x, y = z |
| 40 | + dxdt = p.alpha * x - p.beta * x * y |
| 41 | + dydt = p.delta * x * y - p.gamma * y |
| 42 | + return [dxdt, dydt] |
| 43 | + |
| 44 | +def simulate_lv(params: LVParams, |
| 45 | + z0=(10.0, 5.0), |
| 46 | + t_span=(0.0, 50.0), |
| 47 | + n_points=1000, |
| 48 | + rtol=1e-8, atol=1e-10): |
| 49 | + """ |
| 50 | + Integra numericamente o sistema e retorna tempos e trajetórias (X(t), Y(t)). |
| 51 | + """ |
| 52 | + t_eval = np.linspace(t_span[0], t_span[1], n_points) |
| 53 | + sol = solve_ivp(lambda t, z: lv_ode(t, z, params), |
| 54 | + t_span=t_span, y0=np.array(z0, dtype=float), |
| 55 | + t_eval=t_eval, rtol=rtol, atol=atol, vectorized=False) |
| 56 | + if not sol.success: |
| 57 | + raise RuntimeError(f"Falha na integração: {sol.message}") |
| 58 | + return sol.t, sol.y[0], sol.y[1] |
| 59 | + |
| 60 | +def equilibrium_points(params: LVParams): |
| 61 | + """ |
| 62 | + Retorna os dois pontos de equilíbrio: |
| 63 | + E0 = (0, 0) |
| 64 | + E1 = (gamma/delta, alpha/beta) |
| 65 | + """ |
| 66 | + e0 = (0.0, 0.0) |
| 67 | + e1 = (params.gamma / params.delta, params.alpha / params.beta) |
| 68 | + return e0, e1 |
| 69 | + |
| 70 | +def plot_time_series(t, X, Y, params: LVParams, z0,filename='lv_time_series.png'): |
| 71 | + plt.figure(figsize=(9, 4.5)) |
| 72 | + plt.plot(t, X, label='Presas (X)') |
| 73 | + plt.plot(t, Y, label='Predadores (Y)') |
| 74 | + plt.xlabel('Tempo') |
| 75 | + plt.ylabel('População') |
| 76 | + plt.title('Dinâmica Predador–Presa (Lotka–Volterra)') |
| 77 | + plt.legend(loc='best') |
| 78 | + plt.grid(True, alpha=0.3) |
| 79 | + plt.tight_layout() |
| 80 | + plt.savefig(filename, dpi=160) |
| 81 | + plt.show() |
| 82 | + |
| 83 | +def plot_phase_plane(X, Y, params: LVParams, z0,filename='lv_phase_plane.png'): |
| 84 | + e0, e1 = equilibrium_points(params) |
| 85 | + |
| 86 | + plt.figure(figsize=(5.5, 5.5)) |
| 87 | + plt.plot(X, Y, linewidth=2.0, label='Trajetória') |
| 88 | + plt.scatter([z0[0]], [z0[1]], marker='p', s=120, color='red', label='Condição inicial em $t_0$') |
| 89 | + plt.scatter([e0[0]], [e0[1]], marker='x', s=80, label='Equilíbrio (0,0)') |
| 90 | + plt.scatter([e1[0]], [e1[1]], marker='o', s=60, label=r'Equilíbrio $(\gamma/\delta,\ \alpha/\beta)$') |
| 91 | + plt.xlabel('Presas (X)') |
| 92 | + plt.ylabel('Predadores (Y)') |
| 93 | + plt.title('Plano de Fase') |
| 94 | + plt.grid(True, alpha=0.3) |
| 95 | + plt.legend(loc='best') |
| 96 | + plt.tight_layout() |
| 97 | + plt.savefig(filename, dpi=160) |
| 98 | + plt.show() |
| 99 | + |
| 100 | +#============================================================================== |
| 101 | +# Cell 5 | Section: Gráficos: séries temporais e plano de fase |
| 102 | +#============================================================================== |
| 103 | +params = LVParams(alpha=1.0, beta=0.1, gamma=1.5, delta=0.075) |
| 104 | +z0 = (10.0, 5.0) |
| 105 | +t, X, Y = simulate_lv(params, z0=z0, t_span=(0, 25), n_points=1200) |
| 106 | +plot_time_series(t, X, Y, params, z0,filename='lv_time_series_1.png') |
| 107 | +plot_phase_plane(X, Y, params, z0,filename='lv_phase_plane_1.png') |
| 108 | + |
| 109 | +#============================================================================== |
| 110 | +# Cell 7 | Section: Análise da sensibilidade a variação de parâmetros (variação de $\alpha$ e $\gamma$) |
| 111 | +#============================================================================== |
| 112 | +def sweep_params(base: LVParams, alphas = [0.8, 1.0, 1.2], gammas = [1.2, 1.5, 1.8], z0=(10.0, 5.0), t_span=(0, 40), n_points=800): |
| 113 | + # Variação de alpha (presas) |
| 114 | + plt.figure(figsize=(9, 4.5)) |
| 115 | + for a in alphas: |
| 116 | + p = LVParams(alpha=a, beta=base.beta, gamma=base.gamma, delta=base.delta) |
| 117 | + t, X, Y = simulate_lv(p, z0=z0, t_span=t_span, n_points=n_points) |
| 118 | + plt.plot(t, X, label=fr'$\alpha={a:.1f}$') |
| 119 | + plt.xlabel('Tempo'); plt.ylabel('Presas (X)') |
| 120 | + plt.title('Sensibilidade: variação de α (presas)') |
| 121 | + plt.grid(True, alpha=0.3); plt.legend() |
| 122 | + plt.tight_layout(); plt.savefig('lv_sens_alpha.png', dpi=160); plt.show() |
| 123 | + |
| 124 | + # Plano de fase |
| 125 | + plt.figure(figsize=(5.5, 5.5)) |
| 126 | + for a in alphas: |
| 127 | + p = LVParams(alpha=a, beta=base.beta, gamma=base.gamma, delta=base.delta) |
| 128 | + t, X, Y = simulate_lv(p, z0=z0, t_span=t_span, n_points=n_points) |
| 129 | + plt.plot(X, Y, linewidth=2.0, label=fr'Trajetória: $\alpha={a:.1f}$') |
| 130 | + plt.xlabel('Presas (X)') |
| 131 | + plt.ylabel('Predadores (Y)') |
| 132 | + plt.title('Plano de Fase') |
| 133 | + plt.grid(True, alpha=0.3) |
| 134 | + plt.legend(loc='upper right') |
| 135 | + plt.tight_layout() |
| 136 | + plt.savefig('lv_phase_plane_var_alpha.png', dpi=160) |
| 137 | + plt.show() |
| 138 | + |
| 139 | + # Variação de gamma (predadores) |
| 140 | + plt.figure(figsize=(9, 4.5)) |
| 141 | + for g in gammas: |
| 142 | + p = LVParams(alpha=base.alpha, beta=base.beta, gamma=g, delta=base.delta) |
| 143 | + t, X, Y = simulate_lv(p, z0=z0, t_span=t_span, n_points=n_points) |
| 144 | + plt.plot(t, Y, label=fr'$\gamma={g:.1f}$') |
| 145 | + plt.xlabel('Tempo'); plt.ylabel('Predadores (Y)') |
| 146 | + plt.title('Sensibilidade: variação de γ (predadores)') |
| 147 | + plt.grid(True, alpha=0.3); plt.legend() |
| 148 | + plt.tight_layout(); plt.savefig('lv_sens_gamma.png', dpi=160); plt.show() |
| 149 | + |
| 150 | + # Plano de fase |
| 151 | + plt.figure(figsize=(5.5, 5.5)) |
| 152 | + for g in gammas: |
| 153 | + p = LVParams(alpha=base.alpha, beta=base.beta, gamma=g, delta=base.delta) |
| 154 | + t, X, Y = simulate_lv(p, z0=z0, t_span=t_span, n_points=n_points) |
| 155 | + plt.plot(X, Y, linewidth=2.0, label=fr'Trajetória: $\gamma={g:.1f}$') |
| 156 | + plt.xlabel('Presas (X)') |
| 157 | + plt.ylabel('Predadores (Y)') |
| 158 | + plt.title('Plano de Fase') |
| 159 | + plt.grid(True, alpha=0.3) |
| 160 | + plt.legend(loc='upper right') |
| 161 | + plt.tight_layout() |
| 162 | + plt.savefig('lv_phase_plane_var_gamma.png', dpi=160) |
| 163 | + plt.show() |
| 164 | + |
| 165 | + |
| 166 | +base = LVParams(alpha=1.0, beta=0.1, gamma=1.5, delta=0.075) |
| 167 | +alphas = [0.8, 1.0, 1.2] |
| 168 | +gammas = [1.2, 1.5, 1.8] |
| 169 | +sweep_params(base,alphas,gammas) |
| 170 | + |
| 171 | +#============================================================================== |
| 172 | +# Cell 9 | Section: Validação rápida dos equilíbrios |
| 173 | +#============================================================================== |
| 174 | +params = LVParams() |
| 175 | +e0, e1 = equilibrium_points(params) |
| 176 | +params = LVParams(alpha=1.0, beta=0.1, gamma=1.5, delta=0.075) |
| 177 | + |
| 178 | +# Derivadas em cada equilíbrio (devem ser ~0): |
| 179 | +for name, pt in [('E0', e0), ('E1', e1)]: |
| 180 | + dx, dy = lv_ode(0.0, np.array(pt, dtype=float), params) |
| 181 | + print(f"{name} = {pt} -> (dX/dt, dY/dt) = ({dx:.6f}, {dy:.6f})") |
| 182 | + |
| 183 | + t, X, Y = simulate_lv(params, z0=pt, t_span=(0, 25), n_points=1200) |
| 184 | + |
| 185 | + filename = f'lv_time_series_eq_{pt[0]:.0f}_{pt[1]:.0f}.png' |
| 186 | + plot_time_series(t, X, Y, params, pt,filename=filename) |
| 187 | + filename = f'lv_phase_plane_eq_{pt[0]:.0f}_{pt[1]:.0f}.png' |
| 188 | + plot_phase_plane(X, Y, params, pt,filename=filename) |
| 189 | + |
0 commit comments