Here’s a challenging number theory problem solved using $Python$, with a clear explanation and visualization.

The problem is related to modular arithmetic and the Chinese Remainder Theorem (CRT):

Problem

Find all integers $( x )$ such that:
$$
x \equiv 2 \pmod{3}, \quad x \equiv 3 \pmod{5}, \quad x \equiv 2 \pmod{7}
$$
Then visualize the results within a range, highlighting the congruences visually.

Explanation

This is a system of simultaneous congruences.

The Chinese Remainder Theorem guarantees a unique solution modulo $( n = 3 \cdot 5 \cdot 7 = 105 )$ since the moduli are pairwise coprime.

1. Solve for $( x )$ satisfying the given congruences.
2. Verify solutions by computing $( x \mod 3 )$, $( x \mod 5 )$, and $( x \mod 7 )$.
3. Visualize the congruences within a given range, e.g., $( 0 \leq x < 200 )$.

Python Implementation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import numpy as np
import matplotlib.pyplot as plt
from sympy.ntheory.modular import solve_congruence

# Define the congruences
congruences = [(2, 3), (3, 5), (2, 7)]

# Solve using the Chinese Remainder Theorem
solution, modulus = solve_congruence(*congruences)

print(f"Solution: x ≡ {solution} (mod {modulus})")

