Option Pricing Analysis

Implementing the Black-Scholes Model with Python

I’ll provide you with a financial mathematics example and solve it using $Python$, including proper mathematical notation and visualization.

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from matplotlib import cm
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D

# Black-Scholes Option Pricing Model
def black_scholes(S, K, T, r, sigma, option='call'):
    """
    Calculate Black-Scholes option price for a call or put
    
    Parameters:
    S: Current stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        price = S * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        price = K * np.exp(-r * T) * stats.norm.cdf(-d2) - S * stats.norm.cdf(-d1)
        
    return price

# Example 1: Visualize call option price as a function of stock price and time to maturity
# Parameters
S_range = np.linspace(80, 120, 50)  # Stock price range
T_range = np.linspace(0.1, 1, 50)   # Time to maturity range
K = 100                             # Strike price
r = 0.05                            # Risk-free rate
sigma = 0.2                         # Volatility

# Create meshgrid
S_mesh, T_mesh = np.meshgrid(S_range, T_range)
option_price = np.zeros_like(S_mesh)

# Calculate option price for each combination
for i in range(len(T_range)):
    for j in range(len(S_range)):
        option_price[i, j] = black_scholes(S_mesh[i, j], K, T_mesh[i, j], r, sigma, 'call')

# Plot the 3D surface
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(S_mesh, T_mesh, option_price, cmap=cm.coolwarm, alpha=0.8)

# Add labels
ax.set_xlabel('Stock Price ($)')
ax.set_ylabel('Time to Maturity (years)')
ax.set_zlabel('Call Option Price ($)')
ax.set_title('Black-Scholes Call Option Price')

# Add a color bar
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

plt.savefig('black_scholes_3d.png')
plt.show()

# Example 2: Visualize the Greeks (Delta, Gamma, Theta)
def bs_delta(S, K, T, r, sigma, option='call'):
    """Calculate option Delta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    
    if option == 'call':
        return stats.norm.cdf(d1)
    else:
        return stats.norm.cdf(d1) - 1

