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)$.

Entropy Changes

January 8, 2025

Below is an example problem related to Entropy Changes, its solution in $Python$, and a detailed explanation.

Problem: Entropy Changes in Mixing Two Substances

Two substances, A and B, with masses $( m_A = 1 \text{kg} )$ and $( m_B = 2 \text{kg} )$, are mixed.

Both substances are water with a specific heat capacity $( c = 4184 \text{J/(kg·K)} )$.

Initial temperatures: $( T_A = 80^\circ \text{C} )$ and $( T_B = 20^\circ \text{C} )$.
Final equilibrium temperature $( T_f )$ is reached by energy conservation.

Calculate the total entropy change $( \Delta S )$ of the system during this process.

Solution Approach

Use energy conservation to calculate the final equilibrium temperature $( T_f )$:
$$
m_A c (T_A - T_f) = m_B c (T_f - T_B)
$$
Calculate the entropy change for each substance:
$$
\Delta S = \int_{T_i}^{T_f} \frac{dQ}{T}
$$
Since $( dQ = mc , dT )$:
$$
\Delta S = mc \ln\left(\frac{T_f}{T_i}\right)
$$
Sum up $( \Delta S_A )$ and $( \Delta S_B )$ for the total entropy change.

Python Implementation

import numpy as np
import matplotlib.pyplot as plt

# Constants
c = 4184  # Specific heat capacity (J/(kg·K))
m_A = 1.0  # Mass of substance A (kg)
m_B = 2.0  # Mass of substance B (kg)
T_A = 80 + 273.15  # Initial temperature of A (K)
T_B = 20 + 273.15  # Initial temperature of B (K)

# Calculate final equilibrium temperature (energy conservation)
T_f = (m_A * c * T_A + m_B * c * T_B) / (m_A * c + m_B * c)

# Entropy change calculations
def entropy_change(m, c, T_i, T_f):
    return m * c * np.log(T_f / T_i)

delta_S_A = entropy_change(m_A, c, T_A, T_f)
delta_S_B = entropy_change(m_B, c, T_B, T_f)
total_entropy_change = delta_S_A + delta_S_B

# Display results
print(f"Final equilibrium temperature (T_f): {T_f - 273.15:.2f} °C")
print(f"Entropy change for substance A (ΔS_A): {delta_S_A:.2f} J/K")
print(f"Entropy change for substance B (ΔS_B): {delta_S_B:.2f} J/K")
print(f"Total entropy change (ΔS_total): {total_entropy_change:.2f} J/K")

# Visualization
temperatures = np.linspace(T_B, T_A, 100)  # Range of temperatures
entropy_A = m_A * c * np.log(T_f / temperatures)
entropy_B = m_B * c * np.log(T_f / temperatures[::-1])

plt.figure(figsize=(8, 6))
plt.plot(temperatures - 273.15, entropy_A, label="Entropy Change (A)")
plt.plot(temperatures - 273.15, entropy_B, label="Entropy Change (B)")
plt.axhline(y=total_entropy_change, color='r', linestyle='--', label="Total ΔS")
plt.title("Entropy Changes During Mixing")
plt.xlabel("Temperature (°C)")
plt.ylabel("Entropy Change (J/K)")
plt.legend()
plt.grid()
plt.show()

Explanation of the Code

Constants and Inputs:
- The specific heat capacity $( c )$ is set for water.
- Initial temperatures are converted to Kelvin for calculations.
Final Temperature Calculation:
- Using energy conservation, $( T_f )$ is calculated.
Entropy Change Function:
- $ \Delta S = mc \ln(T_f / T_i) $ is implemented for both substances.
Results:
- The entropy change for substances A and B and the total entropy change are displayed.
Graph:
- The graph shows how the entropy of A and B changes with temperature, along with the total entropy change.

Results Interpretation

Final equilibrium temperature (T_f): 40.00 °C
Entropy change for substance A (ΔS_A): -502.96 J/K
Entropy change for substance B (ΔS_B): 552.27 J/K
Total entropy change (ΔS_total): 49.31 J/K

