Ordinary Differential Equation (ODE)

January 18, 2025

I’ll help you solve an Ordinary Differential Equation (ODE) using $Python$.

Let’s solve a spring-mass-damper system, which is a classic second-order ODE problem.

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

def spring_mass_damper(state, t, m, k, c):
    """
    Spring-mass-damper system
    Second order ODE: m*d²x/dt² + c*dx/dt + k*x = 0
    Converted to system of first order ODEs:
    dx/dt = v
    dv/dt = -(k/m)x - (c/m)v
    """
    x, v = state
    dxdt = v
    dvdt = -(k/m) * x - (c/m) * v
    return [dxdt, dvdt]

# System parameters
m = 1.0    # mass (kg)
k = 10.0   # spring constant (N/m)
c = 0.5    # damping coefficient (Ns/m)

# Initial conditions
x0 = 1.0   # initial displacement (m)
v0 = 0.0   # initial velocity (m/s)
state0 = [x0, v0]

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

# Solve ODE
solution = odeint(spring_mass_damper, state0, t, args=(m, k, c))
x = solution[:, 0]  # displacement
v = solution[:, 1]  # velocity

# Calculate energy
potential_energy = 0.5 * k * x**2
kinetic_energy = 0.5 * m * v**2
total_energy = potential_energy + kinetic_energy

# Create plots
plt.figure(figsize=(15, 10))

# Position vs Time
plt.subplot(2, 2, 1)
plt.plot(t, x, 'b-', label='Position')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position vs Time')
plt.legend()

# Velocity vs Time
plt.subplot(2, 2, 2)
plt.plot(t, v, 'r-', label='Velocity')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Velocity (m/s)')
plt.title('Velocity vs Time')
plt.legend()

# Phase Portrait
plt.subplot(2, 2, 3)
plt.plot(x, v, 'g-', label='Phase Portrait')
plt.grid(True)
plt.xlabel('Position (m)')
plt.ylabel('Velocity (m/s)')
plt.title('Phase Portrait')
plt.legend()

# Energy vs Time
plt.subplot(2, 2, 4)
plt.plot(t, potential_energy, 'b-', label='Potential Energy')
plt.plot(t, kinetic_energy, 'r-', label='Kinetic Energy')
plt.plot(t, total_energy, 'g-', label='Total Energy')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Energy (J)')
plt.title('Energy vs Time')
plt.legend()

plt.tight_layout()
plt.show()

Let me explain the problem and solution in detail:

The Problem: Spring-Mass-Damper System
- Mathematical form: $m(d²x/dt²) + c(dx/dt) + kx = 0$
- Where:
  - $x$ is displacement
  - $m$ is mass
  - $k$ is spring constant
  - $c$ is damping coefficient
Solution Method:
- Convert 2nd order ODE to system of two 1st order ODEs
- Use scipy’s odeint for numerical integration
- System parameters:
  - Mass = $1.0 kg$
  - Spring constant = $10.0 N/m$
  - Damping coefficient = $0.5 Ns/m$
Visualization (four plots):

Position vs Time:
- Shows oscillatory motion with damping
- Amplitude decreases over time due to damping
Velocity vs Time:
- Shows velocity oscillations
- Amplitude decreases faster than position
Phase Portrait:
- Shows system trajectory in position-velocity space
- Spiral pattern indicates damped oscillation
Energy vs Time:
- Shows potential and kinetic energy exchange
- Total energy decreases due to damping

Physical Interpretation:
- System starts at $x = 1m$ with zero velocity
- Oscillates back and forth while losing energy
- Motion gradually decreases due to damping
- System eventually comes to rest at equilibrium ($x = 0$)
Energy Analysis:
- Energy constantly transfers between potential and kinetic
- Total energy decreases exponentially due to damping
- Potential energy is maximum at maximum displacement
- Kinetic energy is maximum at zero displacement

Partial Differential Equation (PDE)

January 17, 2025

I’ll help you solve a Partial Differential Equation (PDE) example using $Python$.

Let’s solve a 2D heat equation, which is a classic PDE problem.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def solve_2d_heat_equation(dx=0.1, dy=0.1, dt=0.0001, t_final=0.01):
    """
    Solve 2D heat equation using finite difference method
    ∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²)
    with initial conditions: u(x,y,0) = sin(πx)sin(πy)
    and boundary conditions: u = 0 at boundaries
    """
    # Thermal diffusivity
    alpha = 1.0
    
    # Create grid
    x = np.arange(0, 1+dx, dx)
    y = np.arange(0, 1+dy, dy)
    t = np.arange(0, t_final+dt, dt)
    
    # Initialize temperature grid
    nx, ny = len(x), len(y)
    u = np.zeros((nx, ny))
    
    # Set initial conditions
    X, Y = np.meshgrid(x, y)
    u = np.sin(np.pi*X) * np.sin(np.pi*Y)
    
    # Create arrays for new and old solutions
    u_new = np.zeros((nx, ny))
    
    # Stability condition check
    stability = alpha * dt * (1/(dx**2) + 1/(dy**2))
    print(f"Stability parameter: {stability:.3f} (should be < 0.5)")
    
    # Time stepping
    for n in range(len(t)-1):
        # Interior points
        for i in range(1, nx-1):
            for j in range(1, ny-1):
                u_new[i,j] = u[i,j] + alpha*dt*(
                    (u[i+1,j] - 2*u[i,j] + u[i-1,j])/dx**2 +
                    (u[i,j+1] - 2*u[i,j] + u[i,j-1])/dy**2
                )
        
        # Update solution
        u = u_new.copy()
        
        # Enforce boundary conditions
        u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 0
    
    return X, Y, u