def bs_gamma(S, K, T, r, sigma):
    """Calculate option Gamma (same for calls and puts)"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    return stats.norm.pdf(d1) / (S * sigma * np.sqrt(T))

def bs_theta(S, K, T, r, sigma, option='call'):
    """Calculate option Theta"""
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    if option == 'call':
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * stats.norm.cdf(d2)
    else:
        theta = -S * stats.norm.pdf(d1) * sigma / (2 * np.sqrt(T)) + r * K * np.exp(-r * T) * stats.norm.cdf(-d2)
        
    return theta / 365  # Convert to daily

# Calculate Greeks for various stock prices
S_range = np.linspace(80, 120, 100)
T = 0.5  # 6 months

delta_call = [bs_delta(S, K, T, r, sigma, 'call') for S in S_range]
delta_put = [bs_delta(S, K, T, r, sigma, 'put') for S in S_range]
gamma = [bs_gamma(S, K, T, r, sigma) for S in S_range]
theta_call = [bs_theta(S, K, T, r, sigma, 'call') for S in S_range]
theta_put = [bs_theta(S, K, T, r, sigma, 'put') for S in S_range]

# Plot the Greeks
fig, axs = plt.subplots(3, 1, figsize=(12, 15))

# Delta
axs[0].plot(S_range, delta_call, 'b-', label='Call Delta')
axs[0].plot(S_range, delta_put, 'r-', label='Put Delta')
axs[0].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[0].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[0].set_xlabel('Stock Price ($)')
axs[0].set_ylabel('Delta')
axs[0].set_title('Option Delta vs Stock Price')
axs[0].grid(True)
axs[0].legend()

# Gamma
axs[1].plot(S_range, gamma, 'g-')
axs[1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[1].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[1].set_xlabel('Stock Price ($)')
axs[1].set_ylabel('Gamma')
axs[1].set_title('Option Gamma vs Stock Price')
axs[1].grid(True)

# Theta
axs[2].plot(S_range, theta_call, 'b-', label='Call Theta')
axs[2].plot(S_range, theta_put, 'r-', label='Put Theta')
axs[2].axhline(y=0, color='k', linestyle='-', alpha=0.3)
axs[2].axvline(x=K, color='k', linestyle='--', alpha=0.3)
axs[2].set_xlabel('Stock Price ($)')
axs[2].set_ylabel('Theta ($ per day)')
axs[2].set_title('Option Theta vs Stock Price')
axs[2].grid(True)
axs[2].legend()

plt.tight_layout()
plt.savefig('option_greeks.png')
plt.show()

# Example 3: Monte Carlo simulation for option pricing
def monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations=10000, option='call'):
    """
    Price options using Monte Carlo simulation
    
    Parameters:
    S0: Initial stock price
    K: Strike price
    T: Time to maturity (in years)
    r: Risk-free interest rate
    sigma: Volatility
    num_simulations: Number of simulation paths
    option: 'call' or 'put'
    
    Returns:
    Option price
    """
    # Generate random simulations of stock price at expiration
    z = np.random.standard_normal(num_simulations)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)
    
    # Calculate payoff at expiration
    if option == 'call':
        payoff = np.maximum(ST - K, 0)
    else:
        payoff = np.maximum(K - ST, 0)
    
    # Calculate present value
    option_price = np.exp(-r * T) * np.mean(payoff)
    
    return option_price, ST

# Run Monte Carlo simulation
np.random.seed(42)  # For reproducibility
S0 = 100
num_simulations = 10000

mc_price, ST = monte_carlo_option_pricing(S0, K, T, r, sigma, num_simulations, 'call')
bs_price = black_scholes(S0, K, T, r, sigma, 'call')

# Create histogram of simulated stock prices at expiration
plt.figure(figsize=(12, 6))
plt.hist(ST, bins=50, alpha=0.5, color='blue', density=True)

# Add a normal distribution curve for comparison
x = np.linspace(min(ST), max(ST), 100)
mu = S0 * np.exp(r * T)  # Expected value under risk-neutral measure
std = S0 * np.exp(r * T) * np.sqrt(np.exp(sigma**2 * T) - 1)
plt.plot(x, stats.norm.pdf(x, mu, std), 'r-', linewidth=2)

plt.axvline(x=K, color='k', linestyle='--', label=f'Strike Price (K={K})')
plt.axvline(x=S0, color='g', linestyle='-', label=f'Initial Price (S0={S0})')

plt.title(f'Monte Carlo Simulation of Stock Prices at Expiration (T={T} years)')
plt.xlabel('Stock Price at Expiration ($)')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.savefig('monte_carlo_simulation.png')
plt.show()

# Compare Monte Carlo with Black-Scholes
print(f"Monte Carlo Call Option Price: ${mc_price:.4f}")
print(f"Black-Scholes Call Option Price: ${bs_price:.4f}")
print(f"Difference: ${abs(mc_price - bs_price):.4f}")

Black-Scholes Option Pricing Model in Financial Mathematics

In this example, I’ll solve and analyze a classic financial mathematics problem: option pricing using the Black-Scholes model.

The Black-Scholes model is fundamental in financial mathematics for pricing European-style options.

Mathematical Background

The Black-Scholes formula for European call and put options is given by:

For a call option:
$$C(S,t) = S \cdot N(d_1) - K e^{-r(T-t)} \cdot N(d_2)$$

For a put option:
$$P(S,t) = K e^{-r(T-t)} \cdot N(-d_2) - S \cdot N(-d_1)$$

Where:
$$d_1 = \frac{\ln(S/K) + (r + \sigma^2/2)(T-t)}{\sigma\sqrt{T-t}}$$
$$d_2 = d_1 - \sigma\sqrt{T-t}$$

And:

  • $S$ is the current stock price
  • $K$ is the strike price
  • $r$ is the risk-free interest rate
  • $T-t$ is the time to maturity
  • $\sigma$ is the volatility of the stock
  • $N(·)$ is the cumulative distribution function of the standard normal distribution

Code Explanation

The code implements the Black-Scholes model and provides three detailed examples:

  1. 3D Visualization of Option Prices: Shows how call option prices change with stock price and time to maturity
  2. Option Greeks Analysis: Calculates and visualizes Delta, Gamma, and Theta
  3. Monte Carlo Simulation: Prices options using simulation and compares with the analytical solution

Key Functions:

  1. black_scholes(): Implements the analytical Black-Scholes formula
  2. bs_delta(), bs_gamma(), bs_theta(): Calculate option Greeks
  3. monte_carlo_option_pricing(): Simulates stock price paths using geometric Brownian motion

Visualization Analysis

1. 3D Surface Plot of Call Option Prices

The 3D plot shows how the call option price (z-axis) varies with:

  • Stock price (x-axis): As stock price increases, call option value increases
  • Time to maturity (y-axis): The time value of the option is visible

This visualization helps understand the non-linear relationship between these variables.

2. Option Greeks Visualization

The Greeks represent sensitivities of option prices to various factors:

  • Delta: Measures the rate of change of option price with respect to changes in the underlying asset’s price

    • Call delta ranges from $0$ to $1$
    • Put delta ranges from $-1$ to $0$
    • At-the-money options have delta near $0.5$ (calls) or $-0.5$ (puts)

  • Gamma: Measures the rate of change of delta with respect to changes in the underlying price

    • Highest at-the-money and decreases as the option moves in or out of the money
    • Same for both calls and puts

  • Theta: Measures the rate of change of option price with respect to the passage of time (time decay)

    • Generally negative for both calls and puts (options lose value as time passes)
    • Most significant for at-the-money options

3. Monte Carlo Simulation

Monte Carlo Call Option Price: $6.8874
Black-Scholes Call Option Price: $6.8887
Difference: $0.0013

The histogram shows the distribution of simulated stock prices at expiration:

  • The red curve represents the theoretical lognormal distribution
  • The vertical lines mark the initial stock price and strike price
  • The simulation validates the Black-Scholes model by producing a similar price

When we compare the Monte Carlo estimate with the analytical Black-Scholes solution, they should be very close, with differences attributable to simulation error.

Conclusion

This example demonstrates several fundamental concepts in financial mathematics:

  1. Analytical solution of the Black-Scholes equation
  2. Option Greeks for risk management
  3. Monte Carlo methods for pricing financial instruments

The visualizations help develop intuition about option pricing behavior and the relationships between different variables in the model.

Visualizing the Prime Number Theorem with Python：A Comparison of \( \pi(x) \) and \( \frac{x}{\log x} \)

Here’s an example from Analytic Number Theory with a complete solution and visualization in $Python$.

📚 Problem: Prime Number Theorem Approximation

The Prime Number Theorem (PNT) states that:

$$
\pi(x) \sim \frac{x}{\log x}
$$

Where:

  • $( \pi(x) )$ is the number of primes less than or equal to $( x )$.
  • $( \frac{x}{\log x} )$ approximates $( \pi(x) )$ for large values of $( x )$.

🎯 Goal:

  1. Calculate $( \pi(x) )$ for values up to $( x = 100,000 )$.
  2. Compare it with the approximation $( \frac{x}{\log x} )$.
  3. Visualize the results using a graph.

📝 Step 1: Define and Calculate Functions

  • pi(x): Count primes less than or equal to $( x )$.
  • Approximation: Use $( \frac{x}{\log x} )$.

🧩 Step 2: Python Code 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
# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
from sympy import primepi
from math import log

# Define the range and step size
x_values = np.linspace(10, 100000, 200)  # Range from 10 to 100,000 with 200 points

# Calculate prime counting function pi(x)
pi_values = [primepi(int(x)) for x in x_values]

# Calculate the Prime Number Theorem (PNT) approximation
pnt_approx = [x / log(x) for x in x_values]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x_values, pi_values, label=r'$\pi(x)$ (Prime Count)', color='blue')
plt.plot(x_values, pnt_approx, label=r'$\frac{x}{\log x}$ (PNT Approximation)', linestyle='--', color='red')

# Add title and labels
plt.title("Prime Number Theorem: $\pi(x)$ vs. $\dfrac{x}{\log x}$", fontsize=14)
plt.xlabel("x", fontsize=12)
plt.ylabel("Value", fontsize=12)
plt.legend()
plt.grid(True)

# Show the plot
plt.show()

📊 Step 3: Explanation and Visualization

✅ Explanation:

  1. pi(x): Using primepi() from sympy to compute the exact number of primes up to $( x )$.
  2. PNT Approximation: Using the formula $( \frac{x}{\log x} )$.
  3. Graph:
    • Blue line: Exact prime counting function.
    • Red dashed line: PNT approximation.

📈 Step 4: Interpretation of Results

  • The graph shows that the approximation $( \frac{x}{\log x} )$ becomes increasingly accurate as $( x )$ increases.
  • For smaller values of $( x )$, there is a noticeable gap, but for large $( x )$, the approximation closely tracks $( \pi(x) )$.

📢 Insights:

  • The Prime Number Theorem provides a remarkably simple yet powerful approximation for the distribution of prime numbers.
  • This visualization helps highlight how the error margin decreases with larger $( x )$.

Visualization and Analysis of Matrix Operations

Linear Transformations in Python

I’ll provide an example problem about linear transformations, solve it using $Python$, and visualize the results with graphs.

Here’s a comprehensive example:

Linear Transformation Example and Solution with Python

Let’s consider a linear transformation $T: \mathbb{R}^2 \rightarrow \mathbb{R}^2$ defined by the matrix:

$$T = \begin{pmatrix} 2 & 1 \ 1 & 3 \end{pmatrix}$$

We’ll solve the following problems:

  1. Find the transformed coordinates of several points
  2. Visualize how this transformation affects a unit square
  3. Compute the eigenvalues and eigenvectors of the transformation
  4. Visualize the effect of the transformation on the eigenvectors
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d import Axes3D

# Define the linear transformation matrix
T = np.array([[2, 1], 
              [1, 3]])

print("Linear Transformation Matrix T:")
print(T)

# Part 1: Transform some example points
points = np.array([[1, 0], [0, 1], [1, 1], [-1, 2]])
transformed_points = np.dot(points, T)

print("\nOriginal Points:")
for p in points:
    print(f"[{p[0]}, {p[1]}]")
    
print("\nTransformed Points:")
for p in transformed_points:
    print(f"[{p[0]}, {p[1]}]")

# Part 2: Visualize the effect on a unit square
def plot_transformation(T):
    # Create a figure with two subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Define the unit square
    square = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
    
    # Transform the square
    transformed_square = np.dot(square, T)
    
    # Plot the original square
    ax1.plot(square[:, 0], square[:, 1], 'b-')
    ax1.fill(square[:, 0], square[:, 1], 'lightblue', alpha=0.5)
    ax1.set_xlim(-0.5, 3.5)
    ax1.set_ylim(-0.5, 3.5)
    ax1.grid(True)
    ax1.set_title('Original Unit Square')
    ax1.set_aspect('equal')
    
    # Plot the transformed square
    ax2.plot(transformed_square[:, 0], transformed_square[:, 1], 'r-')
    ax2.fill(transformed_square[:, 0], transformed_square[:, 1], 'lightcoral', alpha=0.5)
    ax2.set_xlim(-0.5, 3.5)
    ax2.set_ylim(-0.5, 3.5)
    ax2.grid(True)
    ax2.set_title('Transformed Square')
    ax2.set_aspect('equal')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation
plot_transformation(T)

# Part 3: Find eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(T)

print("\nEigenvalues:")
for i, val in enumerate(eigenvalues):
    print(f"λ_{i+1} = {val:.4f}")

print("\nEigenvectors (as columns):")
print(eigenvectors)

# Part 4: Visualize eigenvectors and their transformations
def plot_eigenvectors(T, eigenvalues, eigenvectors):
    # Create a figure
    fig, ax = plt.subplots(figsize=(8, 8))
    
    # Set limits and grid
    ax.set_xlim(-3, 3)
    ax.set_ylim(-3, 3)
    ax.grid(True)
    ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    ax.axvline(x=0, color='k', linestyle='-', alpha=0.3)
    
    # Colors for the eigenvectors
    colors = ['blue', 'green']
    
    # Plot eigenvectors and their transformations
    for i in range(len(eigenvalues)):
        # Get the i-th eigenvector
        v = eigenvectors[:, i]
        
        # Scale eigenvector for better visualization
        v_scaled = v * 2
        
        # Plot the original eigenvector
        ax.arrow(0, 0, v_scaled[0], v_scaled[1], head_width=0.1, head_length=0.1, 
                 fc=colors[i], ec=colors[i], label=f"Eigenvector {i+1}")
        
        # Transform the eigenvector
        tv = np.dot(T, v)
        tv_scaled = tv * 2
        
        # Plot the transformed eigenvector
        ax.arrow(0, 0, tv_scaled[0], tv_scaled[1], head_width=0.1, head_length=0.1, 
                 fc='lightcoral', ec='lightcoral', linestyle='--')
        
        # Add text for eigenvalue
        ax.text(v_scaled[0]*1.1, v_scaled[1]*1.1, 
                f"λ_{i+1}={eigenvalues[i]:.2f}", 
                color=colors[i])
    
    # Add a unit circle and its transformation
    theta = np.linspace(0, 2*np.pi, 100)
    circle_x = np.cos(theta)
    circle_y = np.sin(theta)
    
    # Stack coordinates to create points
    circle_points = np.column_stack((circle_x, circle_y))
    
    # Transform the circle
    transformed_circle = np.dot(circle_points, T)
    
    # Plot the unit circle
    ax.plot(circle_x, circle_y, 'k-', alpha=0.3)
    
    # Plot the transformed circle
    ax.plot(transformed_circle[:, 0], transformed_circle[:, 1], 'r-', alpha=0.3)
    
    ax.set_title('Eigenvectors and Their Transformations')
    ax.legend()
    plt.axis('equal')
    plt.tight_layout()
    plt.show()

# Visualize eigenvectors
plot_eigenvectors(T, eigenvalues, eigenvectors)

# Part 5: Visualize the transformation in 3D
def plot_transformation_3d():
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a grid of points
    x = np.linspace(-2, 2, 10)
    y = np.linspace(-2, 2, 10)
    X, Y = np.meshgrid(x, y)
    
    # Flatten the grid points
    points = np.column_stack((X.flatten(), Y.flatten()))
    
    # Transform the points
    transformed_points = np.dot(points, T)
    
    # Reshape the transformed coordinates
    X_transformed = transformed_points[:, 0].reshape(X.shape)
    Y_transformed = transformed_points[:, 1].reshape(Y.shape)
    
    # Create Z coordinates (all zeros for 2D transformation)
    Z = np.zeros(X.shape)
    Z_transformed = np.ones(X.shape)  # Offset for visualization
    
    # Plot the original grid
    ax.plot_surface(X, Y, Z, alpha=0.5, color='blue', label='Original')
    
    # Plot the transformed grid
    ax.plot_surface(X_transformed, Y_transformed, Z_transformed, alpha=0.5, color='red', label='Transformed')
    
    # Add lines connecting original and transformed points
    for i in range(len(points)):
        ax.plot([points[i, 0], transformed_points[i, 0]],
                [points[i, 1], transformed_points[i, 1]],
                [0, 1], 'k-', alpha=0.2)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z (Visualization Only)')
    ax.set_title('3D Visualization of 2D Linear Transformation')
    
    plt.tight_layout()
    plt.show()

# Visualize the transformation in 3D
plot_transformation_3d()

Explanation of the Solution

1. Defining the Linear Transformation

We define a linear transformation $T$ using a $2 \times 2$ matrix:

$$T = \begin{pmatrix} 2 & 1 \ 1 & 3 \end{pmatrix}$$

This matrix represents the transformation, and we can apply it to any vector $\vec{v} \in \mathbb{R}^2$ by matrix multiplication: $T\vec{v}$.

2. Transforming Points

When we apply the transformation to specific points, we get:

  • The point $(1, 0)$ transforms to $(2, 1)$
  • The point $(0, 1)$ transforms to $(1, 3)$
  • The point $(1, 1)$ transforms to $(3, 4)$
  • The point $(-1, 2)$ transforms to $(1, 5)$

This demonstrates how the linear transformation maps points from the input space to the output space.

3. Visualization of the Transformation

Linear Transformation Matrix T:
[[2 1]
 [1 3]]

Original Points:
[1, 0]
[0, 1]
[1, 1]
[-1, 2]

Transformed Points:
[2, 1]
[1, 3]
[3, 4]
[0, 5]

The code visualizes how the unit square (with vertices at $(0,0)$, $(1,0)$, $(1,1)$, and $(0,1)$) is transformed by $T$.

As you can see in the first visualization, the square is stretched, rotated, and skewed by the transformation.

4. Eigenvalues and Eigenvectors

Eigenvalues:
λ_1 = 1.3820
λ_2 = 3.6180

Eigenvectors (as columns):
[[-0.85065081 -0.52573111]
 [ 0.52573111 -0.85065081]]

The eigenvalues and eigenvectors of the transformation matrix tell us about the fundamental behavior of the transformation:

  • The eigenvalues are approximately $\lambda_1 \approx 1.38$ and $\lambda_2 \approx 3.62$
  • The corresponding eigenvectors are shown in the output

An eigenvector is a special vector that, when transformed by $T$, only gets scaled by its corresponding eigenvalue without changing direction.

The second visualization shows these eigenvectors and how they’re affected by the transformation.

5. 3D Visualization

The final visualization provides a 3D perspective of the transformation, showing how a grid of points in the original space maps to points in the transformed space.

Key Insights

  1. Direction of Maximum Stretching:
    The eigenvector corresponding to the larger eigenvalue ($\lambda_2 \approx 3.62$) indicates the direction in which the transformation stretches vectors the most.

  2. Shape Distortion:
    The transformation turns the square into a parallelogram, demonstrating how linear transformations preserve straight lines but can alter angles and distances.

  3. Area Scaling:
    The determinant of the transformation matrix ($\det(T) = 5$) tells us that the transformation scales areas by a factor of 5, which you can observe in the increased size of the transformed square.

  4. Eigenvector Behavior:
    Note how the eigenvectors are only scaled (not rotated) when the transformation is applied.
    This property makes eigenvectors particularly useful in understanding linear transformations.

This example demonstrates the fundamental concepts of linear transformations and how we can visualize and analyze them using Python’s numerical and visualization libraries.

Exploring Nonlinear Dynamics

The Logistic Map Bifurcation Diagram in Python

Let’s dive into an example problem in nonlinear dynamics and solve it using Python in Google Colaboratory.

We’ll consider the logistic map, a classic example of a nonlinear dynamical system that exhibits rich behavior, including chaos, depending on its parameter values.

The logistic map is defined by the recursive equation:

$$
x_{n+1} = r x_n (1 - x_n)
$$

where:

  • $(x_n)$ is the state variable at iteration $(n)$, constrained between $0$ and $1$ (e.g., representing a normalized population),
  • $(r)$ is a growth rate parameter, typically between $0$ and $4$,
  • The nonlinearity arises from the quadratic term $(x_n (1 - x_n))$.

Problem Statement

We’ll simulate the logistic map for different values of $(r)$ and visualize how the system transitions from stable fixed points to periodic orbits and eventually to chaotic behavior. Specifically, we’ll:

  1. Iterate the map for a range of $(r)$ values (e.g., $2.5$ to $4.0$),
  2. Plot the long-term behavior (steady states or oscillations) as a bifurcation diagram,
  3. Use Python to compute and graph the results.

Python Code

Here’s the complete code to simulate and visualize the logistic map’s bifurcation diagram:

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
import numpy as np
import matplotlib.pyplot as plt

# Define the logistic map function
def logistic_map(r, x):
    return r * x * (1 - x)

# Parameters
n_iterations = 1000  # Total iterations
n_transient = 800    # Discard initial transient iterations
n_points = 200       # Number of points to plot per r value
r_values = np.linspace(2.5, 4.0, 1000)  # Range of r values
x0 = 0.1  # Initial condition (must be between 0 and 1)

# Arrays to store results
r_plot = []
x_plot = []

# Iterate over r values
for r in r_values:
    x = x0
    # Run iterations, discard transients
    for i in range(n_iterations):
        x = logistic_map(r, x)
        # Store the last n_points after transient
        if i >= n_transient:
            r_plot.append(r)
            x_plot.append(x)

# Create the bifurcation diagram
plt.figure(figsize=(10, 6))
plt.plot(r_plot, x_plot, ',k', alpha=0.25)  # Plot as small black dots
plt.xlabel('$r$')
plt.ylabel('$x$')
plt.title('Bifurcation Diagram of the Logistic Map')
plt.grid(True)
plt.show()

Code Explanation

  1. Imports:

    • numpy for numerical computations (e.g., creating arrays of $(r)$ values).
    • matplotlib.pyplot for plotting the bifurcation diagram.

  2. Logistic Map Function:

    • logistic_map(r, x) implements the equation $(x_{n+1} = r x_n (1 - x_n))$.

  3. Parameters:

    • n_iterations: Total number of iterations per $(r)$ value.
    • n_transient: Number of initial iterations to discard (to focus on steady-state behavior).
    • n_points: Number of points to plot after transients.
    • r_values: Array of $(r)$ from $2.5$ to $4.0$ with $1000$ steps.
    • x0: Initial condition set to $0.1$ (any value between $0$ and $1$ works).

  4. Simulation:

    • For each $(r)$, start with $(x = x0)$.
    • Iterate the map n_iterations times.
    • After n_transient iterations, store the subsequent $(x)$ values and corresponding $(r)$ values.

  5. Plotting:

    • Use plt.plot with a comma marker (,) to plot points as tiny dots, creating a dense bifurcation diagram.
    • Add labels and a title with TEX-formatted variables (e.g., $r$).

Results and Graph Interpretation

When you run this code in Google Colaboratory, it generates a bifurcation diagram. Here’s what you’ll see:

  • For $(r < 3)$: The system converges to a single fixed point (e.g., at $(r = 2.5)$, $(x)$ stabilizes around $0.6$). This is a stable equilibrium.
  • For $(3 < r < 3.57)$: Period-doubling occurs. At $(r = 3.2)$, you’ll see two values (a period-2 cycle), then four (period-4), and so on.
  • For $(r > 3.57)$: The system enters chaotic regimes, interspersed with periodic windows (e.g., a clear period-3 window around $(r = 3.83)$).
  • At $(r = 4)$: Fully chaotic behavior, with $(x)$ values filling the interval [0, 1] densely.

The graph shows $(x)$ versus $(r)$.
Initially, you see a single line (fixed point), which splits into two, four, etc., before becoming a scattered cloud of points in the chaotic region.
This beautifully illustrates how a simple nonlinear equation can produce complex dynamics.

Conclusion

This example demonstrates key features of nonlinear dynamics: stability, bifurcations, and chaos.

The Python code is efficient and leverages NumPy for fast computation and Matplotlib for clear visualization.

You can tweak parameters like r_values or n_iterations to explore further!

Geodesic Trajectories on a Torus

Computational Analysis and Visualization in Differential Geometry

I’ll provide an example problem in differential geometry, solve it using $Python$, and visualize the results with graphs.

I’ll include both the code and mathematical expressions.

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# Define the torus parametrization
def torus(u, v, R=2, r=1):
    """
    Parametrize a torus with major radius R and minor radius r
    u, v are the parameters in [0, 2π)
    """
    x = (R + r * np.cos(v)) * np.cos(u)
    y = (R + r * np.cos(v)) * np.sin(u)
    z = r * np.sin(v)
    return x, y, z

# Compute the metric coefficients for the torus
def torus_metric(u, v, R=2, r=1):
    """
    Compute the metric coefficients E, F, G for the torus at point (u, v)
    """
    E = (R + r * np.cos(v))**2  # g_11
    F = 0                       # g_12 = g_21
    G = r**2                    # g_22
    return E, F, G

# Compute the Christoffel symbols for the torus
def christoffel_symbols(u, v, R=2, r=1):
    """
    Compute the Christoffel symbols for the torus at point (u, v)
    Returns a 2x2x2 array where Gamma[i,j,k] = Gamma^i_{jk}
    """
    E, F, G = torus_metric(u, v, R, r)
    
    # Derivatives of metric coefficients
    dE_du = 0
    dE_dv = -2 * (R + r * np.cos(v)) * r * np.sin(v)
    dF_du = dF_dv = 0
    dG_du = 0
    dG_dv = 0
    
    # Initialize Christoffel symbols
    Gamma = np.zeros((2, 2, 2))
    
    # g^{11} = 1/E, g^{12} = g^{21} = 0, g^{22} = 1/G
    g_inv_11 = 1/E
    g_inv_22 = 1/G
    
    # Gamma^1_{11}
    Gamma[0, 0, 0] = 0.5 * g_inv_11 * dE_du
    
    # Gamma^1_{12} = Gamma^1_{21}
    Gamma[0, 0, 1] = Gamma[0, 1, 0] = 0.5 * g_inv_11 * dE_dv
    
    # Gamma^1_{22}
    Gamma[0, 1, 1] = -0.5 * g_inv_11 * dG_du
    
    # Gamma^2_{11}
    Gamma[1, 0, 0] = -0.5 * g_inv_22 * dE_dv
    
    # Gamma^2_{12} = Gamma^2_{21}
    Gamma[1, 0, 1] = Gamma[1, 1, 0] = 0.5 * g_inv_22 * dF_du
    
    # Gamma^2_{22}
    Gamma[1, 1, 1] = 0.5 * g_inv_22 * dG_dv
    
    return Gamma

# Geodesic equation as a first-order system
def geodesic_equation(t, y, R=2, r=1):
    """
    The geodesic equation as a first-order system
    y = [u, v, u_dot, v_dot]
    Returns [u_dot, v_dot, u_ddot, v_ddot]
    """
    u, v, u_dot, v_dot = y
    
    # Compute Christoffel symbols
    Gamma = christoffel_symbols(u, v, R, r)
    
    # Geodesic equations (second derivatives)
    u_ddot = -Gamma[0, 0, 0] * u_dot**2 - 2 * Gamma[0, 0, 1] * u_dot * v_dot - Gamma[0, 1, 1] * v_dot**2
    v_ddot = -Gamma[1, 0, 0] * u_dot**2 - 2 * Gamma[1, 0, 1] * u_dot * v_dot - Gamma[1, 1, 1] * v_dot**2
    
    return [u_dot, v_dot, u_ddot, v_ddot]

# Solve the geodesic equation with given initial conditions
def solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R=2, r=1):
    """
    Solve the geodesic equation with initial conditions
    u0, v0: initial position
    u_dot0, v_dot0: initial velocity
    t_span: time span [t_start, t_end]
    """
    y0 = [u0, v0, u_dot0, v_dot0]
    
    # Solve the ODE system
    sol = solve_ivp(
        lambda t, y: geodesic_equation(t, y, R, r),
        t_span,
        y0,
        method='RK45',
        rtol=1e-8,
        atol=1e-8,
        dense_output=True
    )
    
    return sol

# Plot the torus and geodesic
def plot_torus_and_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Plot the torus surface and the geodesic curve
    """
    # Create figure
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    # Create a mesh grid for the torus
    u = np.linspace(0, 2*np.pi, 100)
    v = np.linspace(0, 2*np.pi, 100)
    u_grid, v_grid = np.meshgrid(u, v)
    
    # Compute the torus coordinates
    x, y, z = torus(u_grid, v_grid, R, r)
    
    # Plot the torus surface with translucent appearance
    ax.plot_surface(x, y, z, cmap=cm.viridis, alpha=0.6, edgecolor='none')
    
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    
    # Convert the geodesic curve to Cartesian coordinates
    x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
    
    # Plot the geodesic curve
    ax.plot(x_geo, y_geo, z_geo, 'r-', linewidth=2, label='Geodesic')
    
    # Mark the start and end points
    ax.plot([x_geo[0]], [y_geo[0]], [z_geo[0]], 'go', markersize=8, label='Start')
    ax.plot([x_geo[-1]], [y_geo[-1]], [z_geo[-1]], 'bo', markersize=8, label='End')
    
    # Set labels and title
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Geodesic on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()
    
    return fig, ax