# Generate values of x within a range
x_range = 200  # Define the range for visualization
x_values = [solution + i * modulus for i in range((x_range // modulus) + 1)]

# Verify the solution
assert all(x % m == r for r, m in congruences for x in x_values)

# Visualization
plt.figure(figsize=(10, 6))

# Plot the values of x
plt.scatter(x_values, [0] * len(x_values), color='blue', label='Solutions (x)')

# Show congruences modulo 3, 5, and 7
moduli = [3, 5, 7]
colors = ['red', 'green', 'orange']
for mod, color in zip(moduli, colors):
    congruent_vals = [x for x in range(x_range) if x % mod == solution % mod]
    plt.scatter(congruent_vals, [mod] * len(congruent_vals), color=color, alpha=0.5, label=f'x ≡ {solution % mod} (mod {mod})')

# Annotations and labels
plt.title('Visualization of Solutions to Simultaneous Congruences', fontsize=14)
plt.xlabel('x', fontsize=12)
plt.ylabel('Modulo Conditions', fontsize=12)
plt.legend()
plt.grid(alpha=0.3)
plt.show()

Explanation of Code

1. Solve_congruence: Utilizes $SymPy$ to solve the system of congruences.
2. x_values: Computes the values satisfying the congruences within a specified range.
3. Visualization:
   • Scatter plot for $( x )$ values.
   • Highlights modulo congruence relationships $( x \mod 3 )$, $( x \mod 5 )$, $( x \mod 7 )$.

Expected Output

• Text Output:

  1	Solution: x ≡ 23 (mod 105)

  This means the smallest solution is $( x = 23 )$, and all solutions are of the form $( x = 23 + 105k )$ for integers $( k )$.

• Graphical Output:

  • Blue points indicate solutions $( x )$.
  • Colored scatter plots at $( y = 3, 5, 7 )$ show congruences modulo $3$, $5$, and $7$ respectively.

Problem: Partition Function and Average Energy in a 2-Level System

In statistical mechanics, the partition function is a fundamental quantity that describes the statistical properties of a system in thermal equilibrium.

Using $Python$, we can compute the partition function and derive physical quantities such as the average energy of a 2-level quantum system.

Problem Statement

Consider a system with two energy levels:

1. Ground state energy $E_0 = 0$,
2. Excited state energy $E_1 = \Delta E$.

The system is in thermal equilibrium at temperature $T$. Solve for:

1. The partition function $Z$:
   $$
   Z = \sum_{i} e^{-E_i / k_B T}
   $$
2. The average energy $ \langle E \rangle$:
   $$
   \langle E \rangle = \frac{\sum_{i} E_i e^{-E_i / k_B T}}{Z}
   $$

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

import numpy as np
import matplotlib.pyplot as plt

# Constants
k_B = 1.38e-23  # Boltzmann constant (J/K)
Delta_E = 1.5e-21  # Energy difference between levels (J)

# Functions for partition function and average energy
def partition_function(T):
    return 1 + np.exp(-Delta_E / (k_B * T))

def average_energy(T):
    Z = partition_function(T)
    numerator = (0 * np.exp(0) + Delta_E * np.exp(-Delta_E / (k_B * T)))
    return numerator / Z

# Temperature range
T = np.linspace(1, 1000, 500)  # Temperature in Kelvin

# Calculate partition function and average energy
Z = partition_function(T)
E_avg = average_energy(T)

# Plot results
plt.figure(figsize=(12, 6))

# Partition function
plt.subplot(1, 2, 1)
plt.plot(T, Z, label="Partition Function (Z)", color='blue')
plt.title("Partition Function vs Temperature")
plt.xlabel("Temperature (K)")
plt.ylabel("Z")
plt.grid(alpha=0.3)

# Average energy
plt.subplot(1, 2, 2)
plt.plot(T, E_avg, label="Average Energy ⟨E⟩", color='red')
plt.title("Average Energy vs Temperature")
plt.xlabel("Temperature (K)")
plt.ylabel("⟨E⟩ (J)")
plt.grid(alpha=0.3)

plt.tight_layout()
plt.show()

Explanation

1. Partition Function:

   • The partition function $Z$ is the sum of the $Boltzmann$ factors $e^{-E_i / k_B T}$.
   • For the 2-level system, $Z = 1 + e^{-\Delta E / k_B T}$.

2. Average Energy:

   • The average energy $\langle E \rangle$ is computed by weighting each energy level by its $Boltzmann$ factor, normalized by $Z$.

3. Temperature Dependence:

   • At low temperatures, most particles are in the ground state, and $\langle E \rangle \approx 0$.
   • At high temperatures, both levels are equally populated, and $\langle E \rangle \to \Delta E / 2$.

Results and Insights

• Partition Function:

  • $Z$ starts close to 1 at low temperatures, where only the ground state is occupied.
  • As temperature increases, $Z$ grows, reflecting the contribution of the excited state.

• Average Energy:

  • At low $T$, $(\langle E \rangle \to 0)$, as particles are mostly in the ground state.
  • At high $(T)$, the system reaches a thermal equilibrium where $(\langle E \rangle \to \Delta E / 2)$, consistent with equal population of both states.

This example demonstrates how statistical mechanics connects microscopic energy states with macroscopic thermodynamic quantities.

Problem: Numerical Solution of a Nonlinear Equation

In mathematical analysis, solving nonlinear equations is essential for understanding systems in engineering, physics, and economics.

Here, we use Newton’s Method to find the root of a nonlinear equation.

Problem Statement

Solve the nonlinear equation:
$$
f(x) = x^3 - 2x^2 - x + 2 = 0
$$
starting with an initial guess $(x_0 = 1.5)$.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

import numpy as np
import matplotlib.pyplot as plt

# Define the function and its derivative
def f(x):
    return x**3 - 2*x**2 - x + 2

def f_prime(x):
    return 3*x**2 - 4*x - 1

# Implement Newton's Method
def newtons_method(f, f_prime, x0, tol=1e-6, max_iter=100):
    x = x0
    history = [x]
    for _ in range(max_iter):
        x_new = x - f(x) / f_prime(x)
        history.append(x_new)
        if abs(x_new - x) < tol:
            break
        x = x_new
    return x_new, history

# Solve the equation
x0 = 1.5  # Initial guess
root, history = newtons_method(f, f_prime, x0)

# Generate points for plotting
x = np.linspace(-2, 3, 500)
y = f(x)

# Plot the function and the iterations
plt.figure(figsize=(8, 5))
plt.plot(x, y, label="f(x) = $x^3 - 2x^2 - x + 2$", color='blue')
plt.axhline(0, color='black', linewidth=0.8, linestyle='--', label="y = 0")
plt.scatter(history, [f(xi) for xi in history], color='red', zorder=5, label="Iterations")
plt.plot(history, [f(xi) for xi in history], color='orange', linestyle='--', label="Convergence Path")

# Highlight the root
plt.scatter(root, 0, color='green', s=100, label=f"Root: x = {root:.6f}")

# Add labels, legend, and grid
plt.title("Numerical Solution of $x^3 - 2x^2 - x + 2 = 0$ Using Newton's Method")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.grid(alpha=0.3)
plt.show()

# Print the result
print(f"The root of the equation is approximately: x = {root:.6f}")

Explanation of the Code

1. Newton’s Method:

   • Newton’s Method is an iterative algorithm for finding roots of nonlinear equations.
   • The iteration formula is:
     $$
     x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)}
     $$
   • The method requires an initial guess $x_0$, the function $f(x)$, and its derivative $f’(x)$.

2. Convergence Tracking:

   • The history of $(x_n)$ values is stored and plotted to show the convergence path.

3. Graphical Representation:

   • The function $f(x)$ is plotted along with the iteration points.
   • The root is highlighted on the graph.

Results and Insights

This graph illustrates the solution of the nonlinear equation $f(x) = x^3 - 2x^2 - x + 2$ using Newton’s Method.

• The blue curve represents the function $f(x)$.
• The black dashed line is $y = 0$, showing the x-axis where the root lies.
• Red dots indicate the iteration points computed during the method.
• The orange dashed line represents the convergence path of the iterations.
• The green dot marks the final root found at $x = -1.000000$, where $f(x) = 0$.

This visualization demonstrates how Newton’s Method refines the guess iteratively, converging to the root of the equation.

Problem: Solving a Differential Equation with Boundary Conditions

In applied mathematics, differential equations model physical systems such as heat transfer, population dynamics, or mechanical vibrations.

This example demonstrates solving a second-order differential equation using boundary conditions.

Problem Statement

Solve the differential equation:
$$
\frac{d^2y}{dx^2} - y = 0
$$
with the boundary conditions:
$$
y(0) = 1, \quad y(1) = 2
$$

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp

# Define the differential equation
def equation(x, y):
    return np.vstack((y[1], y[0]))  # y'[0] = y[1], y'[1] = y[0]

# Define boundary conditions
def boundary_conditions(ya, yb):
    return np.array([ya[0] - 1, yb[0] - 2])  # y(0)=1, y(1)=2

# Initial guess for the solution
x = np.linspace(0, 1, 10)
y_guess = np.zeros((2, x.size))  # Initial guess for y and y'

# Solve the boundary value problem
solution = solve_bvp(equation, boundary_conditions, x, y_guess)

# Generate a finer grid for plotting
x_plot = np.linspace(0, 1, 100)
y_plot = solution.sol(x_plot)[0]

# Plot the solution
plt.figure(figsize=(8, 5))
plt.plot(x_plot, y_plot, label="Solution $y(x)$", color='blue')
plt.scatter([0, 1], [1, 2], color='red', label="Boundary Conditions", zorder=5)
plt.title("Solution to $\\frac{d^2y}{dx^2} - y = 0$")
plt.xlabel("x")
plt.ylabel("y(x)")
plt.axhline(0, color='black', linewidth=0.5, linestyle='--')
plt.axvline(0, color='black', linewidth=0.5, linestyle='--')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

Explanation of the Code

1. Equation Definition:

   • The second-order equation is converted into a system of first-order equations:
     $$
     y_1 = y, \quad y_2 = y’
     $$
     Resulting in:
     $$
     y_1’ = y_2, \quad y_2’ = y_1
     $$

2. Boundary Conditions:

   • The boundary conditions $y(0) = 1$ and $y(1) = 2$ are applied.

3. Numerical Solution:

   • The scipy.integrate.solve_bvp function numerically solves the boundary value problem ($BVP$).

4. Visualization:

   • The computed solution is plotted as a smooth curve.
   • The boundary conditions are highlighted as red dots.

Results and Insights

The plot shows the solution $y(x)$ that satisfies the differential equation and the boundary conditions.

The method used is efficient for solving linear and nonlinear boundary value problems commonly encountered in engineering and applied sciences.

Problem: The Prisoner's Dilemma and Strategy Simulation

The Prisoner’s Dilemma is a classic example in game theory, demonstrating why two rational individuals might not cooperate, even if it appears to be in their best interest.

This example analyzes the payoff matrix and simulates multiple rounds of the dilemma using various strategies to observe long-term outcomes.

Objective

1. Define the payoff matrix for the Prisoner’s Dilemma.
2. Simulate repeated interactions using:
   • Always Defect
   • Always Cooperate
   • Tit-for-Tat (mimic the opponent’s previous move).
3. Visualize the accumulated payoffs over multiple rounds for each strategy.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

import numpy as np
import matplotlib.pyplot as plt

# Define the payoff matrix (Player 1, Player 2)
# Rows: Player 1's actions, Columns: Player 2's actions
# Actions: 0 = Cooperate, 1 = Defect
payoff_matrix = np.array([
    [[3, 3], [0, 5]],  # Player 1 cooperates
    [[5, 0], [1, 1]]   # Player 1 defects
])

# Define strategies
def always_cooperate(history, my_history):
    return 0  # Always cooperate

def always_defect(history, my_history):
    return 1  # Always defect

def tit_for_tat(history, my_history):
    if not history:  # Cooperate in the first round
        return 0
    return history[-1]  # Mimic the opponent's last move

# Simulate the game
def simulate_game(strategy1, strategy2, rounds=10):
    history1, history2 = [], []
    score1, score2 = 0, 0

    for _ in range(rounds):
        move1 = strategy1(history2, history1)
        move2 = strategy2(history1, history2)
        
        payoff1, payoff2 = payoff_matrix[move1, move2]
        score1 += payoff1
        score2 += payoff2
        
        history1.append(move1)
        history2.append(move2)
    
    return score1, score2, history1, history2

# Run simulations
rounds = 20
strategies = [always_cooperate, always_defect, tit_for_tat]
results = {}

for i, strat1 in enumerate(strategies):
    for j, strat2 in enumerate(strategies):
        if i <= j:  # Avoid duplicating reverse matches
            label = f"{strat1.__name__} vs {strat2.__name__}"
            score1, score2, hist1, hist2 = simulate_game(strat1, strat2, rounds)
            results[label] = (score1, score2, hist1, hist2)

# Visualize results
fig, ax = plt.subplots(figsize=(10, 6))
labels = list(results.keys())
player1_scores = [results[label][0] for label in labels]
player2_scores = [results[label][1] for label in labels]

bar_width = 0.35
x_indices = np.arange(len(labels))

# Plot bars for Player 1 and Player 2
ax.bar(x_indices - bar_width/2, player1_scores, bar_width, label='Player 1')
ax.bar(x_indices + bar_width/2, player2_scores, bar_width, label='Player 2')

# Add labels and title
ax.set_xlabel("Strategy Pair", fontsize=12)
ax.set_ylabel("Accumulated Payoff", fontsize=12)
ax.set_title("Prisoner's Dilemma: Payoffs for Different Strategies", fontsize=14)
ax.set_xticks(x_indices)
ax.set_xticklabels(labels, rotation=45, ha='right', fontsize=10)
ax.legend()

# Show grid for clarity
ax.grid(axis='y', linestyle='--', alpha=0.7)

plt.tight_layout()
plt.show()

Explanation of the Code

1. Payoff Matrix:

   • Each cell in the matrix represents the payoff for Player $1$ and Player $2$, given their choices (cooperate/defect).
   • Cooperation leads to mutual benefit, while defection may yield higher short-term rewards.

2. Strategies:

   • Always Cooperate: Always chooses to cooperate.
   • Always Defect: Always chooses to defect.
   • Tit-for-Tat: Mimics the opponent’s previous move, encouraging reciprocity.

3. Simulation:

   • The game is repeated for a fixed number of rounds.
   • Accumulated scores and move histories are recorded for both players.

4. Visualization:

   • Bar plots display the total payoffs for Player $1$ and Player $2$ across all strategy pairs.
   • The results highlight how different strategies perform against each other.

Results and Insights

1. Tit-for-Tat Performance:

   • Tit-for-Tat typically performs well, fostering cooperation when paired with cooperative strategies.

2. Always Defect:

   • While it exploits cooperative players, it leads to mutual punishment when paired with itself or Tit-for-Tat.

3. Always Cooperate:

   • Vulnerable to exploitation but performs well in fully cooperative scenarios.

This simulation illustrates the dynamics of strategic decision-making in repeated games, revealing insights into the stability and outcomes of different strategies.

Problem: Numerical Integration Using the Trapezoidal Rule

Numerical integration is a key topic in computational mathematics.

In this example, we will approximate the integral of a given function using the trapezoidal rule and compare the result to the exact value.

Objective

1. Integrate $ f(x) = \sin(x) $ over the interval $([0, \pi])$ numerically using the trapezoidal rule.
2. Visualize the approximation and compare it with the actual curve of the function.

Mathematical Background

The integral of $ f(x) = \sin(x) $ over $([0, \pi])$ is:
$$
\int_0^\pi \sin(x) , dx = 2
$$

The trapezoidal rule approximates the integral by dividing the interval into $( n ) $subintervals and approximating $ f(x) $ as a straight line over each subinterval.
The formula is:
$$
\int_a^b f(x) , dx \approx \frac{h}{2} \left[ f(x_0) + 2 \sum_{i=1}^{n-1} f(x_i) + f(x_n) \right]
$$
where $( h = \frac{b-a}{n} )$ is the width of each subinterval.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

import numpy as np
import matplotlib.pyplot as plt

# Function to integrate
def f(x):
    return np.sin(x)

# Trapezoidal rule implementation
def trapezoidal_rule(func, a, b, n):
    x = np.linspace(a, b, n + 1)  # Subinterval points
    y = func(x)  # Function values at those points
    h = (b - a) / n  # Width of each subinterval
    integral = (h / 2) * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])
    return integral, x, y

