# Set a dark scientific theme
theme(:dark)
# ------------------------------------------------------------------------------
# 1. Configuration & Parameters
# ------------------------------------------------------------------------------
module Config
# System Parameters (Lotka-Volterra)
# α: prey growth, β: predation rate, δ: predator death, γ: predator growth
const P_SYS = (α = 1.1, β = 0.4, δ = 0.4, γ = 0.1)
# Domain Bounds
const X_LIM = (0.0, 15.0) # Prey population range
const Y_LIM = (0.0, 10.0) # Predator population range
# Grid Density for Vector Field
const GRID_DENSITY = 25 # Number of points per axis
# Visualization Styling
const ARROW_SCALE = 0.35 # visual length of arrows
const TRAJECTORY_TIME = (0.0, 40.0)
const TRAJ_COLOR = :cyan
const FIELD_COLOR_GRAD = :inferno
end
# ------------------------------------------------------------------------------
# 2. ODE Definition
# ------------------------------------------------------------------------------
"""
system_dynamics!(du, u, p, t)
The differential equation definition.
u[1] = x (Prey), u[2] = y (Predator)
"""
function system_dynamics!(du, u, p, t)
x, y = u
α, β, δ, γ = p
# dx/dt = αx - βxy
du[1] = α * x - β * x * y
# dy/dt = -δy + γxy
du[2] = -δ * y + γ * x * y
end
# ------------------------------------------------------------------------------
# 3. Vector Field Computation
# ------------------------------------------------------------------------------
"""
compute_vector_field(x_range, y_range, params)
Calculates the vector field (U, V) and Magnitude (M) over a grid.
"""
function compute_vector_field(x_rng, y_rng, p)
# Create meshgrid
x_grid = repeat(x_rng', length(y_rng), 1)
y_grid = repeat(y_rng, 1, length(x_rng))
u_grid = zeros(size(x_grid))
v_grid = zeros(size(y_grid))
mag_grid = zeros(size(x_grid))
du = zeros(2)
# Compute vectors at each grid point
for i in eachindex(x_grid)
state = [x_grid[i], y_grid[i]]
system_dynamics!(du, state, p, 0.0) # t=0 for autonomous systems
# Calculate magnitude
magnitude = norm(du)
mag_grid[i] = magnitude
# Normalize for visualization (prevent giant arrows)
# We perform a "soft" normalization to keep direction clear but bound length
if magnitude > 0
scaling = Config.ARROW_SCALE / (magnitude^0.4) # non-linear scaling for better visuals
u_grid[i] = du[1] * scaling
v_grid[i] = du[2] * scaling
end
end
return x_grid, y_grid, u_grid, v_grid, mag_grid
end
# ------------------------------------------------------------------------------
# 4. Trajectory Simulation
# ------------------------------------------------------------------------------
"""
solve_trajectory(u0, params)
Solves the ODE for a specific initial condition using DifferentialEquations.jl
"""
function solve_trajectory(u0, p)
prob = ODEProblem(system_dynamics!, u0, Config.TRAJECTORY_TIME, p)
# Tsit5 is a standard efficient solver for non-stiff ODEs
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
return sol
end
# ------------------------------------------------------------------------------
# 5. Main Visualization Routine
# ------------------------------------------------------------------------------
function main()
println("Generating Vector Field Visualization...")
# A. Setup Grid
xs = range(Config.X_LIM[1], Config.X_LIM[2], length=Config.GRID_DENSITY)
ys = range(Config.Y_LIM[1], Config.Y_LIM[2], length=Config.GRID_DENSITY)
params =Values = (Config.P_SYS.α, Config.P_SYS.β, Config.P_SYS.δ, Config.P_SYS.γ)
# B. Compute Field
X, Y, U, V, Mag = compute_vector_field(xs, ys, params)
# C. Initialize Plot
plt = plot(
title = "Phase Portrait: Lotka-Volterra System",
xlabel = "Prey Population (x)",
ylabel = "Predator Population (y)",
xlims = Config.X_LIM,
ylims = Config.Y_LIM,
aspect_ratio = :equal,
legend = :topright,
grid = false,
framestyle = :box,
size = (800, 600)
)
# D. Plot Vector Field (Quiver)
# We flatten the arrays because quiver expects vectors, not matrices
quiver!(
plt,
vec(X), vec(Y),
quiver = (vec(U), vec(V)),
color = :white,
alpha = 0.4,
linewidth = 1.0,
label = "Vector Field"
)
# E. Solve and Plot Trajectories
# Define interesting initial conditions
initial_conditions = [
[2.0, 2.0],
[6.0, 4.0],
[10.0, 8.0]
]
for (i, u0) in enumerate(initial_conditions)
sol = solve_trajectory(u0, params)
plot!(plt, sol, vars=(1, 2),
lw=2.5,
color=Config.TRAJ_COLOR,
alpha=0.8,
label = i==1 ? "Trajectories" : "") # Only label the first one
# Mark start points
scatter!(plt, [u0[1]], [u0[2]],
color=:red, marker=:circle, markersize=4, label=nothing)
end
# F. Add Fixed Point (Equilibrium)
# For L-V: x* = δ/γ, y* = α/β
eq_x = Config.P_SYS.δ / Config.P_SYS.γ
eq_y = Config.P_SYS.α / Config.P_SYS.β
scatter!(plt, [eq_x], [eq_y],
color=:yellow, marker=:star5, markersize=8,
label="Equilibrium")
# Display final plot
display(plt)
println("Visualization Complete.")
end
# Execute
main()