# Analyze the length and curvature of the geodesic
def analyze_geodesic(sol, R=2, r=1, num_points=1000):
    """
    Analyze properties of the geodesic such as length and curvature
    """
    # Extract the solution at evenly spaced points
    t = np.linspace(sol.t[0], sol.t[-1], num_points)
    y_interp = sol.sol(t)
    u_sol, v_sol = y_interp[0], y_interp[1]
    u_dot_sol, v_dot_sol = y_interp[2], y_interp[3]
    
    # Compute the metric coefficients along the geodesic
    E, F, G = [], [], []
    for i in range(len(t)):
        E_i, F_i, G_i = torus_metric(u_sol[i], v_sol[i], R, r)
        E.append(E_i)
        F.append(F_i)
        G.append(G_i)
    
    E = np.array(E)
    F = np.array(F)
    G = np.array(G)
    
    # Compute the speed along the geodesic
    speed = np.sqrt(E * u_dot_sol**2 + 2*F * u_dot_sol * v_dot_sol + G * v_dot_sol**2)
    
    # Compute the length of the geodesic
    dt = t[1] - t[0]
    length = np.sum(speed * dt)
    
    # Plot the speed along the geodesic
    plt.figure(figsize=(10, 6))
    plt.plot(t, speed)
    plt.title('Speed along the Geodesic')
    plt.xlabel('t')
    plt.ylabel('Speed')
    plt.grid(True)
    plt.show()
    
    print(f"Geodesic Length: {length:.6f}")
    
    return length