# Parameters
a, b = 0, np.pi  # Integration limits
n = 10  # Number of subintervals

# Calculate numerical integral
numerical_integral, x_points, y_points = trapezoidal_rule(f, a, b, n)

# Exact integral
exact_integral = 2  # From analytical solution

# Plotting
x_full = np.linspace(a, b, 1000)
y_full = f(x_full)

plt.figure(figsize=(10, 6))
plt.plot(x_full, y_full, label='f(x) = sin(x)', color='blue')
plt.fill_between(x_full, 0, y_full, color='blue', alpha=0.1, label='Exact Area')
plt.plot(x_points, y_points, 'o', color='red', label='Trapezoid Nodes')
for i in range(len(x_points) - 1):
    plt.plot([x_points[i], x_points[i + 1]], [y_points[i], y_points[i + 1]], color='orange')

plt.title("Numerical Integration Using the Trapezoidal Rule")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.7)
plt.show()

# Print results
print(f"Numerical Integral (Trapezoidal Rule): {numerical_integral:.4f}")
print(f"Exact Integral: {exact_integral}")
print(f"Error: {abs(numerical_integral - exact_integral):.4f}")

Explanation of the Code

1. Function and Interval:
   • The function $( f(x) = \sin(x) )$ is defined along with the integration interval $([0, \pi])$.