Final Temperature: The equilibrium temperature is approximately $( 40^\circ \text{C} )$.
Entropy Changes: Both substances experience positive entropy changes, as the system moves toward equilibrium.
Total Entropy Change: The total entropy change is positive, consistent with the second law of thermodynamics.

Wave Propagation

January 7, 2025

I’ll help you create a simulation of wave propagation, including both traveling waves and standing waves.

We’ll simulate and visualize the wave equation.

Wave Equation:
$$
\frac{\partial^2 y}{\partial t^2} = c^2\frac{\partial^2 y}{\partial x^2}
$$

$y(x,t)$ is the displacement
$c$ is the wave speed
$x$ is position
$t$ is time

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

class WaveSimulator:
    def __init__(self, length=10, nx=1000, dt=0.01, c=1):
        """
        Initialize wave simulation parameters
        length: length of medium
        nx: number of spatial points
        dt: time step
        c: wave speed
        """
        self.length = length
        self.nx = nx
        self.dt = dt
        self.c = c
        
        self.x = np.linspace(0, length, nx)
        self.dx = self.x[1] - self.x[0]
        
        # Check stability condition (CFL condition)
        if c * dt / self.dx > 1:
            raise ValueError("Stability condition not met. Decrease dt or increase dx.")
            
    def traveling_wave(self, t, wavelength, amplitude=1):
        """Generate traveling wave"""
        k = 2 * np.pi / wavelength  # wave number
        omega = self.c * k          # angular frequency
        return amplitude * np.sin(k * self.x - omega * t)
    
    def standing_wave(self, t, mode=1, amplitude=1):
        """Generate standing wave"""
        k = mode * np.pi / self.length
        omega = self.c * k
        return amplitude * np.sin(k * self.x) * np.cos(omega * t)
    
    def wave_packet(self, t, x0=0, sigma=1, k0=5):
        """Generate Gaussian wave packet"""
        x = self.x - (x0 + self.c * t)
        return np.exp(-(x**2)/(2*sigma**2)) * np.cos(k0 * x)

# Create simulation instance
sim = WaveSimulator(length=10, nx=1000, dt=0.01, c=1)

# Time points for simulation
t = np.linspace(0, 10, 100)
X, T = np.meshgrid(sim.x, t)

# Generate different types of waves
traveling_wave = np.array([sim.traveling_wave(ti, wavelength=2) for ti in t])
standing_wave = np.array([sim.standing_wave(ti, mode=2) for ti in t])
wave_packet = np.array([sim.wave_packet(ti, x0=0, sigma=1) for ti in t])

# Create visualizations
fig = plt.figure(figsize=(15, 12))

# 3D plot of traveling wave
ax1 = fig.add_subplot(311, projection='3d')
ax1.plot_surface(X, T, traveling_wave, cmap='viridis')
ax1.set_xlabel('Position (x)')
ax1.set_ylabel('Time (t)')
ax1.set_zlabel('Amplitude')
ax1.set_title('Traveling Wave')

# 3D plot of standing wave
ax2 = fig.add_subplot(312, projection='3d')
ax2.plot_surface(X, T, standing_wave, cmap='viridis')
ax2.set_xlabel('Position (x)')
ax2.set_ylabel('Time (t)')
ax2.set_zlabel('Amplitude')
ax2.set_title('Standing Wave (Mode 2)')

# 3D plot of wave packet
ax3 = fig.add_subplot(313, projection='3d')
ax3.plot_surface(X, T, wave_packet, cmap='viridis')
ax3.set_xlabel('Position (x)')
ax3.set_ylabel('Time (t)')
ax3.set_zlabel('Amplitude')
ax3.set_title('Gaussian Wave Packet')

plt.tight_layout()
plt.show()

# Create animation of wave propagation
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12))

# Initialize lines
line1, = ax1.plot([], [], 'b-', label='Traveling Wave')
line2, = ax2.plot([], [], 'r-', label='Standing Wave')
line3, = ax3.plot([], [], 'g-', label='Wave Packet')

