Chapter 2: Master Equation and Stochastic Processes

import numpy as np import matplotlib.pyplot as plt class BasicSolver: """Basic solver for master equations and stochastic processes""" def __init__(self, N=100): self.N = N self.x = np.linspace(0, 10, N) self.dx = self.x[1] - self.x[0] def solve(self, T=1.0, dt=0.01): """Basic time evolution solver""" steps = int(T / dt) solution = np.zeros((steps, self.N)) # Initial conditions solution[0, :] = np.exp(-(self.x - 5)**2) # Time evolution for n in range(steps - 1): # Implement specific algorithm here solution[n+1, :] = solution[n, :] # Placeholder return solution def plot(self, solution): """Result visualization""" fig, ax = plt.subplots(figsize=(10, 6)) im = ax.contourf(self.x, np.arange(len(solution)), solution, levels=20, cmap='viridis') ax.set_xlabel('Spatial coordinate x', fontsize=12) ax.set_ylabel('Time step', fontsize=12) ax.set_title('Master Equation and Stochastic Process Simulation', fontsize=14, fontweight='bold') plt.colorbar(im, ax=ax) return fig # Execution example solver = BasicSolver(N=100) solution = solver.solve(T=1.0, dt=0.01) fig = solver.plot(solution) plt.show()

import numpy as np import matplotlib.pyplot as plt from scipy import sparse from scipy.sparse.linalg import spsolve class AdvancedSolver: """Advanced algorithm implementation""" def __init__(self, N=100, method='implicit'): self.N = N self.method = method self.x = np.linspace(0, 10, N) self.dx = self.x[1] - self.x[0] def build_matrix(self, dt): """Matrix construction (for implicit method)""" N = self.N diag = np.ones(N) off_diag = -0.5 * np.ones(N-1) A = sparse.diags([off_diag, diag, off_diag], [-1, 0, 1], format='csr') return A def solve(self, T=1.0, dt=0.01): """Time evolution solver""" steps = int(T / dt) solution = np.zeros((steps, self.N)) # Initial conditions solution[0, :] = self.initial_condition() if self.method == 'implicit': A = self.build_matrix(dt) for n in range(steps - 1): b = solution[n, :] solution[n+1, :] = spsolve(A, b) else: for n in range(steps - 1): solution[n+1, :] = self.explicit_step(solution[n, :], dt) return solution def initial_condition(self): """Set initial conditions""" return np.exp(-(self.x - 5)**2 / 0.5) def explicit_step(self, u, dt): """One step by explicit method""" u_new = u.copy() # Implementation details return u_new def compute_error(self, numerical, analytical): """Error evaluation""" return np.linalg.norm(numerical - analytical) / np.linalg.norm(analytical) def plot_comparison(self): """Comparison plot of solutions""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) for i, (ax, method) in enumerate(zip(axes.flat, ['explicit', 'implicit', 'crank-nicolson', 'spectral'])): ax.set_title(f'{method} method', fontsize=12, fontweight='bold') ax.set_xlabel('x') ax.set_ylabel('u(x,t)') ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Execution example solver = AdvancedSolver(N=200, method='implicit') solution = solver.solve(T=2.0, dt=0.01) print(f"Calculation complete: {solution.shape} time steps")