2. Trapezoidal Rule:
   • The interval is divided into $( n )$ subintervals, and the trapezoidal rule is applied to approximate the integral.
3. Exact Solution:
   • The exact integral is $( 2 )$, used here for error comparison.
4. Visualization:
   • The function is plotted alongside the trapezoids used for approximation.
   • The filled area under $ f(x) $ represents the exact integral.

Results and Insights

1. Numerical Approximation:
   • The trapezoidal rule provides an approximate integral value close to the exact solution.
     For $( n = 10 )$, the numerical result is accurate within a small error margin.
2. Graph:
   • The blue curve shows $ f(x) = \sin(x) $.
   • Red points indicate the nodes where the function is evaluated.
   • Orange lines represent the trapezoids approximating the area under the curve.
3. Error Analysis:
   • Increasing the number of subintervals $( n )$ reduces the approximation error, demonstrating the method’s convergence.

Applications in Computational Mathematics

1. Integration of Complex Functions:
   • Useful when analytical integration is impossible.
2. Engineering and Physics:
   • Applied to problems involving areas, volumes, and physical quantities.
3. Data Analysis:
   • Approximation of cumulative distributions or areas under discrete data points.

This example highlights how numerical integration works and provides a visual understanding of the trapezoidal rule.

Problem: Predicting House Prices Using Linear Regression