# Set axis limits
for ax in [ax1, ax2, ax3]:
    ax.set_xlim(0, sim.length)
    ax.set_ylim(-1.5, 1.5)
    ax.grid(True)
    ax.legend()

ax1.set_title('Traveling Wave')
ax2.set_title('Standing Wave')
ax3.set_title('Wave Packet')

# Animation function
def animate(frame):
    line1.set_data(sim.x, traveling_wave[frame])
    line2.set_data(sim.x, standing_wave[frame])
    line3.set_data(sim.x, wave_packet[frame])
    return line1, line2, line3

# Create animation
anim = FuncAnimation(fig, animate, frames=len(t), interval=50, blit=True)
plt.tight_layout()
plt.show()

Let me explain the key components of this simulation:

Wave Types Implemented:

a) Traveling Wave:
- Form: $y(x,t) = A\sin(kx - \omega t)$
- Moves with constant speed in one direction
- Maintains shape while propagating
b) Standing Wave:
- Form: $y(x,t) = A\sin(kx)\cos(\omega t)$
- Fixed nodes and antinodes
- Result of interference between waves traveling in opposite directions
c) Gaussian Wave Packet:
- Localized wave group
- Combines carrier wave with Gaussian envelope
- Shows dispersion effects

Pendulum Dynamics

January 6, 2025

I’ll help you create a $Python$ program to simulate and visualize the dynamics of a pendulum, including both simple and damped pendulum cases.

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

class PendulumSimulator:
    def __init__(self, L=1.0, g=9.81, m=1.0, b=0.0):
        """
        Initialize pendulum parameters
        L: length of pendulum (m)
        g: gravitational acceleration (m/s^2)
        m: mass of bob (kg)
        b: damping coefficient (kg/s)
        """
        self.L = L
        self.g = g
        self.m = m
        self.b = b

    def derivatives(self, state, t):
        """
        Calculate derivatives for the pendulum system
        state[0] = theta (angle)
        state[1] = omega (angular velocity)
        """
        theta, omega = state
        
        # Calculate acceleration (second derivative of theta)
        # Include nonlinear pendulum dynamics and damping
        theta_dot = omega
        omega_dot = (-self.b * omega / self.m 
                    - (self.g / self.L) * np.sin(theta))
        
        return [theta_dot, omega_dot]

    def simulate(self, theta0, omega0, t_span):
        """
        Simulate pendulum motion
        theta0: initial angle (radians)
        omega0: initial angular velocity (rad/s)
        t_span: time points for simulation
        """
        initial_state = [theta0, omega0]
        solution = odeint(self.derivatives, initial_state, t_span)
        return solution

# Create simulation instances
simple_pendulum = PendulumSimulator(L=1.0, b=0.0)  # No damping
damped_pendulum = PendulumSimulator(L=1.0, b=0.5)  # With damping

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

# Initial conditions
theta0 = np.pi/3  # 60 degrees
omega0 = 0.0      # starting from rest

# Run simulations
simple_solution = simple_pendulum.simulate(theta0, omega0, t)
damped_solution = damped_pendulum.simulate(theta0, omega0, t)

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

# Plot angular displacement
ax1 = fig.add_subplot(221)
ax1.plot(t, simple_solution[:, 0], 'b-', label='Simple Pendulum')
ax1.plot(t, damped_solution[:, 0], 'r--', label='Damped Pendulum')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Angle (radians)')
ax1.set_title('Angular Displacement vs Time')
ax1.grid(True)
ax1.legend()

# Plot angular velocity
ax2 = fig.add_subplot(222)
ax2.plot(t, simple_solution[:, 1], 'b-', label='Simple Pendulum')
ax2.plot(t, damped_solution[:, 1], 'r--', label='Damped Pendulum')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Angular Velocity (rad/s)')
ax2.set_title('Angular Velocity vs Time')
ax2.grid(True)
ax2.legend()