# Main execution code
if __name__ == "__main__":
    # Parameters for the torus
    R = 2.0  # Major radius
    r = 1.0  # Minor radius
    
    # Initial conditions for the geodesic
    u0 = 0.0
    v0 = 0.0
    u_dot0 = 1.0
    v_dot0 = 0.5
    
    # Time span for integration
    t_span = [0, 10]
    
    # Solve the geodesic equation
    sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
    
    # Plot the results
    plot_torus_and_geodesic(sol, R, r)
    
    # Analyze the geodesic
    geodesic_length = analyze_geodesic(sol, R, r)
    
    # Plot a few different geodesics
    plt.figure(figsize=(12, 10))
    
    initial_conditions = [
        (0.0, 0.0, 1.0, 0.0, 'r', 'Horizontal'),
        (0.0, 0.0, 0.0, 1.0, 'g', 'Vertical'),
        (0.0, 0.0, 1.0, 0.5, 'b', 'Diagonal'),
        (0.0, 0.0, 1.0, 2.0, 'purple', 'Steep Diagonal'),
    ]
    
    ax = plt.axes(projection='3d')
    
    # Plot the torus wireframe
    u = np.linspace(0, 2*np.pi, 30)
    v = np.linspace(0, 2*np.pi, 30)
    u_grid, v_grid = np.meshgrid(u, v)
    x, y, z = torus(u_grid, v_grid, R, r)
    ax.plot_surface(x, y, z, color='lightgray', alpha=0.2)
    
    # Compute and plot different geodesics
    for u0, v0, u_dot0, v_dot0, color, label in initial_conditions:
        sol = solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r)
        
        # Extract the solution
        t = np.linspace(sol.t[0], sol.t[-1], 1000)
        y_interp = sol.sol(t)
        u_sol, v_sol = y_interp[0], y_interp[1]
        
        # Convert to Cartesian coordinates
        x_geo, y_geo, z_geo = torus(u_sol, v_sol, R, r)
        
        # Plot the geodesic
        ax.plot(x_geo, y_geo, z_geo, color=color, linewidth=2, label=label)
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Different Geodesics on a Torus (R={}, r={})'.format(R, r))
    ax.legend()
    
    plt.tight_layout()
    plt.show()