In this example, we analyze a dataset containing house features (e.g., size, number of rooms) to predict house prices using linear regression, a core technique in data science.

Objective

1. Train a linear regression model to predict house prices based on square footage.
2. Visualize the results, including the data points and the regression line.

Dataset

We will simulate a dataset of house sizes (in square feet) and corresponding prices (in ＄1000s).

This ensures reproducibility and simplicity.

Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

# Simulated dataset
np.random.seed(42)
sizes = np.random.uniform(500, 3500, 50)  # House sizes (square feet)
prices = sizes * 0.15 + np.random.normal(0, 50, 50)  # Prices with some noise

# Reshape data for sklearn
X = sizes.reshape(-1, 1)  # Predictor (feature)
y = prices  # Target (price)

# Train a linear regression model
model = LinearRegression()
model.fit(X, y)

# Make predictions
predicted_prices = model.predict(X)

# Evaluate the model
mse = mean_squared_error(y, predicted_prices)
print(f"Mean Squared Error: {mse:.2f}")

# Plot the results
plt.figure(figsize=(10, 6))
plt.scatter(sizes, prices, color="blue", label="Actual Prices", alpha=0.7)
plt.plot(sizes, predicted_prices, color="red", label="Regression Line", linewidth=2)
plt.title("House Prices vs. Size")
plt.xlabel("Size (Square Feet)")
plt.ylabel("Price (in $1000s)")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.7)
plt.show()