# Phase space plot
ax3 = fig.add_subplot(223)
ax3.plot(simple_solution[:, 0], simple_solution[:, 1], 'b-', label='Simple Pendulum')
ax3.plot(damped_solution[:, 0], damped_solution[:, 1], 'r--', label='Damped Pendulum')
ax3.set_xlabel('Angle (radians)')
ax3.set_ylabel('Angular Velocity (rad/s)')
ax3.set_title('Phase Space')
ax3.grid(True)
ax3.legend()

# Energy plot
def calculate_energy(solution, pendulum):
    theta, omega = solution[:, 0], solution[:, 1]
    KE = 0.5 * pendulum.m * (pendulum.L * omega)**2
    PE = pendulum.m * pendulum.g * pendulum.L * (1 - np.cos(theta))
    return KE + PE

ax4 = fig.add_subplot(224)
E_simple = calculate_energy(simple_solution, simple_pendulum)
E_damped = calculate_energy(damped_solution, damped_pendulum)
ax4.plot(t, E_simple, 'b-', label='Simple Pendulum')
ax4.plot(t, E_damped, 'r--', label='Damped Pendulum')
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('Total Energy (J)')
ax4.set_title('Total Energy vs Time')
ax4.grid(True)
ax4.legend()

plt.tight_layout()
plt.show()

Let me explain the code and the physics:

Physical System:
- The pendulum consists of a mass m at the end of a rigid, massless rod of length $L$
- Motion is governed by gravity ($g$) and optional damping ($b$)
- The equation of motion is: $θ̈ + (b/m)θ̇ + (g/L)sin(θ) = 0$
Code Structure:
- PendulumSimulator class handles the physics calculations
- Uses scipy.integrate.odeint for numerical integration of the equations of motion
- Simulates both simple (undamped) and damped pendulum cases
Visualizations (4 plots):

a) Angular Displacement vs Time:
- Shows how the angle changes over time
- Simple pendulum oscillates continuously
- Damped pendulum’s oscillations decay

b) Angular Velocity vs Time:
- Shows how the angular velocity changes
- Simple pendulum maintains amplitude
- Damped pendulum shows decreasing velocity magnitude

c) Phase Space Plot:
- Shows relationship between angle and angular velocity
- Simple pendulum forms closed orbit (conservation of energy)
- Damped pendulum spirals toward center (energy dissipation)

d) Total Energy vs Time:
- Shows sum of kinetic and potential energy
- Simple pendulum maintains constant energy
- Damped pendulum shows energy decay

Parameters Used:
- Length ($L$) = $1.0$ m
- Initial angle = $π/3$ radians ($60$ degrees)
- Initial angular velocity = $0$
- Damping coefficient (b) = $0.5$ kg/s for damped case

The simulation shows key physical principles:

Conservation of energy in the simple pendulum
Energy dissipation in the damped pendulum
Nonlinear behavior at large angles
Phase space representation of the dynamics

Magnetic Field of a Current Loop

January 5, 2025

I’ll help you create a $Python$ program to calculate and visualize the magnetic field of a current loop.

We’ll solve a specific example and create visualizations to better understand the field distribution.

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

def magnetic_field_loop(x, y, z, R, I):
    """
    Calculate the magnetic field components (Bx, By, Bz) at point (x, y, z)
    for a circular current loop of radius R carrying current I
    centered at the origin in the x-y plane.
    
    Parameters:
    x, y, z: coordinates of the point where B-field is calculated (meters)
    R: radius of the current loop (meters)
    I: current in the loop (amperes)
    
    Returns:
    tuple: (Bx, By, Bz) magnetic field components in Tesla
    """
    mu0 = 4 * np.pi * 1e-7  # permeability of free space
    
    # Calculate intermediate values
    r = np.sqrt(x**2 + y**2 + z**2)
    if r == 0:
        return (0, 0, 0)
    
    # Calculate components
    factor = mu0 * I * R**2 / (4 * np.pi * (R**2 + r**2)**(3/2))
    
    Bx = factor * 3 * x * z / r**2
    By = factor * 3 * y * z / r**2
    Bz = factor * (3 * z**2 / r**2 - 1)
    
    return (Bx, By, Bz)

# Example parameters
R = 0.1  # radius of loop = 10 cm
I = 1.0  # current = 1 A