Differential Geometry Example: Geodesics on a Torus

In this example, I’ll solve the geodesic equations on a torus surface and visualize the results.

Geodesics are curves that locally minimize distance between points on a manifold, representing the “straightest possible paths” on the surface.

Mathematical Background

Torus Parameterization

A torus with major radius $R$ and minor radius $r$ can be parameterized as:

$$x(u,v) = (R + r\cos v)\cos u$$
$$y(u,v) = (R + r\cos v)\sin u$$
$$z(u,v) = r\sin v$$

where $u, v \in [0, 2\pi)$.

Metric Tensor

The components of the first fundamental form (metric tensor) are:

$$E = g_{11} = (R + r\cos v)^2$$
$$F = g_{12} = g_{21} = 0$$
$$G = g_{22} = r^2$$

Geodesic Equations

The geodesic equations are given by:

$$\frac{d^2u}{dt^2} + \Gamma^1_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^1_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^1_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

$$\frac{d^2v}{dt^2} + \Gamma^2_{11}\left(\frac{du}{dt}\right)^2 + 2\Gamma^2_{12}\frac{du}{dt}\frac{dv}{dt} + \Gamma^2_{22}\left(\frac{dv}{dt}\right)^2 = 0$$

where $\Gamma^i_{jk}$ are the Christoffel symbols.

Code Explanation

  1. Surface Parameterization: I defined the torus using the standard parameterization with major radius $R$ and minor radius $r$.

  2. Metric Calculation: I computed the components of the metric tensor ($E$, $F$, $G$).

  3. Christoffel Symbols: I calculated the Christoffel symbols from the metric and its derivatives.

  4. Geodesic Equations: I formulated the geodesic equations as a first-order ODE system.

  5. Numerical Integration: I used solve_ivp from $SciPy$ to numerically integrate the geodesic equations. This converts our second-order differential equations into a system of first-order ODEs.

  6. Visualization: I created 3D plots showing the torus surface and the geodesic curves on it.

  7. Analysis: I calculated and plotted the speed along the geodesic and computed its total length.

Key Functions in the Code

  • torus(u, v, R, r): Computes the Cartesian coordinates (x, y, z) for given parameters (u, v)
  • torus_metric(u, v, R, r): Calculates the metric tensor components E, F, G
  • christoffel_symbols(u, v, R, r): Computes the Christoffel symbols for the torus
  • geodesic_equation(t, y, R, r): Implements the geodesic equations as a first-order system
  • solve_geodesic(u0, v0, u_dot0, v_dot0, t_span, R, r): Solves the geodesic initial value problem
  • plot_torus_and_geodesic(sol, R, r, num_points): Visualizes the torus and geodesic
  • analyze_geodesic(sol, R, r, num_points): Computes geodesic properties like length and speed

Results Interpretation

When you run the code, you’ll see several visualizations:

  1. Main Geodesic Visualization:

    Shows a single geodesic curve (in red) on the torus surface, with green and blue markers for the start and end points.

  2. Speed Profile:

    A 2D plot showing how the speed of a particle following the geodesic changes over time.
    For a geodesic with proper parameterization, this speed should be constant, but numerical integration can introduce small variations.

  3. Multiple Geodesics:

    A comparison of different geodesics with various initial velocities:

    • Horizontal geodesic ($u_dot0 = 1, v_dot0 = 0$)
    • Vertical geodesic ($u_dot0 = 0, v_dot0 = 1$)
    • Diagonal geodesic ($u_dot0 = 1, v_dot0 = 0.5$)
    • Steep diagonal geodesic ($u_dot0 = 1, v_dot0 = 2.0$)

Mathematical Insights

  1. Horizontal and Vertical Geodesics: Special cases where the geodesics follow circles on the torus.
    When $u_dot0 = 1$ and $v_dot0 = 0$, the geodesic follows a circle of constant v.
    When $u_dot0 = 0$ and $v_dot0 = 1$, it follows a circle of constant u.

  2. Diagonal Geodesics: More complex paths that wind around the torus, potentially never closing back on themselves (creating dense coverage if extended long enough).

  3. Geodesic Behavior: Notice how geodesics tend to avoid the inner part of the torus and prefer to travel along the outer part where the curvature is less extreme.
    This is a direct consequence of the geodesic equations trying to minimize the path length.

Applications

This type of analysis has applications in:

  • Computer graphics (for path planning on surfaces)
  • General relativity (geodesics represent the paths of free-falling particles)
  • Robotics (finding optimal paths on curved surfaces)
  • Differential geometry (studying the intrinsic properties of surfaces)

You can modify the code to explore different surfaces by changing the parameterization, metric, and Christoffel symbols accordingly.

The Central Limit Theorem

A Visual Exploration with Python

I’ll solve a probability theory problem using $Python$, explain the mathematical concepts, and visualize the results with graphs.

Let’s look at a classic probability problem: the Central Limit Theorem demonstration.

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Set random seed for reproducibility
np.random.seed(42)

# Parameters
n_samples = 10000  # Number of samples to generate
sample_sizes = [1, 2, 5, 10, 30, 50]  # Different sample sizes to demonstrate CLT
distributions = {
    'uniform': stats.uniform(0, 1),
    'exponential': stats.expon(scale=1),
    'poisson': stats.poisson(5)
}

# Function to generate samples and calculate means
def generate_sample_means(distribution, sample_size, n_samples):
    # Generate n_samples sets of sample_size samples and calculate their means
    sample_means = np.array([
        np.mean(distribution.rvs(size=sample_size)) 
        for _ in range(n_samples)
    ])
    return sample_means

# Setup the plots
fig, axes = plt.subplots(len(distributions), len(sample_sizes), figsize=(20, 12))
fig.suptitle('Central Limit Theorem Demonstration', fontsize=16)

# Create plots for each distribution and sample size
for i, (dist_name, distribution) in enumerate(distributions.items()):
    # Calculate the theoretical mean and std of the distribution
    if dist_name == 'uniform':
        true_mean = 0.5
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_mean = 1
        true_std = 1
    else:  # poisson
        true_mean = 5
        true_std = np.sqrt(5)
    
    for j, sample_size in enumerate(sample_sizes):
        # Generate sample means
        means = generate_sample_means(distribution, sample_size, n_samples)
        
        # Calculate expected standard deviation of the sampling distribution
        expected_std = true_std / np.sqrt(sample_size)
        
        # Plot histogram
        ax = axes[i, j]
        ax.hist(means, bins=50, density=True, alpha=0.6, color='skyblue')
        
        # Plot the normal distribution curve for comparison
        x = np.linspace(min(means), max(means), 1000)
        normal_pdf = stats.norm.pdf(x, true_mean, expected_std)
        ax.plot(x, normal_pdf, 'r-', linewidth=2)
        
        # Add theoretical mean line
        ax.axvline(true_mean, color='green', linestyle='--', linewidth=1.5)
        
        # Set title and labels
        if i == 0:
            ax.set_title(f'Sample Size = {sample_size}')
        if j == 0:
            ax.set_ylabel(f'{dist_name.capitalize()}')
        
        # Add some stats
        sample_mean = np.mean(means)
        sample_std = np.std(means)
        ax.text(0.05, 0.95, f'Mean: {sample_mean:.4f}\nSTD: {sample_std:.4f}', 
                transform=ax.transAxes, fontsize=8, 
                verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))

plt.tight_layout()
plt.subplots_adjust(top=0.9)

# Create additional visualization for convergence of variance
plt.figure(figsize=(10, 6))
colors = ['blue', 'red', 'green']
markers = ['o', 's', '^']

for i, (dist_name, distribution) in enumerate(distributions.items()):
    if dist_name == 'uniform':
        true_std = 1/np.sqrt(12)
    elif dist_name == 'exponential':
        true_std = 1
    else:  # poisson
        true_std = np.sqrt(5)
    
    std_ratios = []
    for sample_size in range(1, 51):
        means = generate_sample_means(distribution, sample_size, n_samples)
        sample_std = np.std(means)
        expected_std = true_std / np.sqrt(sample_size)
        std_ratio = sample_std / expected_std
        std_ratios.append(std_ratio)
    
    plt.plot(range(1, 51), std_ratios, color=colors[i], marker=markers[i], 
             linestyle='-', label=f'{dist_name.capitalize()}', markersize=4)