# Solve the equation
X, Y, u_final = solve_2d_heat_equation()

# Create visualization
fig = plt.figure(figsize=(12, 5))

# Surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(X, Y, u_final, cmap=cm.viridis)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('Temperature')
ax1.set_title('2D Heat Distribution (3D View)')
fig.colorbar(surf)

# Contour plot
ax2 = fig.add_subplot(122)
contour = ax2.contourf(X, Y, u_final, levels=20, cmap=cm.viridis)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('2D Heat Distribution (Top View)')
fig.colorbar(contour)

plt.tight_layout()
plt.show()

Let me explain this code and the problem we’re solving:

The Problem: We’re solving the 2D heat equation:
- Equation: $∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²)$
- Initial condition: $u(x,y,0) = sin(πx)sin(πy)$
- Boundary conditions: $u = 0$ at all boundaries
- $α$ is the thermal diffusivity (set to $1.0$ in our case)
Solution Method:
- We use the Finite Difference Method (FDM)
- Space is discretized into a grid with steps $dx$ and $dy$
- Time is discretized with step $dt$
- The solution uses an explicit scheme for time-stepping
Code Structure:
- Creates a mesh grid for $x$ and $y$ coordinates
- Sets initial conditions using sine functions
- Implements time-stepping loop with boundary conditions
- Uses $numpy$ for efficient array operations
Visualization:
- Left plot: 3D surface plot showing temperature distribution
- Right plot: 2D contour plot showing temperature from top view
- Colors indicate temperature (darker blue = cooler, yellow = warmer)

Stability parameter: 0.020 (should be < 0.5)

Results Interpretation:
- The initial sine wave pattern gradually diffuses
- Heat spreads from warmer regions to cooler regions
- The solution maintains zero temperature at boundaries
- The pattern becomes smoother over time due to diffusion

The visualization shows both a 3D surface plot and a 2D contour plot to help understand how heat distributes across the 2D surface.

The color gradient represents temperature variations, with warmer colors indicating higher temperatures.

Advanced Fluid Dynamics Simulation

January 16, 2025

A Python-Based Study of Pipe Flow and Pressure Drop Analysis

I’ll create a example of fluid dynamics, this time focusing on pipe flow analysis using the $Hagen$-$Poiseuille$ equation and pressure drop calculations.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

class PipeFlowAnalyzer:
    def __init__(self):
        # Fluid properties (water at 20°C)
        self.density = 998.2      # kg/m³
        self.viscosity = 0.001    # Pa·s
        
        # Pipe properties
        self.length = 10.0        # m
        self.diameter = 0.05      # m
        self.roughness = 0.0015   # mm (typical for commercial steel)
        
    def reynolds_number(self, velocity):
        """Calculate Reynolds number"""
        return self.density * velocity * self.diameter / self.viscosity
    
    def friction_factor(self, reynolds):
        """
        Calculate Darcy friction factor using Colebrook-White equation
        """
        def colebrook(f):
            return (1/np.sqrt(f) + 2.0*np.log10(self.roughness/(3.7*self.diameter*1000) + 
                                               2.51/(reynolds*np.sqrt(f))))
        
        if reynolds < 2300:
            return 64/reynolds
        else:
            # Initial guess using Swamee-Jain equation
            f_0 = 0.25/(np.log10(self.roughness/(3.7*self.diameter*1000) + 5.74/reynolds**0.9)**2)
            return fsolve(colebrook, f_0)[0]
    
    def pressure_drop(self, velocity):
        """Calculate pressure drop using Darcy-Weisbach equation"""
        re = self.reynolds_number(velocity)
        f = self.friction_factor(re)
        return (f * self.length * self.density * velocity**2) / (2 * self.diameter)
    
    def velocity_profile(self, flow_rate, r_points=50):
        """Calculate velocity profile for laminar flow"""
        radius = self.diameter/2
        r = np.linspace(0, radius, r_points)
        
        # Maximum velocity at center
        v_max = 2 * flow_rate / (np.pi * radius**2)
        
        # Parabolic profile
        v = v_max * (1 - (r/radius)**2)
        return r, v
    
    def analyze_flow_range(self, velocities):
        """Analyze flow characteristics over a range of velocities"""
        reynolds = [self.reynolds_number(v) for v in velocities]
        pressure_drops = [self.pressure_drop(v) for v in velocities]
        friction_factors = [self.friction_factor(re) for re in reynolds]
        
        return reynolds, pressure_drops, friction_factors
    
    def plot_analysis(self, velocities):
        """Create visualizations of flow analysis"""
        reynolds, pressure_drops, friction_factors = self.analyze_flow_range(velocities)
        
        fig = plt.figure(figsize=(15, 10))
        
        # Pressure drop vs Velocity
        ax1 = fig.add_subplot(221)
        ax1.plot(velocities, pressure_drops, 'b-', linewidth=2)
        ax1.set_xlabel('Flow Velocity (m/s)')
        ax1.set_ylabel('Pressure Drop (Pa)')
        ax1.set_title('Pressure Drop vs Flow Velocity')
        ax1.grid(True)
        
        # Friction factor vs Reynolds number
        ax2 = fig.add_subplot(222)
        ax2.semilogx(reynolds, friction_factors, 'r-', linewidth=2)
        ax2.set_xlabel('Reynolds Number')
        ax2.set_ylabel('Friction Factor')
        ax2.set_title('Moody Diagram')
        ax2.grid(True)
        
        # Velocity profile for laminar flow
        ax3 = fig.add_subplot(223)
        flow_rate = 0.001  # m³/s
        r, v = self.velocity_profile(flow_rate)
        ax3.plot(r*1000, v, 'g-', linewidth=2)  # Convert radius to mm
        ax3.set_xlabel('Radial Position (mm)')
        ax3.set_ylabel('Velocity (m/s)')
        ax3.set_title('Laminar Flow Velocity Profile')
        ax3.grid(True)
        
        # Head loss vs Flow rate
        ax4 = fig.add_subplot(224)
        head_losses = [dp/(self.density * 9.81) for dp in pressure_drops]
        flow_rates = [v * np.pi * (self.diameter/2)**2 for v in velocities]
        ax4.plot(flow_rates, head_losses, 'm-', linewidth=2)
        ax4.set_xlabel('Flow Rate (m³/s)')
        ax4.set_ylabel('Head Loss (m)')
        ax4.set_title('Head Loss vs Flow Rate')
        ax4.grid(True)
        
        plt.tight_layout()
        return fig