Explanation of the Code

1. Simulated Data:
   • House sizes are sampled randomly between $500$ and $3500$ square feet.
   • Prices are calculated using a simple linear relationship $( \text{Price} = 0.15 \times \text{Size} + \text{Noise} )$, where the noise simulates real-world variability.
2. Model Training:
   • A LinearRegression model from sklearn is trained using house sizes as the feature and prices as the target.
3. Prediction:
   • The trained model predicts house prices, which are plotted alongside the actual prices.
4. Model Evaluation:
   • The Mean Squared Error (MSE) quantifies the model’s accuracy.

Visualization and Insights

1. Scatter Plot:
   • Blue dots represent the actual prices for different house sizes.
2. Regression Line:
   • The red line shows the best-fit line learned by the linear regression model.
3. Model Performance:
   • The MSE indicates the average squared difference between the predicted and actual prices.

Applications in Data Science

1. Regression Analysis:
   • This approach is foundational in predictive modeling for understanding relationships between variables.
2. Real Estate Analytics:
   • Real-world data would involve more features like location, number of bedrooms, and year built.
3. Scalability:
   • This method easily extends to multivariate regression with multiple predictors.

This example demonstrates the workflow of data science: data preprocessing, modeling, and visualization, providing insights into how features like size influence house prices.

Here’s a fascinating example from statistical mechanics: the Maxwell-Boltzmann distribution of particle speeds in an ideal gas.

This distribution describes the probable speeds of particles at a given temperature, which is central to understanding gas behavior in thermodynamics.

Problem: Plotting the Maxwell-Boltzmann Speed Distribution

The Maxwell-Boltzmann distribution is a probability distribution for the speed $( v )$ of particles in an ideal gas at a certain temperature $( T )$. Given:

• Temperature $( T )$ of the gas,
• Particle mass $( m )$,

we can determine the likelihood of finding particles moving at different speeds.

The probability density function of the Maxwell-Boltzmann speed distribution is given by:
$$
f(v) = 4\pi \left( \frac{m}{2 \pi k_B T} \right)^{3/2} v^2 \exp \left( -\frac{mv^2}{2k_B T} \right)
$$

• $( v )$ is the particle speed,
• $( m )$ is the mass of each particle,
• $( T )$ is the absolute temperature,
• $( k_B )$ is the Boltzmann constant.

Goal

1. Calculate the Maxwell-Boltzmann distribution for speeds.
2. Plot the distribution to observe how it varies with temperature.

Python Code

Below is the $Python$ code to compute and plot the Maxwell-Boltzmann distribution for a given temperature and particle mass:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

import numpy as np
import matplotlib.pyplot as plt

# Constants
k_B = 1.38e-23  # Boltzmann constant in J/K
m = 4.65e-26    # Mass of a nitrogen molecule in kg (example particle)

# Temperature
T = 300  # Temperature in Kelvin

# Maxwell-Boltzmann speed distribution function
def maxwell_boltzmann(v, T, m):
    coefficient = 4 * np.pi * (m / (2 * np.pi * k_B * T)) ** (3/2)
    return coefficient * v**2 * np.exp(-m * v**2 / (2 * k_B * T))

# Speed range
v = np.linspace(0, 2000, 500)