plt.axhline(1, color='black', linestyle='--')
plt.xlabel('Sample Size')
plt.ylabel('Ratio of Sample STD to Expected STD')
plt.title('Convergence of Sampling Distribution Standard Deviation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Let's also visualize the original distributions
plt.figure(figsize=(15, 5))
x_values = {
    'uniform': np.linspace(0, 1, 1000),
    'exponential': np.linspace(0, 5, 1000),
    'poisson': np.arange(0, 20)
}

for i, (dist_name, distribution) in enumerate(distributions.items()):
    plt.subplot(1, 3, i+1)
    
    if dist_name == 'poisson':
        # For discrete distribution
        x = x_values[dist_name]
        plt.bar(x, distribution.pmf(x), alpha=0.6, color=colors[i])
    else:
        # For continuous distributions
        x = x_values[dist_name]
        plt.plot(x, distribution.pdf(x), color=colors[i], linewidth=2)
        plt.fill_between(x, distribution.pdf(x), alpha=0.3, color=colors[i])
    
    plt.title(f'{dist_name.capitalize()} Distribution')
    plt.xlabel('Value')
    plt.ylabel('Density' if dist_name != 'poisson' else 'Probability')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Central Limit Theorem Example

The Central Limit Theorem (CLT) is one of the most important results in probability theory.

It states that when independent random variables are added, their properly normalized sum tends toward a normal distribution even if the original variables are not normally distributed.

Mathematical Formulation

Let $X_1, X_2, \ldots, X_n$ be independent and identically distributed random variables with mean $\mu$ and finite variance $\sigma^2$.

The CLT states that as $n$ approaches infinity, the distribution of:

$$Z_n = \frac{\bar{X}_n - \mu}{\sigma/\sqrt{n}}$$

converges to the standard normal distribution $N(0,1)$,
where is the sample mean.

In other words, the sampling distribution of the mean approaches a normal distribution with mean $\mu$ and standard deviation $\frac{\sigma}{\sqrt{n}}$ as the sample size increases.

Code Explanation

My code demonstrates the Central Limit Theorem by:

  1. Generating samples from three different distributions (uniform, exponential, and Poisson)
  2. Computing the mean of each sample for various sample sizes
  3. Plotting the distribution of these sample means alongside a normal distribution curve
  4. Tracking how quickly the standard deviation of the sample means converges to the theoretical value

Key components:

  • generate_sample_means() creates multiple samples of a specified size and calculates their means
  • The first plot shows histograms of sample means for different distributions and sample sizes
  • The second plot shows how the standard deviation ratio converges to the expected value
  • The third plot visualizes the original distributions to highlight their different shapes

Theoretical Insights

For each distribution, I’ve calculated the true mean and standard deviation:

  1. Uniform distribution on [0,1]:

    • Mean: $\mu = \frac{0+1}{2} = 0.5$
    • Variance: $\sigma^2 = \frac{(1-0)^2}{12} = \frac{1}{12}$

  2. Exponential distribution with scale=1:

    • Mean: $\mu = \lambda = 1$
    • Variance: $\sigma^2 = \lambda^2 = 1$

  3. Poisson distribution with λ=5:

    • Mean: $\mu = \lambda = 5$
    • Variance: $\sigma^2 = \lambda = 5$

The standard deviation of the sampling distribution should be $\frac{\sigma}{\sqrt{n}}$ where n is the sample size.

Results Interpretation

The graphs demonstrate several important aspects of the CLT:



  1. As sample size increases, the distribution of sample means becomes more normal (bell-shaped), regardless of the original distribution shape
  2. The larger the sample size, the smaller the standard deviation of the sampling distribution
  3. The convergence happens relatively quickly - even with sample sizes of $30$, the approximation is quite good
  4. The standard deviation of the sample means scales proportionally to $\frac{1}{\sqrt{n}}$

The bottom plot showing the ratio of observed to theoretical standard deviation confirms that our experimental results match the theory as it approaches $1.0$ for all distributions as sample size increases.

This demonstrates why the CLT is so powerful in statistics - it allows us to make inferences about population parameters regardless of the underlying distribution, as long as we have a sufficiently large sample size.

Solving the Lotka-Volterra Predator-Prey Model

Mathematical Analysis and Python Visualization

I’ll provide an example problem in biomathematics, solve it using $Python$, and include visualizations to illustrate the results.

Let me walk you through a classic example: the Lotka-Volterra predator-prey model.

Lotka-Volterra Predator-Prey Model

The Lotka-Volterra equations describe the dynamics of biological systems in which two species interact, one as a predator and one as prey.

This model is fundamental in mathematical biology.

The basic equations are:

$$\frac{dx}{dt} = \alpha x - \beta xy$$
$$\frac{dy}{dt} = \delta xy - \gamma y$$

Where:

  • $x$ is the prey population
  • $y$ is the predator population
  • $\alpha$ is the prey’s natural growth rate
  • $\beta$ is the predation rate
  • $\delta$ is the predator’s reproduction rate (proportional to food intake)
  • $\gamma$ is the predator’s natural mortality rate

Let’s solve this system of differential equations numerically using $Python$:

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib import cm

# Define the Lotka-Volterra model
def lotka_volterra(state, t, alpha, beta, delta, gamma):
    x, y = state
    dx_dt = alpha * x - beta * x * y
    dy_dt = delta * x * y - gamma * y
    return [dx_dt, dy_dt]

# Set parameter values
alpha = 1.0   # Prey growth rate
beta = 0.1    # Predation rate
delta = 0.02  # Conversion efficiency (predator reproduction)
gamma = 0.4   # Predator mortality rate

# Initial conditions
x0 = 40       # Initial prey population
y0 = 9        # Initial predator population
state0 = [x0, y0]

# Time points to solve for
t = np.linspace(0, 100, 1000)

# Solve the differential equations
solution = odeint(lotka_volterra, state0, t, args=(alpha, beta, delta, gamma))
x_solution = solution[:, 0]  # Prey population over time
y_solution = solution[:, 1]  # Predator population over time

# Create a figure with multiple plots
plt.figure(figsize=(16, 12))

# Plot 1: Population vs. Time
plt.subplot(2, 2, 1)
plt.plot(t, x_solution, 'g-', label='Prey (x)')
plt.plot(t, y_solution, 'r-', label='Predator (y)')
plt.grid(True)
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Population Dynamics Over Time')
plt.legend()

# Plot 2: Phase Plot
plt.subplot(2, 2, 2)
plt.plot(x_solution, y_solution, 'b-')
plt.grid(True)
plt.xlabel('Prey Population (x)')
plt.ylabel('Predator Population (y)')
plt.title('Phase Plot: Predator vs. Prey')

# Mark the initial point
plt.plot(x0, y0, 'ro', label='Initial State')
plt.legend()

# Plot 3: 3D visualization of the dynamics
ax = plt.subplot(2, 2, 3, projection='3d')
ax.plot(x_solution, y_solution, t, color='purple')
ax.set_xlabel('Prey Population (x)')
ax.set_ylabel('Predator Population (y)')
ax.set_zlabel('Time (t)')
ax.set_title('3D Trajectory of the System')

# Plot 4: Vector field showing the dynamics
plt.subplot(2, 2, 4)
x_range = np.linspace(0.5 * min(x_solution), 1.5 * max(x_solution), 20)
y_range = np.linspace(0.5 * min(y_solution), 1.5 * max(y_solution), 20)
X, Y = np.meshgrid(x_range, y_range)

# Calculate derivatives
DX = alpha * X - beta * X * Y
DY = delta * X * Y - gamma * Y

# Normalize the vectors for better visualization
norm = np.sqrt(DX**2 + DY**2)
# Avoid division by zero
norm[norm == 0] = 1
DX = DX / norm
DY = DY / norm

plt.quiver(X, Y, DX, DY, norm, cmap=cm.viridis)
plt.colorbar(label='Vector Magnitude')
plt.xlabel('Prey Population (x)')
plt.ylabel('Predator Population (y)')
plt.title('Vector Field of the System')

# Overlay the solution trajectory
plt.plot(x_solution, y_solution, 'r-', alpha=0.3, label='Trajectory')
plt.legend()

plt.tight_layout()
plt.savefig('lotka_volterra_analysis.png', dpi=300)
plt.show()

# Calculate statistics about the cycles
from scipy.signal import find_peaks

# Find peaks in prey population to analyze cycles
peaks, _ = find_peaks(x_solution, height=np.mean(x_solution))
peak_times = t[peaks]
cycle_periods = np.diff(peak_times)

print(f"Average cycle period: {np.mean(cycle_periods):.2f} time units")
print(f"Maximum prey population: {np.max(x_solution):.2f}")
print(f"Maximum predator population: {np.max(y_solution):.2f}")
print(f"Minimum prey population: {np.min(x_solution):.2f}")
print(f"Minimum predator population: {np.min(y_solution):.2f}")

Code Explanation

  1. Model Definition: I’ve defined the Lotka-Volterra differential equations in a function that returns the rate of change for both populations.

  2. Parameters:

    • $α (alpha) = 1.0$: Growth rate of prey in absence of predators
    • $β (beta) = 0.1$: Rate at which predators kill prey
    • $δ (delta) = 0.02$: Rate at which predators increase by consuming prey
    • $γ (gamma) = 0.4$: Natural death rate of predators

  3. Numerical Solution: Using $SciPy$’s odeint to solve the differential equations over time.

  4. Visualization: The code creates four different plots:

    • Time series of both populations
    • Phase plot showing the predator-prey cycle
    • 3D trajectory showing how the system evolves
    • Vector field that illustrates the dynamics of the system

  5. Analysis: The code calculates statistics about the cycles, such as average period and population extremes.

Results Interpretation

When you run this code, you’ll observe the classic predator-prey relationship:

Average cycle period: 10.30 time units
Maximum prey population: 40.52
Maximum predator population: 15.95
Minimum prey population: 7.95
Minimum predator population: 5.75

  1. Cyclical Behavior: Both populations oscillate, but out of phase with each other.
    When prey increase, predators follow with a delay. More predators reduce the prey population, which then causes predator decline.

  2. Phase Plot: The closed loops in the phase plot indicate periodic behavior - the system returns to similar states over time rather than reaching a stable equilibrium.

  3. Vector Field: The arrows show the direction of change at each point in the state space.
    This helps understand how the system would evolve from any initial condition.

  4. Ecological Insights:

    • The system exhibits natural oscillations without external forcing
    • The predator population peaks after the prey population
    • Neither species goes extinct under these parameters
    • The system is sensitive to initial conditions and parameter values

This model, despite its simplicity, captures fundamental dynamics of predator-prey systems observed in nature, such as the famous case of lynx and snowshoe hare populations in Canada, which show similar oscillatory patterns.

You can modify the parameters to explore different ecological scenarios, such as:

  • Higher predation rates (increase $β$)
  • Lower predator mortality (decrease $γ$)
  • Different initial population ratios

These changes would result in different dynamics, potentially including extinction of one or both species under extreme parameter values.

Integer Programming：Solving the Knapsack Problem in Python with PuLP

I’ll solve an integer programming example using $Python$, showing the implementation, solution, and visualization.

Let’s solve a classic integer programming problem - the Knapsack Problem.

Integer Programming Example: Knapsack Problem

The knapsack problem is a classic optimization problem: Given a set of items, each with a weight and a value, determine which items to include in a collection so that the total weight is less than or equal to a given limit and the total value is maximized.

Mathematically, we can formulate this as:

$$\begin{align}
\text{Maximize} & \sum_{i=1}^{n} v_i x_i \
\text{Subject to} & \sum_{i=1}^{n} w_i x_i \leq W \
& x_i \in {0,1} \quad \forall i \in {1,2,\ldots,n}
\end{align}$$

Where:

  • $v_i$ is the value of item $i$
  • $w_i$ is the weight of item $i$
  • $x_i$ is a binary variable indicating whether item $i$ is selected
  • $W$ is the maximum weight capacity

Let’s implement this using $PuLP$, a popular linear programming library in $Python$:

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pulp import *

# Problem data
items = ["Item_" + str(i) for i in range(1, 11)]
values = [10, 13, 18, 31, 7, 15, 28, 22, 11, 17]
weights = [5, 8, 12, 16, 3, 9, 14, 11, 5, 8]
capacity = 40

# Create a dataframe for better visualization
data = pd.DataFrame({
    'Item': items,
    'Value': values,
    'Weight': weights,
    'Value/Weight': [v/w for v, w in zip(values, weights)]
})

# Create the LP problem
prob = LpProblem("Knapsack_Problem", LpMaximize)

# Create binary variables for each item
x = LpVariable.dicts("Include", items, lowBound=0, upBound=1, cat='Integer')

# Objective function: maximize total value
prob += lpSum([values[i] * x[items[i]] for i in range(len(items))]), "Total Value"

# Constraint: total weight must not exceed capacity
prob += lpSum([weights[i] * x[items[i]] for i in range(len(items))]) <= capacity, "Weight_Constraint"

# Solve the problem
prob.solve()

# Extract results
selected_items = []
selected_values = []
selected_weights = []
total_value = 0
total_weight = 0

for i in range(len(items)):
    if value(x[items[i]]) == 1:
        selected_items.append(items[i])
        selected_values.append(values[i])
        selected_weights.append(weights[i])
        total_value += values[i]
        total_weight += weights[i]

# Create results dataframe
results = pd.DataFrame({
    'Item': selected_items,
    'Value': selected_values,
    'Weight': selected_weights
})

# Print solution information
print("Status:", LpStatus[prob.status])
print("\nSelected Items:")
print(results)
print("\nTotal Value:", total_value)
print("Total Weight:", total_weight)
print("Capacity:", capacity)

# Visualization of the solution
plt.figure(figsize=(12, 10))

# Plot 1: All items with Value vs Weight
plt.subplot(2, 2, 1)
sns.scatterplot(x='Weight', y='Value', data=data, s=100)
for i, row in data.iterrows():
    plt.annotate(row['Item'], (row['Weight'], row['Value']), fontsize=9)
plt.title('All Items: Value vs Weight')
plt.grid(True, alpha=0.3)

# Plot 2: Selected vs Not Selected
plt.subplot(2, 2, 2)
selected_status = ['Selected' if item in selected_items else 'Not Selected' for item in items]
plot_data = data.copy()
plot_data['Status'] = selected_status
sns.scatterplot(x='Weight', y='Value', hue='Status', data=plot_data, s=100)
plt.title('Selected vs Not Selected Items')
plt.grid(True, alpha=0.3)

# Plot 3: Value/Weight Efficiency
plt.subplot(2, 2, 3)
sns.barplot(x='Item', y='Value/Weight', data=data)
plt.xticks(rotation=45)
plt.title('Value/Weight Ratio of Items')
plt.grid(True, alpha=0.3)

# Plot 4: Contribution to Solution
plt.subplot(2, 2, 4)
if not results.empty:
    results['Percentage'] = results['Value'] / total_value * 100
    sns.barplot(x='Item', y='Percentage', data=results)
    plt.xticks(rotation=45)
    plt.title('Contribution to Total Value (%)')
    plt.ylabel('Percentage of Total Value')
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Code Explanation

  1. Setting up the problem:

    • We defined $10$ items with corresponding values and weights
    • The knapsack capacity is set to $40$ units

  2. Formulating the integer program with PuLP:

    • Created binary decision variables for each item (x)
    • Set the objective function to maximize total value
    • Added a constraint to ensure the total weight doesn’t exceed capacity

  3. Solving the problem:

    • $PuLP$ uses branch and bound algorithms to solve the integer program
    • Results are extracted and organized into a pandas DataFrame

  4. Visualization:

    • Four plots are created to analyze different aspects of the solution:
      • All items plotted by weight vs. value
      • Selected vs. not selected items
      • Value/weight efficiency ratio for each item
      • Each selected item’s contribution to the total value

Solution Analysis

When executed, this code will select the optimal combination of items that maximizes value while staying within the weight limit.

The optimization model makes trade-offs based on both the absolute value of items and their weight efficiency (value per unit weight).

Items with high values and relatively low weights are preferred.

The visualizations help us understand:

  1. The distribution of all items in the value-weight space
  2. Which items were selected by the optimizer
  3. The efficiency (value/weight) of each item
  4. How much each selected item contributes to the total solution value

This example demonstrates how integer programming can be applied to combinatorial optimization problems, finding the optimal solution among a finite but potentially very large set of possibilities.

Result

Status: Optimal

Selected Items:
     Item  Value  Weight
0  Item_1     10       5
1  Item_4     31      16
2  Item_5      7       3
3  Item_8     22      11
4  Item_9     11       5

Total Value: 81
Total Weight: 40
Capacity: 40

Calculating the Expected Present Value of a Life Insurance Policy Using Python

Problem: Expected Present Value of an Insurance Policy

A whole life insurance policy pays a benefit of $100,000 upon the policyholder’s death.

The probability of death at age $( x )$ is given by the mortality table.

The force of interest is $( \delta = 0.05 )$.

Compute the expected present value (EPV) of the policy for a policyholder aged $30$ using a simplified mortality model.

Mathematical Formulation

The expected present value (EPV) is given by:

$$
EPV = \sum_{t=1}^{\infty} P(T_x = t) \cdot 100,000 \cdot e^{-\delta t}
$$

where:

  • $( T_x )$ is the future lifetime of the insured aged $( x )$.
  • $ P(T_x = t) $ represents the probability that the insured dies exactly at time $( t )$.
  • $ \delta $ is the force of interest.

For this problem, we assume a simplified mortality model where the probability of surviving from age $( x )$ to age $( x+t )$ follows an exponential distribution:

$$
P(T_x > t) = e^{-\mu t}
$$

where $ \mu = 0.02 $ is the mortality rate.

The probability of dying at time $( t )$ is:

$$
P(T_x = t) = \mu e^{-\mu t}
$$

Now, we implement this in $Python$ and visualize the expected present value.

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
import numpy as np
import matplotlib.pyplot as plt

# Parameters
mu = 0.02  # Mortality rate
delta = 0.05  # Force of interest
benefit = 100000  # Death benefit

# Time range for computation
t = np.arange(1, 100, 1)  # Assume a maximum lifetime of 100 years

# Probability of dying at time t
prob_death = mu * np.exp(-mu * t)

# Discount factor
discount_factor = np.exp(-delta * t)

# Expected Present Value Calculation
EPV = np.sum(prob_death * benefit * discount_factor)

print(f"Expected Present Value (EPV) of the policy: ${EPV:,.2f}")

# Visualization
plt.figure(figsize=(8, 5))
plt.plot(t, prob_death * benefit * discount_factor, label='Present Value of Payout at t')
plt.xlabel('Years')
plt.ylabel('Present Value of Expected Payout')
plt.title('Expected Present Value of Insurance Payout Over Time')
plt.legend()
plt.grid()
plt.show()

Explanation of the Code

  1. We define the mortality rate $( \mu )$, force of interest $( \delta )$, and death benefit.
  2. We create an array of years from $1$ to $100$, assuming the policyholder can live up to $130$ years.
  3. We compute the probability of death at each year using $( P(T_x = t) = \mu e^{-\mu t} )$.
  4. We compute the discount factor $( e^{-\delta t} )$.
  5. We calculate the EPV as the sum of discounted expected payouts.
  6. Finally, we plot the present value of payouts over time to visualize their impact.

Interpretation of the Result

[Result]

Expected Present Value (EPV) of the policy: $27,556.12

  • The EPV represents the fair price an insurer should charge for this policy (ignoring expenses and profit).
  • The graph shows how the present value of potential payouts decreases over time due to discounting.

Differential Geometry

Example Problem in Differential Geometry

Let’s consider a classic problem in differential geometry: finding the geodesic (shortest path) on a unit sphere.

A geodesic on a sphere is a great circle, and we’ll derive its equation using the Euler-Lagrange equations from the calculus of variations.

Then, we’ll solve it numerically in $Python$ and visualize the result.

Problem Statement

On a unit sphere parameterized by spherical coordinates $(\theta, \phi)$, where $\theta$ is the polar angle and $\phi$ is the azimuthal angle, the metric is given by:
$$
ds^2 = d\theta^2 + \sin^2\theta , d\phi^2
$$
We want to find the geodesic between two points, say $(\theta_1, \phi_1) = (0.5, 0)$ and $(\theta_2, \phi_2) = (1.0, 1.0)$, by minimizing the arc length functional:
$$
S = \int_{t_1}^{t_2} \sqrt{\dot{\theta}^2 + \sin^2\theta , \dot{\phi}^2} , dt
$$
where $\dot{\theta} = \frac{d\theta}{dt}$ and $\dot{\phi} = \frac{d\phi}{dt}$.

Step 1: Euler-Lagrange Equations

The Lagrangian is:
$$
L = \sqrt{\dot{\theta}^2 + \sin^2\theta , \dot{\phi}^2}
$$
The Euler-Lagrange equations for $\theta$ and $\phi$ are:
$$
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}} \right) = \frac{\partial L}{\partial \theta}, \quad \frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\phi}} \right) = \frac{\partial L}{\partial \phi}
$$
Computing the partial derivatives:

  • $\frac{\partial L}{\partial \dot{\theta}} = \frac{\dot{\theta}}{\sqrt{\dot{\theta}^2 + \sin^2\theta , \dot{\phi}^2}}$
  • $\frac{\partial L}{\partial \dot{\phi}} = \frac{\sin^2\theta , \dot{\phi}}{\sqrt{\dot{\theta}^2 + \sin^2\theta , \dot{\phi}^2}}$
  • $\frac{\partial L}{\partial \theta} = \frac{\sin\theta \cos\theta , \dot{\phi}^2}{\sqrt{\dot{\theta}^2 + \sin^2\theta , \dot{\phi}^2}}$
  • $\frac{\partial L}{\partial \phi} = 0$ (since $L$ does not depend explicitly on $\phi$)