# Create analyzer and run analysis
analyzer = PipeFlowAnalyzer()

# Define velocity range for analysis
velocities = np.linspace(0.1, 5, 100)

# Plot results
analyzer.plot_analysis(velocities)

# Print analysis for specific velocity
test_velocity = 1.0  # m/s
re = analyzer.reynolds_number(test_velocity)
dp = analyzer.pressure_drop(test_velocity)
f = analyzer.friction_factor(re)

print("Pipe Flow Analysis Results:")
print(f"\nTest conditions at velocity = {test_velocity} m/s:")
print(f"Reynolds Number: {re:.0f}")
print(f"Flow Regime: {'Laminar' if re < 2300 else 'Turbulent'}")
print(f"Friction Factor: {f:.4f}")
print(f"Pressure Drop: {dp/1000:.2f} kPa")
print(f"Head Loss: {dp/(analyzer.density * 9.81):.2f} m")

Let me explain this pipe flow analysis example:

Core Fluid Dynamics Concepts:
- $Hagen$-$Poiseuille$ flow
- Reynolds number
- Friction factor
- Pressure drop
- Head loss
Mathematical Components:
- $Darcy$-$Weisbach$ equation
- $Colebrook$-$White$ equation
- Velocity profile calculations
- Head loss computations
Visualization Elements:
- Top left: Pressure drop vs. flow velocity
- Top right: Moody diagram (friction factor vs. Reynolds number)
- Bottom left: Velocity profile in pipe
- Bottom right: Head loss vs. flow rate
Key Analysis Features:
- Flow regime determination
- Friction factor calculation
- Pressure loss prediction
- Velocity distribution
Engineering Applications:
- Pipe system design
- Pump selection
- Energy loss calculations
- Flow rate optimization

Output

Pipe Flow Analysis Results:

Test conditions at velocity = 1.0 m/s:
Reynolds Number: 49910
Flow Regime: Turbulent
Friction Factor: 0.0210
Pressure Drop: 2.10 kPa
Head Loss: 0.21 m

The visualizations help understand:

How pressure drop varies with flow velocity
The relationship between friction factor and Reynolds number
The shape of the velocity profile in the pipe
Energy losses in pipe systems

Mathematical Sociology

January 15, 2025

Example: Modeling Social Influence in a Network

We will simulate how opinions spread in a social network.

The idea is based on the DeGroot model of consensus, where each node in the network updates its opinion based on a weighted average of its neighbors’ opinions.

Problem Setup

A network of individuals is represented as a graph.
Each node has an initial opinion (randomized).
Opinions are updated iteratively using the weighted average of neighboring opinions.
We will observe how opinions converge over time.

Python Code Implementation

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

# Step 1: Create a social network graph
G = nx.erdos_renyi_graph(n=20, p=0.2, seed=42)  # 20 nodes, 20% connection probability

# Step 2: Initialize random opinions
opinions = np.random.rand(len(G.nodes))
initial_opinions = opinions.copy()

# Step 3: Define the adjacency matrix and normalize weights
adj_matrix = nx.to_numpy_array(G)
weights = adj_matrix / adj_matrix.sum(axis=1, keepdims=True)

# Step 4: Simulate opinion updates
iterations = 20
opinion_history = [opinions.copy()]

for _ in range(iterations):
    opinions = weights @ opinions  # Weighted average of neighbors
    opinion_history.append(opinions.copy())

# Step 5: Visualization of the results
plt.figure(figsize=(12, 6))