# Calculate distribution
f_v = maxwell_boltzmann(v, T, m)

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(v, f_v, color="blue", label=f"T = {T} K")
plt.xlabel("Speed (m/s)")
plt.ylabel("Probability Density f(v)")
plt.title("Maxwell-Boltzmann Speed Distribution")
plt.legend()
plt.grid(True, linestyle="--", alpha=0.7)

# Display the plot
plt.show()

Explanation of the Code

1. Function Definition:
   • We define a function maxwell_boltzmann(v, T, m) to calculate the probability density $ f(v) $ at different speeds $( v )$ for a given temperature $( T )$ and particle mass $( m )$.
2. Speed Range:
   • We create an array v representing a range of speeds up to around $2000$ m/s, which covers the typical speeds of nitrogen molecules at room temperature.
3. Distribution Calculation:
   • Using the maxwell_boltzmann function, we compute $ f(v) $ for each speed in v.
4. Plotting:
   • We plot $( v )$ on the $x$-axis and $ f(v) $ on the $y$-axis, giving a clear visualization of the speed distribution.

Visualization and Interpretation

The graph shows:

• Peak: The curve peaks at a specific speed, which represents the most probable speed of particles at this temperature.
• Spread: The distribution has a long tail, indicating some particles have higher speeds, though with decreasing probability.
• Temperature Dependence: As temperature increases (not shown here but if varied), the peak shifts to higher speeds, indicating faster-moving particles.

Insights from the Maxwell-Boltzmann Distribution

1. Temperature Effect: Higher temperatures lead to broader, more spread-out distributions, as particles move faster on average.
2. Statistical Behavior: The distribution embodies the statistical nature of molecular motion, providing a basis for understanding macroscopic gas properties, like pressure and temperature.
3. Applications: Maxwell-Boltzmann statistics are crucial in fields like chemistry, thermodynamics, and kinetic theory, helping explain diffusion rates, reaction kinetics, and energy transfer.

This visualization brings out the fundamental principles of statistical mechanics and demonstrates how particle motion distributions are governed by temperature.

Here’s a visually engaging number theory problem that explores Fermat’s Little Theorem through cyclic patterns in modular exponentiation.

Problem: Visualizing Cyclic Patterns in Modular Exponentiation

Fermat’s Little Theorem states that if $( p )$ is a prime number and $( a )$ is an integer not divisible by $( p )$, then:
$$
a^{p-1} \equiv 1 \pmod{p}
$$
This theorem implies that for each $( a )$, the sequence of powers $( a^1, a^2, \dots, a^{p-1} )$ modulo $( p )$ will eventually repeat, forming a cyclic pattern.

In this example, we will:

1. Choose a modulus $( p )$, a prime number (e.g., $23$).
2. Compute $( a^k \mod p )$ for each integer $( a )$ from $1$ to $( p-1 )$ and for $( k )$ from $1$ to $( p-1 )$.
3. Plot these powers to observe any patterns and cycles.

Approach

We’ll calculate $( a^k \mod p )$ for each value of $( a )$ from $1$ to $( p-1 )$ and plot the results as points, where the x-axis represents the base $( a )$, the y-axis represents the result of $( a^k \mod p )$, and different colors show different powers $( k )$.

Python Code

Here’s the $Python$ code to find and plot the cyclic patterns:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import matplotlib.pyplot as plt
import numpy as np

# Define a prime modulus
p = 23

# Generate data for powers
data = []
for a in range(1, p):
    for k in range(1, p):
        value = (a ** k) % p
        data.append((a, value, k))

# Separate values for plotting
x_values = [item[0] for item in data]
y_values = [item[1] for item in data]
colors = [item[2] for item in data]  # Different colors for each power k

# Plotting
plt.figure(figsize=(12, 8))
scatter = plt.scatter(x_values, y_values, c=colors, cmap="viridis", s=50, alpha=0.7, edgecolor="black")
plt.colorbar(scatter, label="Exponent k")
plt.xlabel("Base (a)")
plt.ylabel("Result of a^k mod 23")
plt.title("Cyclic Patterns in Modular Exponentiation (mod 23)")
plt.grid(True, linestyle="--", alpha=0.6)

# Show the plot
plt.show()

Explanation of the Code