# Create grid of points
x = np.linspace(-0.2, 0.2, 20)
y = np.linspace(-0.2, 0.2, 20)
z = np.linspace(-0.2, 0.2, 20)
X, Y, Z = np.meshgrid(x, y, z)

# Calculate B-field components
Bx = np.zeros_like(X)
By = np.zeros_like(Y)
Bz = np.zeros_like(Z)

for i in range(len(x)):
    for j in range(len(y)):
        for k in range(len(z)):
            Bx[i,j,k], By[i,j,k], Bz[i,j,k] = magnetic_field_loop(X[i,j,k], Y[i,j,k], Z[i,j,k], R, I)

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

# XZ-plane plot (y = 0)
ax1 = fig.add_subplot(131)
y_mid = len(y)//2
magnitude = np.sqrt(Bx[:,y_mid,:]**2 + Bz[:,y_mid,:]**2)
ax1.streamplot(x, z, Bx[:,y_mid,:], Bz[:,y_mid,:], 
              color=magnitude, cmap='viridis')
ax1.set_title('Magnetic field lines in XZ-plane (y = 0)')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('z (m)')
ax1.add_patch(plt.Circle((0,0), R, fill=False, color='red'))

# Field strength along z-axis
ax2 = fig.add_subplot(132)
z_axis = np.linspace(-0.2, 0.2, 100)
Bz_axis = [magnetic_field_loop(0, 0, z_val, R, I)[2] for z_val in z_axis]
ax2.plot(z_axis, Bz_axis)
ax2.set_title('Bz along z-axis')
ax2.set_xlabel('z (m)')
ax2.set_ylabel('Bz (T)')
ax2.grid(True)

# 3D vector plot
ax3 = fig.add_subplot(133, projection='3d')
skip = 3
ax3.quiver(X[::skip,::skip,::skip], 
           Y[::skip,::skip,::skip], 
           Z[::skip,::skip,::skip],
           Bx[::skip,::skip,::skip], 
           By[::skip,::skip,::skip], 
           Bz[::skip,::skip,::skip],
           length=0.05, normalize=True)
ax3.set_title('3D magnetic field vectors')
ax3.set_xlabel('x (m)')
ax3.set_ylabel('y (m)')
ax3.set_zlabel('z (m)')

plt.tight_layout()
plt.show()

Let me explain the code and the physics behind it:

Physics Background:
- A circular current loop generates a magnetic field described by the $Biot$-$Savart$ $law$
- The field has axial symmetry around the $z$-axis
- The field is strongest near the loop and decreases with distance
- Inside the loop, the field is mostly uniform and parallel to the $z$-axis
Code Structure:
- The magnetic_field_loop function calculates the magnetic field components (Bx, By, Bz) at any point in space
- We use the analytical solution for the magnetic field of a current loop
- The calculation includes the permeability of free space (μ₀)
- We create a $3D$ grid of points to evaluate the field
Visualizations:
- Left plot: Magnetic field lines in the $xz$-plane ($y = 0$)
  - The red circle represents the current loop
  - The streamlines show the direction of the magnetic field
  - Colors indicate field strength
- Middle plot: Bz component along the $z$-axis
  - Shows how the field strength varies along the axis of the loop
  - The field is strongest near the loop and decreases with distance
- Right plot: $3D$ vector representation
  - Arrows show the direction of the magnetic field in $3D$ space
  - The length of the arrows is normalized for better visualization
Parameters Used:
- Loop radius ($R$) = $0.1$ meters ($10$ cm)
- Current ($I$) = $1$ ampere
- Calculated field is in Tesla units

The visualizations show that:

The field is symmetric around the $z$-axis
Field lines emerge from one side of the loop and curve around to enter the other side
The field is strongest near the loop and decreases with distance
Inside the loop, the field is relatively uniform and points along the $z$-axis

Spring-Mass System

January 4, 2025

I’ll create a simulation of a Spring-Mass System that demonstrates both simple harmonic motion and damped oscillation with external forcing.

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