# Plot the convergence of opinions over time
for i, opinion_track in enumerate(np.array(opinion_history).T):
    plt.plot(opinion_track, label=f'Node {i}', alpha=0.6)

plt.title("Opinion Dynamics in a Social Network (DeGroot Model)")
plt.xlabel("Iteration")
plt.ylabel("Opinion Value")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize='small')
plt.grid(True)
plt.show()

# Plot the initial and final opinions on the graph
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Initial opinions
nx.draw(
    G, 
    with_labels=True, 
    node_color=initial_opinions, 
    cmap=plt.cm.viridis, 
    node_size=500, 
    ax=axes[0]
)
axes[0].set_title("Initial Opinions")

# Final opinions
nx.draw(
    G, 
    with_labels=True, 
    node_color=opinions, 
    cmap=plt.cm.viridis, 
    node_size=500, 
    ax=axes[1]
)
axes[1].set_title("Final Opinions")

plt.show()

Explanation

Graph Creation:
A random graph is created using the $Erdős–Rényi model$ to simulate a social network.
Initialization:
Each node is assigned a random opinion value between $0$ and $1$.
Opinion Update Rule:
Opinions are updated iteratively as a weighted average of neighbors’ opinions, governed by the adjacency matrix of the graph.
Visualization:
- The convergence of opinions over iterations is plotted.
- The initial and final opinions are visualized on the network.

Results

Convergence Dynamics:
- The plot of opinion dynamics shows how individual nodes’ opinions evolve and converge over iterations.
- Nodes with many connections (high degree) influence others more strongly.

Graph Visualization:
- Initial opinions are randomly distributed.
- Final opinions converge to a consensus or remain close to a weighted average of initial opinions depending on the network structure.

This example demonstrates how mathematical sociology models can capture social influence processes using $Python$.

Mathematical Biology

January 14, 2025

Example: Lotka-Volterra Predator-Prey Model

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

This is a classic example in mathematical biology.

Problem Statement

Define the Lotka-Volterra equations for predator-prey dynamics.
Solve the equations numerically using $Python$.
Visualize the population dynamics over time and as a phase plot.

Lotka-Volterra Equations

The equations are:
$$
\frac{dx}{dt} = \alpha x - \beta x y
$$
$$
\frac{dy}{dt} = \delta x y - \gamma y
$$

$x(t)$: Prey population
$y(t)$: Predator population
$\alpha$: Prey birth rate
$\beta$: Predation rate
$\delta$: Predator reproduction rate
$\gamma$: Predator death rate

Python Code

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

# Define the Lotka-Volterra equations
def lotka_volterra(t, z, alpha, beta, delta, gamma):
    x, y = z
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# Parameters
alpha = 0.1   # Prey birth rate
beta = 0.02   # Predation rate
delta = 0.01  # Predator reproduction rate
gamma = 0.1   # Predator death rate

# Initial populations
x0 = 40  # Initial prey population
y0 = 9   # Initial predator population

# Time span for simulation
t_span = (0, 200)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the equations
solution = solve_ivp(
    lotka_volterra, t_span, [x0, y0], args=(alpha, beta, delta, gamma), t_eval=t_eval
)

# Extract results
time = solution.t
prey = solution.y[0]
predators = solution.y[1]

# Plot population dynamics over time
plt.figure(figsize=(10, 6))
plt.plot(time, prey, label="Prey Population", color="blue")
plt.plot(time, predators, label="Predator Population", color="red")
plt.title("Lotka-Volterra Predator-Prey Dynamics")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend()
plt.grid()
plt.show()

# Phase plot
plt.figure(figsize=(8, 6))
plt.plot(prey, predators, color="purple")
plt.title("Phase Plot: Predator vs. Prey")
plt.xlabel("Prey Population")
plt.ylabel("Predator Population")
plt.grid()
plt.show()

Code Explanation

Lotka-Volterra Function:
- Defines the differential equations for the prey and predator populations.
Parameters:
- $ \alpha, \beta, \delta, \gamma $: Control the growth and interaction rates of prey and predators.
Numerical Solution:
- The solve_ivp function from scipy.integrate solves the equations over time.
Visualization:
- Population Dynamics Plot: Shows how prey and predator populations change over time.
- Phase Plot: Illustrates the predator-prey relationship in the phase space.

Results

Ecosystem Analysis Results:

Model Parameters:
Prey growth rate (α): 1.0
Predation rate (β): 0.1
Predator reproduction rate (δ): 0.075
Predator death rate (γ): 0.5

Equilibrium Point:
Prey population: 6.67
Predator population: 10.00

Stability Analysis:
Eigenvalues: 0.000+0.707j, 0.000-0.707j
System is neutrally stable

Population Dynamics:
- Prey population oscillates due to predation.
- Predator population follows the prey’s oscillation with a phase lag.
Phase Plot:
- A closed-loop trajectory represents the cyclic nature of the predator-prey interaction.

This simulation provides insight into the oscillatory behavior of ecological systems and the delicate balance between predator and prey populations.

Financial Mathematics

January 13, 2025

Black-Scholes Option Pricing Model

The Black-Scholes model is a widely used framework in financial mathematics for pricing European options.

In this example, we will calculate the price of a European call option using the Black-Scholes formula.

Problem

Given:

