main.pytxcode

=>PYTHONTEX#py#default#default#0#code#####52#
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from scipy.integrate import odeint, solve_ivp
from scipy.special import factorial
from pathlib import Path

# Set consistent style and random seed
plt.style.use('seaborn-v0_8-whitegrid')
np.random.seed(42)

# Figure settings
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['figure.dpi'] = 150
plt.rcParams['savefig.bbox'] = 'tight'
plt.rcParams['savefig.pad_inches'] = 0.1

# Ensure figures directory exists
Path('figures').mkdir(parents=True, exist_ok=True)
=>PYTHONTEX#py#default#default#1#code#####153#
# First-order ODE examples and solutions
print("Ordinary Differential Equations Analysis:")

# Example 1: dy/dt = -2y + sin(t), y(0) = 1
def ode_example1(y, t):
    """First-order linear ODE: dy/dt = -2y + sin(t)"""
    return -2*y + np.sin(t)

# Analytical solution for comparison
def analytical_solution1(t):
    """Analytical solution: y = (1/5)(sin(t) - 2*cos(t)) + (7/5)*exp(-2*t)"""
    return (1/5)*(np.sin(t) - 2*np.cos(t)) + (7/5)*np.exp(-2*t)

# Time points
t = np.linspace(0, 3, 100)
y0 = 1  # Initial condition

# Numerical solution
y_numerical = odeint(ode_example1, y0, t).flatten()

# Analytical solution
y_analytical = analytical_solution1(t)

# Create visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot solutions comparison
ax1.plot(t, y_numerical, 'b-', linewidth=2, label='Numerical Solution')
ax1.plot(t, y_analytical, 'r--', linewidth=2, label='Analytical Solution')
ax1.set_xlabel('Time t', fontsize=12)
ax1.set_ylabel('y(t)', fontsize=12)
ax1.set_title('First-Order Linear ODE Solution', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot error
error = np.abs(y_numerical - y_analytical)
ax2.semilogy(t, error, 'g-', linewidth=2)
ax2.set_xlabel('Time t', fontsize=12)
ax2.set_ylabel('Absolute Error', fontsize=12)
ax2.set_title('Numerical vs Analytical Error', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('figures/ode_first_order.pdf', dpi=300, bbox_inches='tight')
plt.close()

print(f"Maximum error: {np.max(error):.2e}")
print(r"First-order ODE analysis saved to figures/ode\_first\_order.pdf")
=>PYTHONTEX#py#default#default#2#code#####222#
# Second-order ODE: Damped harmonic oscillator
def harmonic_oscillator(y, t, gamma, omega0):
    """
    Damped harmonic oscillator: x'' + 2*gamma*x' + omega0^2*x = 0
    State vector: y = [x, x']
    """
    x, x_dot = y
    x_ddot = -2*gamma*x_dot - omega0**2*x
    return [x_dot, x_ddot]

# Parameters for different damping regimes
omega0 = 1.0  # Natural frequency
gammas = [0.0, 0.2, 0.5, 1.0, 2.0]  # Different damping coefficients
labels = ['Undamped', 'Underdamped', 'Underdamped', 'Critical', 'Overdamped']

# Time array
t = np.linspace(0, 10, 1000)

# Initial conditions: x(0) = 1, x'(0) = 0
y0 = [1, 0]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Solve and plot for different damping values
for gamma, label in zip(gammas, labels):
    sol = odeint(harmonic_oscillator, y0, t, args=(gamma, omega0))
    x = sol[:, 0]  # Position
    x_dot = sol[:, 1]  # Velocity

    ax1.plot(t, x, linewidth=2, label=f'{label} (γ={gamma})')

    # Phase space plot (position vs velocity)
    ax2.plot(x, x_dot, linewidth=2, label=f'{label} (γ={gamma})')

ax1.set_xlabel('Time t', fontsize=12)
ax1.set_ylabel('Position x(t)', fontsize=12)
ax1.set_title('Damped Harmonic Oscillator', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Velocity dx/dt', fontsize=12)
ax2.set_title('Phase Space Trajectories', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('figures/harmonic_oscillator.pdf', dpi=300, bbox_inches='tight')
plt.close()

print(r"Harmonic oscillator analysis saved to figures/harmonic\_oscillator.pdf")
=>PYTHONTEX#py#default#default#3#code#####299#
# Heat equation numerical solution using finite differences
def heat_equation_fd(nx, nt, alpha, L, T):
    """
    Solve 1D heat equation using finite difference method
    nx: number of spatial points
    nt: number of time steps
    alpha: thermal diffusivity
    L: spatial domain length
    T: total time
    """

    # Grid setup
    dx = L / (nx - 1)
    dt = T / (nt - 1)
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)

    # Stability condition for explicit scheme
    r = alpha * dt / dx**2
    print(f"Stability parameter r = {r:.3f} (should be $\\leq$ 0.5)")

    # Initialize temperature array
    u = np.zeros((nt, nx))

    # Initial condition: Gaussian temperature distribution
    u[0, :] = np.exp(-((x - L/2) / (L/10))**2)

    # Boundary conditions: u(0,t) = u(L,t) = 0
    u[:, 0] = 0
    u[:, -1] = 0

    # Time stepping using explicit finite difference
    for n in range(0, nt-1):
        for i in range(1, nx-1):
            u[n+1, i] = u[n, i] + r * (u[n, i+1] - 2*u[n, i] + u[n, i-1])

    return x, t, u

# Parameters
nx = 50  # Number of spatial points
nt = 500  # Number of time steps
alpha = 0.01  # Thermal diffusivity
L = 1.0  # Domain length
T = 0.5  # Total time

# Solve heat equation
x, t, u = heat_equation_fd(nx, nt, alpha, L, T)

# Create visualization
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot temperature profiles at different times
time_indices = [0, nt//4, nt//2, 3*nt//4, nt-1]
for idx in time_indices:
    ax1.plot(x, u[idx, :], linewidth=2, label=f't = {t[idx]:.3f}')

ax1.set_xlabel('Position x', fontsize=12)
ax1.set_ylabel('Temperature u(x,t)', fontsize=12)
ax1.set_title('Heat Equation: Temperature Profiles', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Contour plot showing temperature evolution
X, T_mesh = np.meshgrid(x, t)
contour = ax2.contourf(X, T_mesh, u, levels=20, cmap='hot')
ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Time t', fontsize=12)
ax2.set_title('Heat Equation: Temperature Evolution', fontsize=14, fontweight='bold')
plt.colorbar(contour, ax=ax2, label='Temperature')

plt.tight_layout()
plt.savefig('figures/heat_equation.pdf', dpi=300, bbox_inches='tight')
plt.close()

print(r"Heat equation analysis saved to figures/heat\_equation.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|