ubuntu2404
=>PYTHONTEX#py#default#default#0#code#####57#
import matplotlib
matplotlib.use("Agg") # non-GUI backend; set before importing pyplot
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
np.random.seed(42)
def lotka_volterra_system(alpha=1.0, beta=0.5, gamma=0.2, delta=0.8, y0=None, t_end=40, n_t=2000):
"""Lotka–Volterra predator–prey model with analysis and plotting helpers."""
def f(t, z):
x, y = z
dxdt = alpha * x - beta * x * y
dydt = gamma * x * y - delta * y
return [dxdt, dydt]
# Equilibrium
x_eq = delta / gamma
y_eq = alpha / beta
# Pick IC away from equilibrium to show oscillations
if y0 is None:
y0 = [1.25 * x_eq, 0.75 * y_eq] # off-equilibrium
# Solve
t_span = (0, t_end)
t_eval = np.linspace(0, t_end, n_t)
sol = solve_ivp(f, t_span, y0, t_eval=t_eval, method='RK45', rtol=1e-8, atol=1e-10)
# Peak timing (prey leads predator)
x_t, y_t = sol.y
peaks_x, _ = find_peaks(x_t)
peaks_y, _ = find_peaks(y_t)
t_x_peak = sol.t[peaks_x[0]] if len(peaks_x) else np.nan
t_y_peak = sol.t[peaks_y[0]] if len(peaks_y) else np.nan
lead = t_y_peak - t_x_peak # >0 means prey peaks first
# Period estimate from prey peaks
T = sol.t[peaks_x[1]] - sol.t[peaks_x[0]] if len(peaks_x) > 1 else np.nan
print("Lotka–Volterra Model (predator–prey):")
print(f" Parameters: alpha={alpha:.3f}, beta={beta:.3f}, gamma={gamma:.3f}, delta={delta:.3f}")
print(f" Equilibrium: $(x^*,y^*)=({x_eq:.3f},{y_eq:.3f})$")
print(f" Initial conditions: $x_0={y0[0]:.3f}$, $y_0={y0[1]:.3f}$")
if np.isfinite(lead):
who_leads = "Prey lead predator" if lead > 0 else "Predator lead prey"
print(f" {who_leads} by about {abs(lead):.2f} time units; estimated period $T\\approx{T:.2f}$")
return sol.t, sol.y, (x_eq, y_eq), (t_x_peak, t_y_peak), T
def sir_epidemic_model():
"""SIR epidemic model (kept for comparison)."""
def sir_system(t, y, beta, gamma, N):
S, I, R = y
dSdt = -beta * S * I / N
dIdt = beta * S * I / N - gamma * I
dRdt = gamma * I
return [dSdt, dIdt, dRdt]
beta = 0.5
gamma = 0.1
N = 1000
R0 = beta / gamma
S0, I0, R0_init = 999, 1, 0
y0 = [S0, I0, R0_init]
t_span = (0, 100)
t_eval = np.linspace(0, 100, 1000)
sol = solve_ivp(lambda t, y: sir_system(t, y, beta, gamma, N),
t_span, y0, t_eval=t_eval, method='RK45')
peak_idx = np.argmax(sol.y[1])
peak_time = sol.t[peak_idx]
peak_infections = sol.y[1][peak_idx]
print("\nSIR Epidemic Model:")
print(f" beta={beta}, gamma={gamma}, N={N}, $R_0={R0:.2f}$")
print(f" Peak infections: {peak_infections:.0f} at $t={peak_time:.1f}$")
return sol.t, sol.y, R0, peak_time, peak_infections
# Plot style
try:
plt.style.use('seaborn-v0_8-colorblind')
except Exception:
plt.style.use('seaborn-colorblind')
plt.rcParams.update({
"figure.dpi": 120,
"axes.grid": True,
"grid.alpha": 0.3,
"axes.spines.top": False,
"axes.spines.right": False,
})
# Analyze LV and SIR
lv_time, lv_solution, (x_eq, y_eq), (t_x_peak, t_y_peak), T = lotka_volterra_system()
sir_time, sir_solution, R0, peak_t, peak_I = sir_epidemic_model()
# Improved LV plots: time series with peak annotations
x_lv, y_lv = lv_solution
fig1, ax = plt.subplots(1, 1, figsize=(9, 4), constrained_layout=True)
ax.plot(lv_time, x_lv, label='Prey $x(t)$', color='#1f77b4', lw=2)
ax.plot(lv_time, y_lv, label='Predator $y(t)$', color='#ff7f0e', lw=2)
ax.axhline(x_eq, color='gray', ls='--', lw=1, label='$x^*$ nullcline ($y=\\alpha/\\beta$)')
ax.axhline(y_eq, color='gray', ls=':', lw=1, label='$y^*$ nullcline ($x=\\delta/\\gamma$)')
# Mark first peaks and lead–lag
from math import isfinite
if isfinite(t_x_peak):
ax.axvline(t_x_peak, color='#1f77b4', ls='--', lw=1, alpha=0.8)
ax.plot([t_x_peak], [np.interp(t_x_peak, lv_time, x_lv)], 'o', color='#1f77b4', ms=4)
if isfinite(t_y_peak):
ax.axvline(t_y_peak, color='#ff7f0e', ls='--', lw=1, alpha=0.8)
ax.plot([t_y_peak], [np.interp(t_y_peak, lv_time, y_lv)], 'o', color='#ff7f0e', ms=4)
if isfinite(t_x_peak) and isfinite(t_y_peak):
lead = t_y_peak - t_x_peak
txt = f"Prey peak leads predator by {lead:.2f}"
ax.annotate(txt, xy=(0.02, 0.95), xycoords='axes fraction',
ha='left', va='top', fontsize=10,
bbox=dict(boxstyle='round,pad=0.25', fc='white', ec='0.7', alpha=0.9))
ax.set_xlabel('$t$')
ax.set_ylabel('Population')
ax.set_title('Lotka–Volterra predator–prey: oscillations and phase lead')
ax.legend(loc='best')
# LV phase plane with vector field, nullclines, and multiple trajectories
fig2, ax2 = plt.subplots(1, 1, figsize=(6.5, 5.5), constrained_layout=True)
# Vector field (streamplot)
x_max = 2.5 * x_eq
y_max = 2.5 * y_eq
X, Y = np.meshgrid(np.linspace(0, x_max, 32), np.linspace(0, y_max, 32))
U = (1.0 * X) - (0.5 * X * Y) # using alpha=1.0, beta=0.5
V = (0.2 * X * Y) - (0.8 * Y) # using gamma=0.2, delta=0.8
speed = np.hypot(U, V)
ax2.streamplot(X, Y, U, V, density=1.2, color='#cccccc', arrowsize=1.2, linewidth=1)
# Nullclines
ax2.axhline(y_eq, color='gray', ls=':', lw=1.5, label='$dx/dt=0: y=\\alpha/\\beta$')
ax2.axvline(x_eq, color='gray', ls='--', lw=1.5, label='$dy/dt=0: x=\\delta/\\gamma$')
# Multiple trajectories around equilibrium
ics = [
[1.15 * x_eq, 0.85 * y_eq],
[0.75 * x_eq, 1.35 * y_eq],
[1.60 * x_eq, 0.55 * y_eq],
]
colors = ['#2ca02c', '#9467bd', '#8c564b']
t_eval = np.linspace(0, 60, 4000)
def solve_ic(ic):
def f(t, z):
x, y = z
return [1.0 * x - 0.5 * x * y, 0.2 * x * y - 0.8 * y]
return solve_ivp(f, (0, t_eval[-1]), ic, t_eval=t_eval, rtol=1e-8, atol=1e-10)
for (ic, c) in zip(ics, colors):
sol_ic = solve_ic(ic)
ax2.plot(sol_ic.y[0], sol_ic.y[1], color=c, lw=2, label=f'Traj from $(x_0,y_0)=({ic[0]:.1f},{ic[1]:.1f})$')
# Highlight equilibrium
ax2.plot([x_eq], [y_eq], 'r*', ms=12, label='Equilibrium $(x^*,y^*)$')
ax2.set_xlim(0, x_max)
ax2.set_ylim(0, y_max)
ax2.set_xlabel('Prey $x$')
ax2.set_ylabel('Predator $y$')
ax2.set_title('Lotka–Volterra: phase plane, nullclines, and trajectories')
ax2.legend(loc='best', fontsize=9)
# SIR model plot (unchanged logic, improved styling)
S, I, R = sir_solution
fig3, ax3 = plt.subplots(figsize=(8, 4), constrained_layout=True)
ax3.plot(sir_time, S, label='$S(t)$', color='#1f77b4', lw=2)
ax3.plot(sir_time, I, label='$I(t)$', color='#d62728', lw=2)
ax3.plot(sir_time, R, label='$R(t)$', color='#2ca02c', lw=2)
ax3.axvline(peak_t, color='k', ls='--', lw=1, label=f'Peak $I$ at $t={peak_t:.1f}$')
ax3.plot([peak_t], [peak_I], 'ko', ms=4)
ax3.set_xlabel('$t$')
ax3.set_ylabel('Population')
ax3.set_title(f'SIR model ($R_0={R0:.2f}$)')
ax3.legend(loc='best')
# Save figures instead of showing
fig1.savefig("lv_timeseries.png", bbox_inches="tight")
fig2.savefig("lv_phaseplane.png", bbox_inches="tight")
fig3.savefig("sir_model.png", bbox_inches="tight")
=>PYTHONTEX#py#default#default#1#code#####265#
def monte_carlo_integration():
"""Demonstrate Monte Carlo integration methods"""
def integrand(x):
"""Test function: f(x) = exp(-x^2)*sin(x)"""
return np.exp(-x**2) * np.sin(x)
# Integration domain
a, b = 0, 2
# Analytical approximation for comparison
from scipy.integrate import quad
analytical_result, _ = quad(integrand, a, b)
# Monte Carlo integration with different sample sizes
sample_sizes = [100, 500, 1000, 5000, 10000, 50000]
mc_estimates = []
mc_errors = []
for n in sample_sizes:
# Generate random samples
x_samples = np.random.uniform(a, b, n)
f_samples = integrand(x_samples)
# Monte Carlo estimate
mc_estimate = (b - a) * np.mean(f_samples)
mc_estimates.append(mc_estimate)
# Error estimate
error = abs(mc_estimate - analytical_result)
mc_errors.append(error)
print(f"\nMonte Carlo Integration:")
print(f" Integrand: exp(-x squared)*sin(x) on [0, 2]")
print(f" Analytical result: {analytical_result:.6f}")
print(f" Sample sizes: {sample_sizes}")
return sample_sizes, mc_estimates, mc_errors, analytical_result
def random_walk_simulation():
"""Simulate 2D random walk and analyze properties"""
n_steps = 10000
n_walks = 100
# Store final positions
final_positions = np.zeros((n_walks, 2))
# Generate random walks
for walk in range(n_walks):
# Random steps: +1 or -1 in each direction
steps_x = np.random.choice([-1, 1], n_steps)
steps_y = np.random.choice([-1, 1], n_steps)
# Cumulative position
position_x = np.cumsum(steps_x)
position_y = np.cumsum(steps_y)
final_positions[walk] = [position_x[-1], position_y[-1]]
# Store one example trajectory
if walk == 0:
example_trajectory = np.column_stack([position_x, position_y])
# Analyze displacement statistics
displacements = np.linalg.norm(final_positions, axis=1)
mean_displacement = np.mean(displacements)
std_displacement = np.std(displacements)
# Theoretical expectation: E[|X_n|] ~ sqrt(n)
theoretical_rms = np.sqrt(n_steps)
print(f"\nRandom Walk Simulation:")
print(f" Steps per walk: {n_steps}")
print(f" Number of walks: {n_walks}")
print(f" Mean displacement: {mean_displacement:.1f}")
print(f" Theoretical RMS: {theoretical_rms:.1f}")
print(f" Displacement std: {std_displacement:.1f}")
return example_trajectory, final_positions, displacements
# Perform Monte Carlo simulations
sample_sizes, mc_est, mc_err, analytical = monte_carlo_integration()
trajectory, final_pos, displacements = random_walk_simulation()
=>PYTHONTEX#py#default#default#2#code#####369#
def flocking_simulation():
"""Implement simplified boids flocking model"""
class Boid:
def __init__(self, x, y, vx, vy):
self.x = x
self.y = y
self.vx = vx
self.vy = vy
def update(self, boids, width, height, dt=0.1):
"""Update boid position and velocity based on flocking rules"""
# Flocking parameters
separation_radius = 5.0
alignment_radius = 10.0
cohesion_radius = 15.0
max_speed = 2.0
separation_force = np.array([0.0, 0.0])
alignment_force = np.array([0.0, 0.0])
cohesion_force = np.array([0.0, 0.0])
neighbors = []
for other in boids:
if other is not self:
dx = other.x - self.x
dy = other.y - self.y
distance = np.sqrt(dx**2 + dy**2)
if distance < cohesion_radius:
neighbors.append(other)
# Separation: steer away from neighbors
if distance < separation_radius and distance > 0:
separation_force[0] -= dx / distance
separation_force[1] -= dy / distance
# Alignment: match neighbor velocities
if distance < alignment_radius:
alignment_force[0] += other.vx
alignment_force[1] += other.vy
# Cohesion: steer toward center of neighbors
cohesion_force[0] += dx
cohesion_force[1] += dy
if len(neighbors) > 0:
# Normalize forces
separation_force *= 1.5
alignment_force /= len(neighbors)
cohesion_force /= len(neighbors)
cohesion_force *= 0.01
# Apply forces
self.vx += separation_force[0] * dt + alignment_force[0] * dt + cohesion_force[0] * dt
self.vy += separation_force[1] * dt + alignment_force[1] * dt + cohesion_force[1] * dt
# Limit speed
speed = np.sqrt(self.vx**2 + self.vy**2)
if speed > max_speed:
self.vx = self.vx / speed * max_speed
self.vy = self.vy / speed * max_speed
# Update position
self.x += self.vx * dt
self.y += self.vy * dt
# Periodic boundary conditions
self.x = self.x % width
self.y = self.y % height
# Initialize simulation
n_boids = 50
width, height = 100, 100
boids = []
for _ in range(n_boids):
x = np.random.uniform(0, width)
y = np.random.uniform(0, height)
vx = np.random.uniform(-1, 1)
vy = np.random.uniform(-1, 1)
boids.append(Boid(x, y, vx, vy))
# Simulate for several time steps
n_steps = 200
positions_history = []
for step in range(n_steps):
# Update all boids
for boid in boids:
boid.update(boids, width, height)
# Record positions every 10 steps
if step % 10 == 0:
positions = [(boid.x, boid.y) for boid in boids]
positions_history.append(positions)
print(f"\nFlocking Simulation:")
print(f" Number of boids: {n_boids}")
print(f" Domain size: {width}×{height}")
print(f" Simulation steps: {n_steps}")
print(f" Recorded snapshots: {len(positions_history)}")
return positions_history, width, height
# Run agent-based modeling simulation
flock_history, domain_width, domain_height = flocking_simulation()
=>PYTHONTEX#py#default#default#3#code#####481#
# Create comprehensive modeling and simulation visualization
import os
os.makedirs('assets', exist_ok=True)
fig, axes = plt.subplots(3, 3, figsize=(18, 15))
fig.suptitle('Mathematical Modeling and Simulation', fontsize=16, fontweight='bold')
# Lotka-Volterra time series
ax1 = axes[0, 0]
ax1.plot(lv_time, lv_solution[0], 'b-', linewidth=2, label='Prey')
ax1.plot(lv_time, lv_solution[1], 'r-', linewidth=2, label='Predator')
ax1.plot([0, lv_time[-1]], [x_eq, x_eq], 'b--', alpha=0.7, label='Prey equilibrium')
ax1.plot([0, lv_time[-1]], [y_eq, y_eq], 'r--', alpha=0.7, label='Predator equilibrium')
ax1.set_xlabel('Time')
ax1.set_ylabel('Population')
ax1.set_title('Lotka-Volterra System')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Lotka-Volterra phase portrait
ax2 = axes[0, 1]
ax2.plot(lv_solution[0], lv_solution[1], 'g-', linewidth=2)
ax2.plot(x_eq, y_eq, 'ko', markersize=8, label='Equilibrium')
ax2.plot(lv_solution[0][0], lv_solution[1][0], 'go', markersize=8, label='Initial')
ax2.set_xlabel('Prey Population')
ax2.set_ylabel('Predator Population')
ax2.set_title('Phase Portrait')
ax2.legend()
ax2.grid(True, alpha=0.3)
# SIR model
ax3 = axes[0, 2]
ax3.plot(sir_time, sir_solution[0], 'b-', linewidth=2, label='Susceptible')
ax3.plot(sir_time, sir_solution[1], 'r-', linewidth=2, label='Infected')
ax3.plot(sir_time, sir_solution[2], 'g-', linewidth=2, label='Recovered')
ax3.axvline(peak_t, color='red', linestyle='--', alpha=0.7, label=f'Peak at t={peak_t:.1f}')
ax3.set_xlabel('Time (days)')
ax3.set_ylabel('Population')
ax3.set_title(f'SIR Epidemic Model (R0={R0:.1f})')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Monte Carlo convergence
ax4 = axes[1, 0]
ax4.loglog(sample_sizes, mc_err, 'bo-', linewidth=2, markersize=6)
ax4.loglog(sample_sizes, [1/np.sqrt(n) for n in sample_sizes], 'r--', linewidth=2, label='1/√n')
ax4.set_xlabel('Sample Size')
ax4.set_ylabel('Integration Error')
ax4.set_title('Monte Carlo Convergence')
ax4.legend()
ax4.grid(True, alpha=0.3)
# Random walk trajectory
ax5 = axes[1, 1]
ax5.plot(trajectory[:1000, 0], trajectory[:1000, 1], 'b-', linewidth=1, alpha=0.7)
ax5.plot(trajectory[0, 0], trajectory[0, 1], 'go', markersize=8, label='Start')
ax5.plot(trajectory[-1, 0], trajectory[-1, 1], 'ro', markersize=8, label='End')
ax5.set_xlabel('x-position')
ax5.set_ylabel('y-position')
ax5.set_title('Random Walk Trajectory')
ax5.legend()
ax5.grid(True, alpha=0.3)
# Displacement distribution
ax6 = axes[1, 2]
ax6.hist(displacements, bins=20, alpha=0.7, density=True, color='skyblue')
ax6.axvline(np.mean(displacements), color='red', linestyle='--', linewidth=2, label='Mean')
ax6.axvline(np.sqrt(len(trajectory)), color='green', linestyle='--', linewidth=2, label='Theoretical')
ax6.set_xlabel('Final Displacement')
ax6.set_ylabel('Density')
ax6.set_title('Displacement Distribution')
ax6.legend()
ax6.grid(True, alpha=0.3)
# Flocking simulation snapshots
ax7 = axes[2, 0]
# Show initial positions
if len(flock_history) > 0:
initial_pos = np.array(flock_history[0])
ax7.scatter(initial_pos[:, 0], initial_pos[:, 1], c='blue', s=20, alpha=0.6)
ax7.set_xlim(0, domain_width)
ax7.set_ylim(0, domain_height)
ax7.set_xlabel('x')
ax7.set_ylabel('y')
ax7.set_title('Flocking: Initial State')
ax7.grid(True, alpha=0.3)
ax8 = axes[2, 1]
# Show final positions
if len(flock_history) > 0:
final_pos = np.array(flock_history[-1])
ax8.scatter(final_pos[:, 0], final_pos[:, 1], c='red', s=20, alpha=0.6)
ax8.set_xlim(0, domain_width)
ax8.set_ylim(0, domain_height)
ax8.set_xlabel('x')
ax8.set_ylabel('y')
ax8.set_title('Flocking: Final State')
ax8.grid(True, alpha=0.3)
# Model comparison summary
ax9 = axes[2, 2]
model_types = ['Deterministic\n(ODE)', 'Stochastic\n(Monte Carlo)', 'Agent-Based\n(Flocking)']
y_pos = np.arange(len(model_types))
# Dummy data for illustration
complexity = [2, 3, 4]
bars = ax9.barh(y_pos, complexity, color=['lightblue', 'lightgreen', 'lightcoral'])
ax9.set_yticks(y_pos)
ax9.set_yticklabels(model_types)
ax9.set_xlabel('Computational Complexity')
ax9.set_title('Modeling Approaches')
plt.tight_layout()
plt.savefig('assets/mathematical_modeling.pdf', dpi=300, bbox_inches='tight')
plt.close()
print("Mathematical modeling analysis saved to assets/mathematical\\_modeling.pdf")
=>PYTHONTEX:SETTINGS#
version=0.18
outputdir=pythontex-files-main
workingdir=.
workingdirset=false
gobble=none
rerun=default
hashdependencies=default
makestderr=false
stderrfilename=full
keeptemps=none
pyfuture=default
pyconfuture=none
pygments=true
pygglobal=:GLOBAL||
fvextfile=-1
pyconbanner=none
pyconfilename=stdin
depythontex=false
pygfamily=py|python3|
pygfamily=pycon|pycon|
pygfamily=sympy|python3|
pygfamily=sympycon|pycon|
pygfamily=pylab|python3|
pygfamily=pylabcon|pycon|