Stock price $(S)$
Strike price $(K)$
Risk-free interest rate $(r)$
Time to maturity $(T)$
Volatility $(\sigma)$ of the stock
Calculate the price of a European call option using the formula:

$$
C = S \cdot N(d_1) - K \cdot e^{-rT} \cdot N(d_2)
$$

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

Here, $N(x)$ is the cumulative distribution function ($CDF$) of the standard normal distribution.

Python Code

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

def black_scholes_call(S, K, r, T, sigma):
    """
    Calculate the price of a European call option using the Black-Scholes formula.
    """
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    call_price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    return call_price

def visualize_black_scholes():
    # Parameters
    S = np.linspace(50, 150, 100)  # Stock price range
    K = 100                       # Strike price
    r = 0.05                      # Risk-free interest rate
    T = 1                         # Time to maturity (1 year)
    sigma = 0.2                   # Volatility (20%)
    
    # Calculate call option prices for each stock price
    call_prices = [black_scholes_call(s, K, r, T, sigma) for s in S]
    
    # Plot the results
    plt.figure(figsize=(10, 6))
    plt.plot(S, call_prices, label="Call Option Price")
    plt.axvline(K, color="red", linestyle="--", label="Strike Price (K)")
    plt.title("European Call Option Price (Black-Scholes Model)")
    plt.xlabel("Stock Price (S)")
    plt.ylabel("Call Option Price")
    plt.legend()
    plt.grid()
    plt.show()

if __name__ == "__main__":
    visualize_black_scholes()

Explanation

Parameters:
- S: The current stock price varies from $50$ to $150$.
- K: The strike price is set at $100.$
- r: The risk-free interest rate is $5$%.
- T: Time to maturity is $1$ year.
- sigma: Volatility is $20$%.
Black-Scholes Formula:
- d1 and d2 are intermediate terms used to calculate the option price.
- norm.cdf(d1) and norm.cdf(d2) give the probabilities of the option being in-the-money.
Visualization:
- The $x$-axis represents the stock price.
- The $y$-axis represents the calculated call option price.
- A red dashed line indicates the strike price $(K = 100)$.

Results

When the stock price $(S)$ is much lower than the strike price $(K)$, the call option price approaches zero because the option is out-of-the-money.
As $(S)$ increases beyond $(K)$, the call option price increases, reflecting the intrinsic value of the option.
The graph provides an intuitive understanding of how the call option price varies with the stock price, according to the Black-Scholes model.

Function Theory

January 12, 2025

I’ll create an example demonstrating function theory using $Python$, focusing on analyzing a complex function with visualization.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def analyze_complex_function(z):
    """
    Analyzes the complex function f(z) = z^2 + 2z + 1
    Args:
        z (complex): Input complex number
    Returns:
        complex: Result of f(z)
    """
    return z**2 + 2*z + 1

# Create a grid of points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = X + Y*1j

# Calculate function values
W = analyze_complex_function(Z)

# Extract real and imaginary parts
real_part = W.real
imag_part = W.imag
abs_values = np.abs(W)

def plot_function_analysis():
    """
    Creates three subplots showing different aspects of the complex function
    """
    fig = plt.figure(figsize=(15, 5))
    
    # Real part
    ax1 = fig.add_subplot(131, projection='3d')
    surf1 = ax1.plot_surface(X, Y, real_part, cmap='viridis')
    ax1.set_title('Real Part')
    ax1.set_xlabel('Re(z)')
    ax1.set_ylabel('Im(z)')
    ax1.set_zlabel('Re(f(z))')
    
    # Imaginary part
    ax2 = fig.add_subplot(132, projection='3d')
    surf2 = ax2.plot_surface(X, Y, imag_part, cmap='viridis')
    ax2.set_title('Imaginary Part')
    ax2.set_xlabel('Re(z)')
    ax2.set_ylabel('Im(z)')
    ax2.set_zlabel('Im(f(z))')
    
    # Absolute value
    ax3 = fig.add_subplot(133, projection='3d')
    surf3 = ax3.plot_surface(X, Y, abs_values, cmap='viridis')
    ax3.set_title('Absolute Value')
    ax3.set_xlabel('Re(z)')
    ax3.set_ylabel('Im(z)')
    ax3.set_zlabel('|f(z)|')
    
    plt.tight_layout()
    return fig

# Plot the results
plot_function_analysis()

# Find critical points (where derivative is zero)
# For f(z) = z^2 + 2z + 1, f'(z) = 2z + 2
# Critical point is at z = -1

critical_point = -1
critical_value = analyze_complex_function(critical_point)
print(f"Critical point: z = {critical_point}")
print(f"Value at critical point: f({critical_point}) = {critical_value}")

Let me explain this example of function theory:

The Complex Function
We’re analyzing the function $f(z) = z² + 2z + 1$, which is a quadratic complex function.
Code Components:
- We define the complex function using $NumPy$ for efficient computation
- Create a grid of complex numbers using meshgrid
- Calculate function values across the complex plane
- Generate $3D$ visualizations showing different aspects of the function
Visualizations:
The code creates three $3D$ plots:
- Real part: Shows how the real component of $f(z)$ varies
- Imaginary part: Shows how the imaginary component varies
- Absolute value: Shows the magnitude of $f(z)$
Mathematical Analysis:
- The function has one critical point at $z = -1$ (where $f’(z) = 0$)
- The function is analytic everywhere (infinitely differentiable)
- It’s a second-degree polynomial in the complex plane