import numpy as np import matplotlib.pyplot as plt class StabilityAnalyzer: """Stability analysis tool""" def __init__(self): self.k_values = np.linspace(0, np.pi, 100) def amplification_factor(self, k, dt, dx, method='FTCS'): """Compute amplification factor""" r = dt / dx**2 if method == 'FTCS': # Forward Time Centered Space g = 1 - 4*r*np.sin(k/2)**2 elif method == 'BTCS': # Backward Time Centered Space g = 1 / (1 + 4*r*np.sin(k/2)**2) elif method == 'Crank-Nicolson': g = (1 - 2*r*np.sin(k/2)**2) / (1 + 2*r*np.sin(k/2)**2) else: g = np.ones_like(k) return g def plot_stability_regions(self): """Plot stability regions""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) methods = ['FTCS', 'BTCS', 'Crank-Nicolson', 'Upwind'] r_values = [0.1, 0.3, 0.5, 0.7] for ax, method in zip(axes.flat, methods): for r in r_values: g = self.amplification_factor(self.k_values, r, 1.0, method) ax.plot(self.k_values, np.abs(g), label=f'r={r}') ax.axhline(y=1, color='k', linestyle='--', alpha=0.5) ax.set_xlabel('Wavenumber k', fontsize=12) ax.set_ylabel('|Amplification factor|', fontsize=12) ax.set_title(f'{method} Method Stability', fontsize=12, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig def von_neumann_analysis(self, dt, dx): """von Neumann stability analysis""" r = dt / dx**2 # CFL condition check if r > 0.5: print(f"Warning: CFL condition violated (r={r:.3f} > 0.5)") return False else: print(f"Stable: r={r:.3f} ≤ 0.5") return True def convergence_test(self): """Convergence test""" dx_values = [0.1, 0.05, 0.025, 0.0125] errors = [] for dx in dx_values: # Error calculation between numerical and theoretical solutions error = dx**2 # Assuming second-order accuracy errors.append(error) # Estimate convergence order fig, ax = plt.subplots(figsize=(10, 6)) ax.loglog(dx_values, errors, 'bo-', linewidth=2, markersize=8, label='Numerical error') ax.loglog(dx_values, [dx**2 for dx in dx_values], 'r--', label='O(Δx²)') ax.set_xlabel('Grid spacing Δx', fontsize=12) ax.set_ylabel('Error', fontsize=12) ax.set_title('Convergence Test', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) return fig # Execution example analyzer = StabilityAnalyzer() fig = analyzer.plot_stability_regions() plt.show() analyzer.von_neumann_analysis(dt=0.001, dx=0.1)

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D class Solver2D: """Solver for 2D problems""" def __init__(self, Nx=50, Ny=50, Lx=1.0, Ly=1.0): self.Nx = Nx self.Ny = Ny self.Lx = Lx self.Ly = Ly self.x = np.linspace(0, Lx, Nx) self.y = np.linspace(0, Ly, Ny) self.X, self.Y = np.meshgrid(self.x, self.y) self.dx = self.x[1] - self.x[0] self.dy = self.y[1] - self.y[0] self.u = np.zeros((Nx, Ny)) def initialize_gaussian(self, x0=0.5, y0=0.5, sigma=0.1): """Gaussian initial conditions""" self.u = np.exp(-((self.X - x0)**2 + (self.Y - y0)**2) / (2*sigma**2)) def laplacian_2d(self, u): """2D Laplacian""" laplacian = np.zeros_like(u) # Interior points laplacian[1:-1, 1:-1] = ( (u[2:, 1:-1] - 2*u[1:-1, 1:-1] + u[:-2, 1:-1]) / self.dx**2 + (u[1:-1, 2:] - 2*u[1:-1, 1:-1] + u[1:-1, :-2]) / self.dy**2 ) return laplacian def step(self, dt): """Time evolution (one step)""" lap = self.laplacian_2d(self.u) self.u += dt * lap def solve(self, T=0.1, dt=0.001): """Time evolution solver""" steps = int(T / dt) for n in range(steps): self.step(dt) if n % 10 == 0: print(f"Step {n}/{steps}") return self.u def plot_solution(self): """Solution visualization""" fig = plt.figure(figsize=(16, 5)) # 2D contour plot ax1 = fig.add_subplot(131) im = ax1.contourf(self.X, self.Y, self.u, levels=20, cmap='viridis') ax1.set_xlabel('x', fontsize=12) ax1.set_ylabel('y', fontsize=12) ax1.set_title('Contour Plot', fontsize=12, fontweight='bold') plt.colorbar(im, ax=ax1) # 3D surface ax2 = fig.add_subplot(132, projection='3d') surf = ax2.plot_surface(self.X, self.Y, self.u, cmap='plasma', alpha=0.8) ax2.set_xlabel('x') ax2.set_ylabel('y') ax2.set_zlabel('u(x,y)') ax2.set_title('3D Surface', fontsize=12, fontweight='bold') # Cross section ax3 = fig.add_subplot(133) mid_y = self.Ny // 2 ax3.plot(self.x, self.u[:, mid_y], 'b-', linewidth=2, label='y=0.5 cross section') ax3.set_xlabel('x', fontsize=12) ax3.set_ylabel('u(x, y=0.5)', fontsize=12) ax3.set_title('Cross Section Profile', fontsize=12, fontweight='bold') ax3.legend() ax3.grid(True, alpha=0.3) plt.tight_layout() return fig # Execution example solver_2d = Solver2D(Nx=100, Ny=100) solver_2d.initialize_gaussian(x0=0.5, y0=0.5, sigma=0.1) solver_2d.solve(T=0.1, dt=0.0005) fig = solver_2d.plot_solution() plt.show()