1. Modular Exponentiation Calculation:
   • For each integer $( a )$ from $1$ to $( p-1 )$, we compute $( a^k \mod p )$ for $( k = 1 )$ to $( p-1 )$.
   • These results reveal the pattern of powers under modulo $( p )$, showing the periodic behavior predicted by Fermat’s Little Theorem.
2. Plotting:
   • Each point $(a, a^k \mod p)$ represents the result of raising $( a )$ to power $( k )$ modulo $( p )$.
   • The color of each point indicates the power $( k )$, which helps in visualizing different cycles.

Visualization

The scatter plot shows:

• Cyclic Behavior: Each base $( a )$ follows a distinct cycle when raised to powers $( k )$ modulo $( p )$.
• Color Coding: Different colors represent different exponents $( k )$, making it easy to observe how the results cycle back, eventually repeating after certain intervals.

Interpretation

1. Cyclic Patterns: This visualization shows how powers of numbers “wrap around” in modular arithmetic, creating cycles based on Fermat’s Little Theorem.
2. Applications: Modular exponentiation is central to cryptography, particularly in $RSA$ and Diffie-Hellman key exchange, where large powers modulo a prime number are computed.
3. Visual Insight: Observing the cyclic patterns provides an intuitive grasp of Fermat’s Little Theorem and the repeating nature of modular exponentiation.

This graph effectively brings out the cyclic beauty in modular exponentiation and demonstrates the predictability and structure found in number theory.

Let’s look at a problem in probability theory related to binomial distributions, which are used to model the probability of a certain number of successes in a series of independent trials.

Problem: Binomial Probability of Flipping a Coin

Suppose you flip a fair coin ($50$% chance of heads and $50$% chance of tails) $20$ times.
We want to know:

1. The probability of getting exactly $(k)$ heads for various values of $(k)$ ($0$ to $20$).
2. How the probabilities are distributed across different numbers of heads.

Solution Outline

1. We’ll use the binomial distribution. The probability of getting exactly $(k)$ heads in $(n)$ flips is given by:
   $$
   P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}
   $$

   • $(n)$ = number of trials ($20$),
   • $(p)$ = probability of success in a single trial ($0.5$ for heads),
   • $(k)$ = number of successes (number of heads).

2. We’ll calculate the probabilities for each value of $(k)$ from $0$ to $20$ and visualize the distribution of these probabilities.

Python Code

Here’s the $Python$ code to calculate and plot the binomial distribution for this problem:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom

# Define parameters
n = 20  # number of trials (coin flips)
p = 0.5  # probability of getting heads in a single flip

# Values of k (number of heads) from 0 to 20
k_values = np.arange(0, n + 1)

# Calculate binomial probabilities for each k
probabilities = binom.pmf(k_values, n, p)

# Plotting the binomial distribution
plt.figure(figsize=(10, 6))
plt.bar(k_values, probabilities, color="skyblue", edgecolor="black")

# Add labels and title
plt.xlabel("Number of Heads (k)")
plt.ylabel("Probability")
plt.title("Binomial Distribution: Probability of Getting k Heads in 20 Coin Flips")
plt.grid(axis="y", linestyle="--", alpha=0.7)

# Show the plot
plt.show()

Explanation of the Code

1. Binomial Probability Calculation:
   • We use binom.pmf(k, n, p) from scipy.stats to calculate the probability mass function (PMF) of a binomial distribution, which gives the probability of getting exactly $( k )$ heads in $( n )$ trials with probability $( p )$ of heads.
2. Plotting:
   • We plot the probabilities for each possible value of $( k )$ (from $0$ to $20$ heads).
   • A bar plot is used to visualize how probabilities are distributed across different numbers of heads.

Visualization

The bar chart shows:

• The probability of each possible outcome (number of heads) when flipping a coin $20$ times.
• The distribution is centered around $10$ heads, since we expect about half of the flips to result in heads due to the $50$% probability of success per flip.

Interpretation

1. Peak at 10 Heads:
   The distribution peaks at $10$, which is the expected number of heads when flipping a coin $20$ times with a $50$% chance for heads each time.
2. Symmetry:
   The binomial distribution here is symmetric around $10$ due to the equal probability of heads and tails.
3. Applications:
   Binomial distributions model scenarios with a fixed number of independent trials and a constant probability of success in each trial, such as quality control testing, marketing response rates, and predicting outcomes in sports or games of chance.