Since $\frac{\partial L}{\partial \phi} = 0$, the second equation implies that $\frac{\partial L}{\partial \dot{\phi}}$ is a constant, which is a conserved quantity (related to angular momentum).

This simplifies to Clairaut’s relation for the sphere:
$$
\sin^2\theta , \dot{\phi} = \text{constant}
$$
However, solving these analytically is complex, so we’ll use numerical integration in $Python$.

Step 2: Python Code

We’ll use scipy.integrate.odeint to solve the differential equations numerically, then plot the geodesic on the sphere.

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
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Define the system of ODEs for the geodesic
def geodesic(y, t):
    theta, phi, dtheta_dt, dphi_dt = y
    # Lagrangian derivatives lead to these second-order ODEs
    d2theta_dt2 = np.sin(theta) * np.cos(theta) * dphi_dt**2
    d2phi_dt2 = -2 * dtheta_dt * dphi_dt * np.cos(theta) / np.sin(theta)
    return [dtheta_dt, dphi_dt, d2theta_dt2, d2phi_dt2]

# Initial conditions
theta1, phi1 = 0.5, 0.0  # Starting point
theta2, phi2 = 1.0, 1.0  # Ending point
t = np.linspace(0, 1, 100)  # Parameter t from 0 to 1

# Guess initial velocities (tuned to satisfy boundary conditions)
dtheta_dt0 = 0.5
dphi_dt0 = 1.0
y0 = [theta1, phi1, dtheta_dt0, dphi_dt0]