import numpy as np import matplotlib.pyplot as plt class MaterialsApplication: """Application simulation to materials science""" def __init__(self, N=200, L=10.0): self.N = N self.L = L self.x = np.linspace(0, L, N) self.dx = self.x[1] - self.x[0] # Material properties self.diffusivity = 1.0 self.reaction_rate = 0.1 def reaction_diffusion(self, u, v): """Reaction-diffusion equations""" # Laplacian lap_u = np.zeros_like(u) lap_u[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / self.dx**2 lap_v = np.zeros_like(v) lap_v[1:-1] = (v[2:] - 2*v[1:-1] + v[:-2]) / self.dx**2 # Reaction term f = u * v**2 du_dt = self.diffusivity * lap_u - f dv_dt = 0.5 * self.diffusivity * lap_v + f return du_dt, dv_dt def simulate_process(self, T=50.0, dt=0.01): """Materials process simulation""" steps = int(T / dt) # Initial conditions u = np.ones(self.N) v = np.zeros(self.N) v[self.N//4:3*self.N//4] = 1.0 # Time evolution u_history = [u.copy()] v_history = [v.copy()] for n in range(steps): du, dv = self.reaction_diffusion(u, v) u += dt * du v += dt * dv # Boundary conditions u[0] = u[1] u[-1] = u[-2] v[0] = v[1] v[-1] = v[-2] if n % 100 == 0: u_history.append(u.copy()) v_history.append(v.copy()) return u_history, v_history def plot_process(self, u_history, v_history): """Process visualization""" fig, axes = plt.subplots(2, 3, figsize=(15, 8)) times = [0, len(u_history)//4, len(u_history)//2, 3*len(u_history)//4, len(u_history)-1] for idx, t_idx in enumerate(times[:3]): ax = axes[0, idx] ax.plot(self.x, u_history[t_idx], 'b-', linewidth=2, label='Component A') ax.plot(self.x, v_history[t_idx], 'r-', linewidth=2, label='Component B') ax.set_xlabel('Position x', fontsize=10) ax.set_ylabel('Concentration', fontsize=10) ax.set_title(f't = {t_idx*10:.1f}', fontsize=11, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) for idx, t_idx in enumerate(times[2:]): ax = axes[1, idx] ax.plot(self.x, u_history[t_idx], 'b-', linewidth=2, label='Component A') ax.plot(self.x, v_history[t_idx], 'r-', linewidth=2, label='Component B') ax.set_xlabel('Position x', fontsize=10) ax.set_ylabel('Concentration', fontsize=10) ax.set_title(f't = {t_idx*10:.1f}', fontsize=11, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Execution example app = MaterialsApplication(N=200, L=10.0) u_hist, v_hist = app.simulate_process(T=100.0, dt=0.01) fig = app.plot_process(u_hist, v_hist) plt.show() print("Materials process simulation complete")