The $3D$ visualizations help understand:

The behavior of the function in different regions of the complex plane
How real and imaginary parts interact
Where the function grows or decreases rapidly
The location and nature of critical points

Environmental Mathematics

January 11, 2025

Example: Modeling Population Growth with Logistic Growth Equation

In Environmental Mathematics, one of the fundamental models used to describe the growth of populations, such as animal populations or human populations, is the Logistic Growth Model.

This model accounts for the carrying capacity of the environment, which limits indefinite exponential growth.

Logistic Growth Equation

The Logistic Growth Model is represented by the following differential equation:

$$
\frac{dP}{dt} = rP \left(1 - \frac{P}{K}\right)
$$

$( P(t) )$ = Population size at time $( t )$
$( r )$ = Intrinsic growth rate (how quickly the population grows)
$( K )$ = Carrying capacity (maximum population size that the environment can support)
$( t )$ = Time
$( P(0) )$ = Initial population size at $( t = 0 )$

The solution to this equation models the population growth over time and shows an S-shaped curve, where the population grows exponentially when small, then slows down as it approaches the carrying capacity $( K )$.

Problem

We will use the Logistic Growth Model to model the population growth of a species over time.
We will visualize how the population changes as a function of time.
We will explore the effect of different intrinsic growth rates and carrying capacities on the population size.

Steps

Define the logistic growth equation and its solution.
Simulate the population growth using $Python$.
Plot the population size over time for different parameters.

Python Code

import numpy as np
import matplotlib.pyplot as plt

# Logistic growth model function
def logistic_growth(t, P0, r, K):
    """
    Computes the population at time t based on the logistic growth model.
    
    Parameters:
    t: time (array or scalar)
    P0: initial population size
    r: intrinsic growth rate
    K: carrying capacity
    
    Returns:
    Population size at time t
    """
    return K / (1 + (K - P0) / P0 * np.exp(-r * t))

# Parameters
P0 = 10       # Initial population size
r = 0.5       # Intrinsic growth rate
K = 1000      # Carrying capacity
t = np.linspace(0, 20, 200)  # Time from 0 to 20 with 200 points

# Compute population over time
population = logistic_growth(t, P0, r, K)

# Plot the logistic growth curve
plt.figure(figsize=(8, 6))
plt.plot(t, population, label=f'Logistic Growth\nP0={P0}, r={r}, K={K}', color='b')

# Labels and title
plt.title('Logistic Growth Model: Population vs Time')
plt.xlabel('Time (t)')
plt.ylabel('Population (P)')
plt.grid(True)
plt.legend()
plt.show()

# Effect of different intrinsic growth rates (r) on population growth
r_values = [0.1, 0.3, 0.5, 0.7]
plt.figure(figsize=(8, 6))

for r in r_values:
    population_r = logistic_growth(t, P0, r, K)
    plt.plot(t, population_r, label=f'r = {r}')

plt.title('Effect of Intrinsic Growth Rate on Logistic Growth')
plt.xlabel('Time (t)')
plt.ylabel('Population (P)')
plt.legend()
plt.grid(True)
plt.show()

# Effect of different carrying capacities (K) on population growth
K_values = [500, 1000, 1500, 2000]
plt.figure(figsize=(8, 6))

for K_val in K_values:
    population_K = logistic_growth(t, P0, r, K_val)
    plt.plot(t, population_K, label=f'K = {K_val}')

plt.title('Effect of Carrying Capacity on Logistic Growth')
plt.xlabel('Time (t)')
plt.ylabel('Population (P)')
plt.legend()
plt.grid(True)
plt.show()

Explanation of the Code

Logistic Growth Function:
- The function logistic_growth(t, P0, r, K) computes the population at time $( t )$ using the logistic growth equation.
  It takes:
  - t: Time (which can be an array for continuous simulation)
  - P0: The initial population size at time $( t = 0 )$
  - r: The intrinsic growth rate of the population
  - K: The carrying capacity (maximum population the environment can support)
- The function calculates the population at time $( t )$ using the formula for logistic growth.
Parameters:
- $( P_0 = 10 )$: Initial population size.
- $( r = 0.5 )$: Intrinsic growth rate (the speed of population growth).
- $( K = 1000 )$: Carrying capacity (the environment can support up to $1000$ individuals).
- t = np.linspace(0, 20, 200): We simulate the population over $20$ time units, with $200$ points to generate a smooth curve.
Plotting the Results:
- The first plot visualizes the population growth over time using the logistic model.
- The second plot shows how different intrinsic growth rates (r = 0.1, 0.3, 0.5, 0.7) affect the population growth.
  A higher growth rate results in faster population growth.
- The third plot shows how varying the carrying capacity (K = 500, 1000, 1500, 2000) influences the population.
  A higher carrying capacity allows the population to grow to a higher final size.

Expected Outcome