def spring_mass_system(state, t, m, k, c, F0, omega_f):
    """
    Define the system of first-order ODEs for the spring-mass system
    state[0] = x (position)
    state[1] = v (velocity)
    """
    x, v = state
    
    # External forcing function
    F_external = F0 * np.sin(omega_f * t)
    
    # Acceleration from Newton's second law
    acceleration = (-k * x - c * v + F_external) / m
    
    return [v, acceleration]

# System parameters
m = 1.0  # mass (kg)
k = 10.0  # spring constant (N/m)
c = 0.5   # damping coefficient (Ns/m)
F0 = 1.0  # forcing amplitude (N)
omega_n = np.sqrt(k/m)  # natural frequency
omega_f = omega_n  # forcing frequency (at resonance)

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

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

# Solve ODEs for three cases
# 1. Simple harmonic motion (no damping, no forcing)
solution_shm = odeint(spring_mass_system, state0, t, args=(m, k, 0, 0, omega_f))

# 2. Damped oscillation (with damping, no forcing)
solution_damped = odeint(spring_mass_system, state0, t, args=(m, k, c, 0, omega_f))

# 3. Forced oscillation (with damping and forcing)
solution_forced = odeint(spring_mass_system, state0, t, args=(m, k, c, F0, omega_f))

# Create plots
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 12))

# Plot 1: Simple Harmonic Motion
ax1.plot(t, solution_shm[:, 0], 'b-', label='Position')
ax1.plot(t, solution_shm[:, 1], 'r--', label='Velocity')
ax1.set_title('Simple Harmonic Motion')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Position (m) / Velocity (m/s)')
ax1.grid(True)
ax1.legend()

# Plot 2: Damped Oscillation
ax2.plot(t, solution_damped[:, 0], 'b-', label='Position')
ax2.plot(t, solution_damped[:, 1], 'r--', label='Velocity')
ax2.set_title('Damped Oscillation')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Position (m) / Velocity (m/s)')
ax2.grid(True)
ax2.legend()

# Plot 3: Forced Oscillation
ax3.plot(t, solution_forced[:, 0], 'b-', label='Position')
ax3.plot(t, solution_forced[:, 1], 'r--', label='Velocity')
ax3.set_title('Forced Oscillation with Damping')
ax3.set_xlabel('Time (s)')
ax3.set_ylabel('Position (m) / Velocity (m/s)')
ax3.grid(True)
ax3.legend()

# Calculate and display system parameters
print(f"Natural frequency: {omega_n:.2f} rad/s")
print(f"Damping ratio: {c/(2*np.sqrt(m*k)):.2f}")

plt.tight_layout()

Let me explain the key components of this simulation:

System Parameters:
- Mass $(m) = 1.0 kg$
- Spring constant $(k) = 10.0 N/m$
- Damping coefficient $(c) = 0.5 Ns/m$
- Natural frequency $(ωn) = √(k/m) ≈ 3.16 rad/s$
- Initial displacement = $1.0 m$
- Initial velocity = $0.0 m/s$
Three Cases Simulated:
- Simple Harmonic Motion (no damping, no forcing)
- Damped Oscillation (with damping, no forcing)
- Forced Oscillation (with damping and external force)
Mathematical Model:
- Based on Newton’s Second Law: F = ma
- Differential equation: $mẍ + cẋ + kx = F₀sin(ωt)$
- Solved using scipy’s odeint function
Plot Analysis:
- Top Plot (Simple Harmonic Motion):
  - Shows perfect sinusoidal oscillation
  - Energy continuously transfers between potential and kinetic
- Middle Plot (Damped Oscillation):
  - Amplitude decreases exponentially due to damping
  - System eventually comes to rest
- Bottom Plot (Forced Oscillation):
  - Initial transient response followed by steady-state oscillation
  - Shows resonance effects as forcing frequency matches natural frequency
Physical Insights:
- Energy conservation in simple harmonic motion
- Energy dissipation in damped systems
- Resonance phenomenon in forced oscillation
- Phase relationships between position and velocity

Result

Natural frequency: 3.16 rad/s
Damping ratio: 0.08