import numpy as np import matplotlib.pyplot as plt import time from numba import jit class PerformanceOptimizer: """Performance optimization and benchmarking""" @staticmethod def naive_implementation(N, steps): """Implementation without optimization""" x = np.linspace(0, 10, N) u = np.exp(-x**2) start_time = time.time() for _ in range(steps): u_new = u.copy() for i in range(1, N-1): u_new[i] = 0.25 * (u[i-1] + 2*u[i] + u[i+1]) u = u_new elapsed = time.time() - start_time return u, elapsed @staticmethod def vectorized_implementation(N, steps): """Vectorized implementation""" x = np.linspace(0, 10, N) u = np.exp(-x**2) start_time = time.time() for _ in range(steps): u[1:-1] = 0.25 * (u[:-2] + 2*u[1:-1] + u[2:]) elapsed = time.time() - start_time return u, elapsed @staticmethod @jit(nopython=True) def numba_kernel(u, N, steps): """Numba JIT optimization kernel""" for _ in range(steps): u_new = u.copy() for i in range(1, N-1): u_new[i] = 0.25 * (u[i-1] + 2*u[i] + u[i+1]) u = u_new return u @staticmethod def numba_implementation(N, steps): """Numba JIT implementation""" x = np.linspace(0, 10, N) u = np.exp(-x**2) start_time = time.time() u = PerformanceOptimizer.numba_kernel(u, N, steps) elapsed = time.time() - start_time return u, elapsed def benchmark(self): """Execute benchmark""" N_values = [100, 500, 1000, 2000] steps = 1000 results = { 'naive': [], 'vectorized': [], 'numba': [] } for N in N_values: print(f"Benchmarking N={N}...") # Naive _, t_naive = self.naive_implementation(N, steps) results['naive'].append(t_naive) # Vectorized _, t_vec = self.vectorized_implementation(N, steps) results['vectorized'].append(t_vec) # Numba (warm-up) self.numba_implementation(10, 10) _, t_numba = self.numba_implementation(N, steps) results['numba'].append(t_numba) return N_values, results def plot_benchmark(self, N_values, results): """Plot benchmark results""" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # Execution time ax1.plot(N_values, results['naive'], 'o-', linewidth=2, label='Naive', markersize=8) ax1.plot(N_values, results['vectorized'], 's-', linewidth=2, label='Vectorized', markersize=8) ax1.plot(N_values, results['numba'], '^-', linewidth=2, label='Numba JIT', markersize=8) ax1.set_xlabel('Number of grid points N', fontsize=12) ax1.set_ylabel('Execution time [s]', fontsize=12) ax1.set_title('Performance Comparison', fontsize=14, fontweight='bold') ax1.legend() ax1.grid(True, alpha=0.3) # Speedup factor speedup_vec = np.array(results['naive']) / np.array(results['vectorized']) speedup_numba = np.array(results['naive']) / np.array(results['numba']) ax2.bar(np.array(N_values) - 50, speedup_vec, width=80, alpha=0.7, label='Vectorized') ax2.bar(np.array(N_values) + 50, speedup_numba, width=80, alpha=0.7, label='Numba') ax2.set_xlabel('Number of grid points N', fontsize=12) ax2.set_ylabel('Speedup factor', fontsize=12) ax2.set_title('Speedup (compared to Naive)', fontsize=14, fontweight='bold') ax2.legend() ax2.grid(True, alpha=0.3, axis='y') plt.tight_layout() return fig # Execution example optimizer = PerformanceOptimizer() N_vals, bench_results = optimizer.benchmark() fig = optimizer.plot_benchmark(N_vals, bench_results) plt.show() print("\n=== Benchmark Results ===") for N, t_naive, t_vec, t_numba in zip(N_vals, bench_results['naive'], bench_results['vectorized'], bench_results['numba']): print(f"N={N}: Naive={t_naive:.4f}s, Vectorized={t_vec:.4f}s, " f"Numba={t_numba:.4f}s")

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class ComprehensiveProject: """Comprehensive exercise project: Practical applications of master equations and stochastic processes""" def __init__(self, Nx=150, Ny=150): self.Nx = Nx self.Ny = Ny self.x = np.linspace(0, 10, Nx) self.y = np.linspace(0, 10, Ny) self.X, self.Y = np.meshgrid(self.x, self.y) self.dx = self.x[1] - self.x[0] self.dy = self.y[1] - self.y[0] # State variables self.field = np.zeros((Nx, Ny)) self.auxiliary = np.zeros((Nx, Ny)) def initialize_complex_condition(self): """Complex initial conditions""" # Superposition of multiple Gaussian distributions centers = [(2, 2), (5, 5), (8, 8), (2, 8), (8, 2)] for (cx, cy) in centers: self.field += np.exp(-((self.X - cx)**2 + (self.Y - cy)**2) / 0.3) # Add noise self.field += 0.1 * np.random.randn(self.Nx, self.Ny) def coupled_evolution(self, dt): """Coupled evolution equations""" # Laplacian computation lap_field = np.zeros_like(self.field) lap_field[1:-1, 1:-1] = ( (self.field[2:, 1:-1] - 2*self.field[1:-1, 1:-1] + self.field[:-2, 1:-1]) / self.dx**2 + (self.field[1:-1, 2:] - 2*self.field[1:-1, 1:-1] + self.field[1:-1, :-2]) / self.dy**2 ) # Nonlinear term nonlinear_term = self.field**2 - self.field**3 # Time evolution self.field += dt * (lap_field + nonlinear_term) # Boundary conditions (Neumann) self.field[0, :] = self.field[1, :] self.field[-1, :] = self.field[-2, :] self.field[:, 0] = self.field[:, 1] self.field[:, -1] = self.field[:, -2] def compute_statistics(self): """Compute statistical quantities""" mean = np.mean(self.field) std = np.std(self.field) total_energy = np.sum(self.field**2) * self.dx * self.dy return { 'mean': mean, 'std': std, 'energy': total_energy, 'min': np.min(self.field), 'max': np.max(self.field) } def run_simulation(self, T=5.0, dt=0.01): """Run simulation""" steps = int(T / dt) # Record statistical information stats_history = [] for n in range(steps): self.coupled_evolution(dt) if n % 10 == 0: stats = self.compute_statistics() stats['time'] = n * dt stats_history.append(stats) print(f"Step {n}/{steps}: Energy={stats['energy']:.4f}") return stats_history def plot_final_results(self, stats_history): """Comprehensive plot of final results""" fig = plt.figure(figsize=(16, 12)) gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3) # Final field (contour) ax1 = fig.add_subplot(gs[0, :2]) im1 = ax1.contourf(self.X, self.Y, self.field, levels=30, cmap='RdBu_r') ax1.set_xlabel('x', fontsize=12) ax1.set_ylabel('y', fontsize=12) ax1.set_title('Final Field Distribution', fontsize=14, fontweight='bold') plt.colorbar(im1, ax=ax1) # 3D surface ax2 = fig.add_subplot(gs[0, 2], projection='3d') surf = ax2.plot_surface(self.X[::3, ::3], self.Y[::3, ::3], self.field[::3, ::3], cmap='viridis') ax2.set_xlabel('x') ax2.set_ylabel('y') ax2.set_zlabel('field') ax2.set_title('3D Visualization', fontsize=12, fontweight='bold') # Energy time evolution ax3 = fig.add_subplot(gs[1, 0]) times = [s['time'] for s in stats_history] energies = [s['energy'] for s in stats_history] ax3.plot(times, energies, 'b-', linewidth=2) ax3.set_xlabel('Time t', fontsize=12) ax3.set_ylabel('Total energy', fontsize=12) ax3.set_title('Energy Conservation', fontsize=12, fontweight='bold') ax3.grid(True, alpha=0.3) # Statistical quantities time evolution ax4 = fig.add_subplot(gs[1, 1]) means = [s['mean'] for s in stats_history] stds = [s['std'] for s in stats_history] ax4.plot(times, means, 'r-', linewidth=2, label='Mean') ax4.plot(times, stds, 'g-', linewidth=2, label='Standard deviation') ax4.set_xlabel('Time t', fontsize=12) ax4.set_ylabel('Statistics', fontsize=12) ax4.set_title('Statistical Properties', fontsize=12, fontweight='bold') ax4.legend() ax4.grid(True, alpha=0.3) # Min/max values ax5 = fig.add_subplot(gs[1, 2]) mins = [s['min'] for s in stats_history] maxs = [s['max'] for s in stats_history] ax5.plot(times, mins, 'b-', linewidth=2, label='Minimum') ax5.plot(times, maxs, 'r-', linewidth=2, label='Maximum') ax5.set_xlabel('Time t', fontsize=12) ax5.set_ylabel('Value', fontsize=12) ax5.set_title('Min/Max Values', fontsize=12, fontweight='bold') ax5.legend() ax5.grid(True, alpha=0.3) # Histogram ax6 = fig.add_subplot(gs[2, 0]) ax6.hist(self.field.flatten(), bins=50, alpha=0.7, color='blue') ax6.set_xlabel('Field value', fontsize=12) ax6.set_ylabel('Frequency', fontsize=12) ax6.set_title('Value Distribution', fontsize=12, fontweight='bold') ax6.grid(True, alpha=0.3, axis='y') # Cross section profile ax7 = fig.add_subplot(gs[2, 1]) mid_y = self.Ny // 2 ax7.plot(self.x, self.field[:, mid_y], 'b-', linewidth=2) ax7.set_xlabel('x', fontsize=12) ax7.set_ylabel('field(x, y=5)', fontsize=12) ax7.set_title('Central Cross Section', fontsize=12, fontweight='bold') ax7.grid(True, alpha=0.3) # Power spectrum ax8 = fig.add_subplot(gs[2, 2]) fft_field = np.fft.fft2(self.field) power_spectrum = np.abs(fft_field)**2 power_spectrum_1d = np.mean(power_spectrum, axis=0) freq = np.fft.fftfreq(self.Nx, self.dx) ax8.loglog(freq[1:self.Nx//2], power_spectrum_1d[1:self.Nx//2], 'b-') ax8.set_xlabel('Wavenumber', fontsize=12) ax8.set_ylabel('Power', fontsize=12) ax8.set_title('Power Spectrum', fontsize=12, fontweight='bold') ax8.grid(True, alpha=0.3) return fig # Main execution project = ComprehensiveProject(Nx=150, Ny=150) project.initialize_complex_condition() print("Starting comprehensive simulation...") stats_history = project.run_simulation(T=5.0, dt=0.01) fig = project.plot_final_results(stats_history) plt.show() print("\n=== Simulation Complete ===") print(f"Final energy: {stats_history[-1]['energy']:.4f}") print(f"Final mean: {stats_history[-1]['mean']:.4f}") print(f"Final standard deviation: {stats_history[-1]['std']:.4f}")