Logistic Growth Curve:
The first plot will show the characteristic S-shaped curve.
The population starts growing exponentially, but as it approaches the carrying capacity $( K )$, it slows down and levels off.
Effect of Intrinsic Growth Rate:
In the second plot, as the intrinsic growth rate $( r )$ increases, the population grows faster at the beginning, but the population ultimately approaches the same carrying capacity in each case.
The time it takes to reach near the carrying capacity is shorter for higher growth rates.
Effect of Carrying Capacity:
In the third plot, increasing the carrying capacity $( K )$ leads to a higher maximum population.
However, the growth curve still follows the same logistic pattern, only the “$ceiling$” is higher.

Result Analysis

The S-shaped curve is a key feature of the logistic growth model, showing how the population grows rapidly at first, but as resources become limited, the growth slows and eventually stabilizes near the carrying capacity.
The intrinsic growth rate $( r )$ affects the speed of growth.
A higher growth rate leads to a quicker approach to the carrying capacity.
The carrying capacity $( K )$ defines the maximum population that the environment can support.
If $( K )$ is increased, the population can grow to a higher value before stabilizing.

Conclusion

This example demonstrates how Environmental Mathematics uses mathematical models like logistic growth to understand population dynamics in a given environment.

By adjusting parameters such as the intrinsic growth rate and carrying capacity, we can predict how populations will evolve over time and make informed decisions for environmental management, conservation, and resource allocation.

Tunneling Effect in Physics

January 10, 2025

An Example with Python Code

The Tunneling Effect in physics refers to the phenomenon where a particle has a probability of crossing a potential barrier, even if its energy is lower than the height of the barrier.

This effect is most commonly discussed in quantum mechanics and has significant applications in areas like nuclear fusion, semiconductor physics, and scanning tunneling microscopes.

One way to visualize this phenomenon is to consider a quantum particle in a potential well with a barrier.

Classically, if the particle’s energy is lower than the barrier height, it should not be able to pass through.

However, quantum mechanically, there is a finite probability for the particle to “tunnel” through the barrier.

Example Problem

We will model a particle in a one-dimensional potential well with a rectangular potential barrier.

The wavefunction of the particle is described by the $Schrödinger$ $equation$, and the probability of tunneling can be computed using a simplified approach.

Potential setup:
A rectangular potential barrier with height $( V_0 )$ and width $( L )$.
Particle energy:
$( E )$, where $( E < V_0 )$, so the particle does not have enough energy to overcome the barrier classically.

We can compute the transmission probability $ T(E) $, which gives the likelihood that the particle will tunnel through the barrier.

For a rectangular potential barrier, the transmission probability is given by:

$$
T(E) = e^{-2 \gamma L}
$$

Where $( \gamma )$ is given by:

$$
\gamma = \frac{\sqrt{2m(V_0 - E)}}{\hbar}
$$

$( m )$ is the mass of the particle.
$( V_0 )$ is the potential barrier height.
$( E )$ is the energy of the particle.
$( L )$ is the width of the potential barrier.
$( \hbar )$ is the reduced Planck’s constant.

Python Implementation:

We will write $Python$ code to compute the transmission probability for various particle energies and plot the result.

import numpy as np
import matplotlib.pyplot as plt

# Constants
hbar = 1.0545718e-34  # Reduced Planck's constant in Joules * seconds
m = 9.11e-31          # Mass of the electron in kg (approximate)
V0 = 1.0e-18          # Potential barrier height in Joules (approx 6 eV)
L = 1.0e-10           # Width of the potential barrier in meters (1 Å)

# Function to calculate tunneling transmission probability
def tunneling_probability(E, V0, L, m, hbar):
    if E >= V0:
        return 1  # No tunneling if energy is higher than or equal to barrier height
    gamma = np.sqrt(2 * m * (V0 - E)) / hbar
    return np.exp(-2 * gamma * L)

# Energies of the particle (from 0 to V0)
E_values = np.linspace(0, V0, 500)

# Calculate tunneling probabilities for each energy
T_values = [tunneling_probability(E, V0, L, m, hbar) for E in E_values]

# Plotting the results
plt.figure(figsize=(8, 6))
plt.plot(E_values / 1.60218e-19, T_values, label='Transmission Probability T(E)', color='b')
plt.xlabel("Energy (eV)")
plt.ylabel("Transmission Probability T(E)")
plt.title("Quantum Tunneling Effect: Transmission Probability vs Energy")
plt.grid(True)
plt.axhline(0, color='black',linewidth=1)
plt.axvline(V0 / 1.60218e-19, color='red', linestyle='--', label='V0 (Barrier Height)')
plt.legend()
plt.show()

Code Explanation:

Constants: We define the constants:
- hbar: Reduced Planck’s constant (in Joules·seconds).
- m: The mass of an electron in $kg$.
- V0: The potential barrier height in Joules ($1$ $eV$ is approximately $( 1.60218 \times 10^{-19} )$ J).
- L: The width of the potential barrier (in meters, $1$ $Ångström$).
Tunneling Probability Function:
- The function tunneling_probability(E, V0, L, m, hbar) computes the transmission probability $ T(E) $ using the formula derived above.
- If the particle’s energy $( E )$ is greater than or equal to the barrier height $( V_0 )$, tunneling is not possible, so $( T(E) = 1 )$ (i.e., no tunneling).
- Otherwise, the transmission probability is computed using the exponential decay formula.
Energy Range:
- We generate a range of energy values from $0$ to $( V_0 )$.
  These energies represent possible particle energies below the barrier height.