# Solve the ODE
sol = odeint(geodesic, y0, t)

# Extract solutions
theta = sol[:, 0]
phi = sol[:, 1]

# Convert to Cartesian coordinates for plotting
x = np.sin(theta) * np.cos(phi)
y = np.sin(theta) * np.sin(phi)
z = np.cos(theta)

# Plot the sphere and geodesic
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')

# Plot the unit sphere
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x_sphere = np.outer(np.cos(u), np.sin(v))
y_sphere = np.outer(np.sin(u), np.sin(v))
z_sphere = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x_sphere, y_sphere, z_sphere, color='b', alpha=0.1)

# Plot the geodesic
ax.plot(x, y, z, color='r', linewidth=2, label='Geodesic')
ax.scatter([np.sin(theta1) * np.cos(phi1)], [np.sin(theta1) * np.sin(phi1)], [np.cos(theta1)], color='g', s=100, label='Start')
ax.scatter([np.sin(theta2) * np.cos(phi2)], [np.sin(theta2) * np.sin(phi2)], [np.cos(theta2)], color='y', s=100, label='End')

# Labels and legend
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.legend()
plt.title('Geodesic on a Unit Sphere')
plt.show()

Code Explanation

  1. Imports: We use numpy for numerical operations, scipy.integrate.odeint for solving ODEs, and matplotlib for 3D plotting.
  2. Geodesic Function: Defines the system of first-order ODEs derived from the Euler-Lagrange equations.
    We convert the second-order equations into four first-order equations: $\theta$, $\phi$, $\dot{\theta}$, and $\dot{\phi}$.
  3. Initial Conditions: We set the starting point $(\theta_1, \phi_1) = (0.5, 0)$ and guess initial velocities.
    In practice, these velocities would be adjusted to hit the target point $(\theta_2, \phi_2) = (1.0, 1.0)$, but for simplicity, we demonstrate the method.
  4. Numerical Solution: odeint solves the ODE system over the parameter $t$.
  5. Visualization: We convert spherical coordinates to Cartesian coordinates $(x = \sin\theta \cos\phi$, $y = \sin\theta \sin\phi$, $z = \cos\theta$) and plot the geodesic on a semi-transparent unit sphere.

Result and Visualization Explanation

  • The plot shows a unit sphere (blue, semi-transparent) with the geodesic (red curve) connecting the start point (green dot) and an approximate end point (yellow dot).
  • The geodesic is a segment of a great circle, the shortest path on the sphere.
  • Note: The initial velocities here are illustrative.
    To precisely hit $(\theta_2, \phi_2)$, you’d need a boundary value problem solver (e.g., scipy.integrate.solve_bvp), but this demonstrates the concept.

This visualization makes it clear how geodesics behave on curved surfaces, a fundamental concept in differential geometry!