Plotting:
- The plot shows the transmission probability $ T(E) $ as a function of energy.
- The $x$-axis represents the energy of the particle in electron volts ($eV$).
- The $y$-axis shows the transmission probability.
- We highlight the potential barrier height $( V_0 )$ using a vertical dashed red line.

Graphical Interpretation

As the energy of the particle approaches the barrier height $( V_0 )$, the transmission probability increases.
For energies much lower than the barrier height, the probability decreases exponentially, indicating that tunneling becomes less likely.
At higher energies (but still below $( V_0 )$), the probability of tunneling increases slightly, though it is still low.

Graph Analysis

Transmission at low energies:
For energies significantly lower than the barrier height, tunneling is highly unlikely, and $ T(E) $ is close to $0$.
Near the barrier height:
As the energy gets closer to $( V_0 )$, the tunneling probability increases due to the reduced decay rate $( \gamma )$.
Beyond the barrier height:
If $( E )$ were greater than or equal to $( V_0 )$, the transmission probability would be $1$, indicating that the particle can easily pass through the barrier.

This model demonstrates how quantum tunneling occurs in a system where a particle has insufficient classical energy to overcome a potential barrier but still has a non-zero probability of passing through it due to the quantum nature of particles.

Phase Transitions

January 9, 2025

Here’s an example of a Phase Transition in physics, demonstrated using the Ising Model in $Python$.

The Ising Model is a mathematical model used to describe ferromagnetism in statistical mechanics.

Problem

We want to simulate the 2D Ising Model to understand phase transitions.
Specifically, we will:

Simulate a 2D lattice of spins that can either be $+1$ or $-1$.
Study how the magnetization of the system changes with temperature.
Observe the critical temperature at which a phase transition occurs.

Python Code Implementation

import numpy as np
import matplotlib.pyplot as plt

def initialize_lattice(size):
    """Initialize a lattice with random spins (+1 or -1)."""
    return np.random.choice([-1, 1], size=(size, size))

def calculate_energy(lattice):
    """Calculate the total energy of the lattice."""
    energy = 0
    size = lattice.shape[0]
    for i in range(size):
        for j in range(size):
            spin = lattice[i, j]
            neighbors = lattice[(i + 1) % size, j] + lattice[i, (j + 1) % size] + \
                        lattice[(i - 1) % size, j] + lattice[i, (j - 1) % size]
            energy -= spin * neighbors
    return energy / 2  # Each pair is counted twice

def metropolis_step(lattice, beta):
    """Perform a single Metropolis step."""
    size = lattice.shape[0]
    for _ in range(size ** 2):
        i, j = np.random.randint(0, size, size=2)
        spin = lattice[i, j]
        neighbors = lattice[(i + 1) % size, j] + lattice[i, (j + 1) % size] + \
                    lattice[(i - 1) % size, j] + lattice[i, (j - 1) % size]
        delta_energy = 2 * spin * neighbors
        if delta_energy < 0 or np.random.rand() < np.exp(-beta * delta_energy):
            lattice[i, j] *= -1
    return lattice

def calculate_magnetization(lattice):
    """Calculate the magnetization of the lattice."""
    return np.abs(np.sum(lattice)) / lattice.size

def simulate_ising(size, temperatures, steps):
    """Simulate the Ising model for a range of temperatures."""
    lattice = initialize_lattice(size)
    magnetizations = []
    
    for temp in temperatures:
        beta = 1 / temp
        for _ in range(steps):  # Equilibrate
            lattice = metropolis_step(lattice, beta)
        magnetization = calculate_magnetization(lattice)
        magnetizations.append(magnetization)
    
    return magnetizations

# Parameters
size = 20
temperatures = np.linspace(1.5, 3.5, 20)
steps = 1000

# Simulation
magnetizations = simulate_ising(size, temperatures, steps)

# Plotting Results
plt.figure(figsize=(8, 6))
plt.plot(temperatures, magnetizations, marker='o', color='b', label="Magnetization")
plt.axvline(x=2.27, color='r', linestyle='--', label="Critical Temperature (Tc)")
plt.title("Phase Transition in 2D Ising Model")
plt.xlabel("Temperature (T)")
plt.ylabel("Magnetization")
plt.legend()
plt.grid()
plt.show()

Explanation

Initialization:
The lattice is a grid of spins, initialized randomly.
Metropolis Algorithm:
For each spin, we calculate the energy change (ΔE) if the spin is flipped.
The flip is accepted based on the Metropolis criterion, which ensures proper statistical weighting.
Magnetization:
The total magnetization of the system is recorded as the absolute sum of spins normalized by the number of spins.
Temperature Sweep:
By sweeping through a range of temperatures, we observe how the magnetization changes, showing a sharp drop near the critical temperature $(( T_c \approx 2.27 )$ for the 2D Ising Model).

Result

The graph shows magnetization as a function of temperature.

Below the critical temperature $(T_c)$, the system exhibits spontaneous magnetization, indicating a ferromagnetic phase.

Above $(T_c)$, the magnetization drops to zero, showing a transition to the paramagnetic phase.

The red vertical line indicates $(T_c)$.