The double-slit experiment is a classic demonstration of the wave-particle duality of light and matter.

Below, I’ll present an example simulating the interference pattern observed in the experiment, solved using $Python$.

Problem Description

Simulate the intensity distribution of light on a screen in the double-slit experiment.

Assume:

1. The wavelength of light $( \lambda )$ is $500$ nm.
2. The distance between the slits $( d )$ is $0.1$ mm.
3. The screen is $1$ m away from the slits.
4. The slits emit coherent light.

Python Implementation

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

import numpy as np
import matplotlib.pyplot as plt

# Constants
wavelength = 500e-9  # Wavelength of light (m)
slit_distance = 0.1e-3  # Distance between slits (m)
screen_distance = 1.0  # Distance to the screen (m)
screen_width = 0.02  # Width of the screen (m)
num_points = 1000  # Number of points on the screen

# Create a range of screen positions
x = np.linspace(-screen_width / 2, screen_width / 2, num_points)

# Calculate the path difference
path_difference = (slit_distance * x) / screen_distance

# Calculate the phase difference
phase_difference = (2 * np.pi / wavelength) * path_difference

# Calculate the intensity
intensity = (np.cos(phase_difference / 2)) ** 2

# Normalize intensity
intensity /= np.max(intensity)

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(x * 1e3, intensity, color='blue', label='Interference Pattern')
plt.title("Double-Slit Experiment Simulation")
plt.xlabel("Position on Screen (mm)")
plt.ylabel("Normalized Intensity")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants: Define physical constants such as the wavelength of light, the distance between the slits, and the screen distance.
2. Screen Positions: Create an array of screen positions (x) to simulate the continuous nature of the screen.
3. Path Difference: Calculate the path difference for light reaching each point on the screen.
4. Phase Difference: Compute the corresponding phase difference.
5. Intensity Calculation: Use the formula $( I = \cos^2(\Delta \phi / 2) )$, normalized to a maximum of $1$.
6. Plotting: Plot the intensity distribution across the screen.

Results and Visualization

The graph shows the intensity distribution as bright and dark fringes.

Bright fringes occur where the waves constructively interfere, while dark fringes correspond to destructive interference.
The central maximum is the brightest fringe, gradually diminishing for higher-order fringes.

The quantum harmonic oscillator is a foundational problem in quantum mechanics.

The time-independent Schrödinger equation for this system is:

$$
-\frac{\hbar^2}{2m} \frac{d^2\psi(x)}{dx^2} + \frac{1}{2} m \omega^2 x^2 \psi(x) = E \psi(x)
$$

• $( \hbar )$: Reduced Planck’s constant
• $( m )$: Mass of the particle
• $( \omega )$: Angular frequency
• $( \psi(x) )$: Wavefunction
• $( E )$: Energy eigenvalue

We’ll solve this equation numerically using the finite difference method and plot the wavefunctions $ \psi(x) $ for the first few energy levels.

Python Code

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

# Constants
hbar = 1.0545718e-34  # Reduced Planck's constant (J·s)
m = 9.10938356e-31  # Mass of electron (kg)
omega = 1e15  # Angular frequency (rad/s)

# Grid setup
x_max = 1e-9  # Maximum x value (meters)
n_points = 1000  # Number of grid points
x = np.linspace(-x_max, x_max, n_points)
dx = x[1] - x[0]  # Spacing

# Kinetic energy operator (finite difference approximation)
kinetic = -(hbar**2 / (2 * m * dx**2)) * (
    np.diag(np.ones(n_points - 1), -1) - 2 * np.diag(np.ones(n_points), 0) + np.diag(np.ones(n_points - 1), 1)
)

# Potential energy operator
potential = 0.5 * m * omega**2 * np.diag(x**2)

# Hamiltonian
hamiltonian = kinetic + potential

# Solving for eigenvalues and eigenvectors
eigenvalues, eigenvectors = eigh(hamiltonian)

# Convert eigenvalues to energy levels (in eV)
eigenvalues_eV = eigenvalues / 1.60218e-19

# Plotting the wavefunctions
plt.figure(figsize=(10, 6))
for n in range(5):  # Plot first 5 wavefunctions
    plt.plot(x * 1e9, eigenvectors[:, n]**2, label=f"n={n}, E={eigenvalues_eV[n]:.2f} eV")
plt.title("Wavefunctions of Quantum Harmonic Oscillator", fontsize=14)
plt.xlabel("x (nm)", fontsize=12)
plt.ylabel("Probability Density", fontsize=12)
plt.legend(fontsize=10)
plt.grid(alpha=0.4)
plt.show()

Explanation of the Code

1. Grid Setup:

   • A finite grid is used to represent the position $( x )$ over a symmetric range.
   • The number of points $( n_{points} )$ determines the resolution.

2. Hamiltonian:

   • The kinetic energy is approximated using a finite difference method.
   • The potential energy is represented as $( \frac{1}{2} m \omega^2 x^2 )$.

3. Solving Eigenvalue Problem:

   • The eigenvalues correspond to energy levels.
   • The eigenvectors are the wavefunctions $( \psi(x) )$.

4. Wavefunction Plot:

   • The first $5$ energy levels’ wavefunctions are plotted as probability densities $( |\psi(x)|^2 )$.

Let’s execute the code to see the results.

Graph Explanation

1. Wavefunctions:

   • The plotted curves represent the probability densities ($(|\psi(x)|^2)$) of the quantum harmonic oscillator for the first $5$ energy levels ($n = 0$ to $n = 4$).
   • Each wavefunction corresponds to a distinct energy eigenvalue $(E_n)$, labeled in electron volts (eV).

2. Nodes:

   • Higher energy levels $(n > 0)$ have more nodes (points where $(|\psi(x)|^2 = 0)$) in their wavefunctions, a characteristic feature of quantum systems.

3. Bound States:

   • The wavefunctions are localized, indicating bound states of the particle within the potential well.

This result aligns with theoretical expectations and provides a clear visualization of quantum mechanical states in a harmonic potential.

The motion of a test particle around a black hole can be modeled using the $Schwarzschild$ metric.

The equations governing the orbit are derived from the effective potential:

$$
V_{\text{eff}}(r) = -\frac{GM}{r} + \frac{L^2}{2r^2} - \frac{GM L^2}{c^2 r^3}
$$

• $( G )$: Gravitational constant
• $( M )$: Mass of the black hole
• $( c )$: Speed of light
• $( L )$: Angular momentum of the test particle
• $( r )$: Radial distance from the black hole

We’ll calculate and visualize the effective potential and simulate the orbit of a particle in the $Schwarzschild$ spacetime.

Python Code and Explanation

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

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

# Constants
G = 6.67430e-11  # Gravitational constant, m^3 kg^-1 s^-2
M = 1.989e30  # Mass of black hole (solar mass), kg
c = 3e8  # Speed of light, m/s
L = 4e10  # Angular momentum, m^2/s

# Effective potential function
def effective_potential(r):
    term1 = -G * M / r
    term2 = L**2 / (2 * r**2)
    term3 = -G * M * L**2 / (c**2 * r**3)
    return term1 + term2 + term3

# Radial range
r = np.linspace(2 * G * M / c**2, 1e8, 1000)

# Plot effective potential
plt.figure(figsize=(8, 5))
plt.plot(r, effective_potential(r), label="Effective Potential", color="blue")
plt.axhline(0, color="black", linestyle="--", linewidth=0.8)
plt.title("Effective Potential for a Black Hole", fontsize=14)
plt.xlabel("Radial Distance $r$ (m)", fontsize=12)
plt.ylabel("Potential $V_{eff}$ (J/kg)", fontsize=12)
plt.legend(fontsize=12)
plt.grid(alpha=0.4)
plt.show()

# Simulating orbit using Schwarzschild metric
def schwarzschild_orbit(t, y):
    r, phi, drdt, dphidt = y
    d2rdt2 = L**2 / r**3 - G * M / r**2 + 3 * G * M * L**2 / (c**2 * r**4)
    d2phidt2 = 0  # Angular momentum is conserved
    return [drdt, dphidt, d2rdt2, 0]

# Initial conditions: r0, phi0, dr0/dt, dphi0/dt
r0 = 10 * G * M / c**2  # Start 10 times the Schwarzschild radius
phi0 = 0
dr0_dt = 0  # Initial radial velocity
dphi0_dt = L / r0**2  # Initial angular velocity

initial_conditions = [r0, phi0, dr0_dt, dphi0_dt]

# Time range
t_span = (0, 1e4)  # seconds
t_eval = np.linspace(*t_span, 1000)

# Solve differential equations
solution = solve_ivp(schwarzschild_orbit, t_span, initial_conditions, t_eval=t_eval)

# Extract results
r_sol = solution.y[0]
phi_sol = solution.y[1]

# Convert to Cartesian coordinates for plotting
x = r_sol * np.cos(phi_sol)
y = r_sol * np.sin(phi_sol)

# Plot orbit
plt.figure(figsize=(8, 8))
plt.plot(x, y, label="Orbit", color="red")
plt.scatter(0, 0, color="black", label="Black Hole")  # Black hole at the origin
plt.title("Orbit of a Test Particle Around a Black Hole", fontsize=14)
plt.xlabel("x (m)", fontsize=12)
plt.ylabel("y (m)", fontsize=12)
plt.legend(fontsize=12)
plt.axis("equal")
plt.grid(alpha=0.4)
plt.show()

Explanation

1. Effective Potential:

   • The plot shows how the effective potential depends on radial distance $( r )$.
     This determines stable and unstable orbits.
   • The $Schwarzschild$ radius $( r_s = 2GM/c^2 )$ acts as a boundary; inside it, no stable orbits exist.

2. Orbit Simulation:

   • We solve the equations of motion using the $Schwarzschild$ metric.
   • The orbit is plotted in Cartesian coordinates, showing the trajectory of the test particle.

Let’s execute the code to visualize the results.

Graphs Interpretation

1. Effective Potential:

   • The first graph shows the effective potential $( V_{\text{eff}} )$ as a function of the radial distance $( r )$.
   • At smaller $( r )$, the potential increases rapidly due to the $Schwarzschild$ radius, showing the unstable nature of orbits near the black hole.

2. Orbit of a Test Particle:

   • The second graph illustrates the path of a test particle orbiting a black hole.
   • The black hole is at the origin, and the particle’s trajectory demonstrates its curved path due to the strong gravitational field.
   • This orbit represents a bound system where the angular momentum and gravitational forces balance.

These visualizations help us understand the dynamics of particles in the curved spacetime around a $Schwarzschild$ black hole.

We will simulate the effects of monetary policy on an economy using a simplified model.

The main variables are:

1. Interest Rate (r):
   Set by the central bank.
2. Inflation Rate (π):
   Determined by economic conditions and influenced by monetary policy.
3. Output Gap (Y):
   The difference between actual output and potential output, affected by interest rate changes.

This simulation will follow the logic of the Taylor Rule, a common monetary policy guideline:
$$
r = r^* + \phi_\pi (\pi - \pi^*) + \phi_Y Y
$$

• $(r^*)$ is the neutral interest rate.
• $(\pi^*)$ is the target inflation rate.
• $(\phi_\pi) and (\phi_Y)$ are policy reaction coefficients for inflation and output gap.

Python Implementation

Below is the $Python$ code to simulate the impact of monetary policy adjustments on inflation and the output gap over time.

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

import numpy as np
import matplotlib.pyplot as plt

# Parameters for the model
r_star = 2.0  # Neutral interest rate (% per year)
pi_star = 2.0  # Target inflation rate (% per year)
phi_pi = 1.5   # Reaction to inflation deviation
phi_Y = 0.5    # Reaction to output gap

# Initial conditions
T = 20  # Number of periods
time = np.arange(1, T + 1)

# Arrays to store values
inflation = np.zeros(T)
output_gap = np.zeros(T)
interest_rate = np.zeros(T)

# Initial values
inflation[0] = 3.0  # Initial inflation rate
output_gap[0] = -1.0  # Initial output gap

# Simulation loop
for t in range(1, T):
    # Determine interest rate using the Taylor Rule
    interest_rate[t] = (
        r_star + phi_pi * (inflation[t-1] - pi_star) + phi_Y * output_gap[t-1]
    )
    
    # Update output gap and inflation based on interest rate
    output_gap[t] = output_gap[t-1] - 0.5 * (interest_rate[t] - r_star)
    inflation[t] = inflation[t-1] - 0.3 * output_gap[t]

# Plotting the results
plt.figure(figsize=(12, 8))

# Plot inflation
plt.subplot(3, 1, 1)
plt.plot(time, inflation, label='Inflation Rate', color='red')
plt.axhline(pi_star, linestyle='--', color='black', label='Target Inflation Rate')
plt.title('Monetary Policy Simulation')
plt.ylabel('Inflation Rate (%)')
plt.legend()

# Plot output gap
plt.subplot(3, 1, 2)
plt.plot(time, output_gap, label='Output Gap', color='blue')
plt.axhline(0, linestyle='--', color='black', label='Zero Output Gap')
plt.ylabel('Output Gap (%)')
plt.legend()

# Plot interest rate
plt.subplot(3, 1, 3)
plt.plot(time, interest_rate, label='Interest Rate', color='green')
plt.axhline(r_star, linestyle='--', color='black', label='Neutral Interest Rate')
plt.xlabel('Time (Periods)')
plt.ylabel('Interest Rate (%)')
plt.legend()

plt.tight_layout()
plt.show()

Explanation of the Code

1. Model Setup:
   • The neutral interest rate $(r^*)$, target inflation rate $(\pi^*)$, and policy reaction coefficients $(\phi_\pi$, $\phi_Y)$ are defined.
2. Simulation Loop:
   • The interest rate is updated using the Taylor Rule.
   • The output gap and inflation rate are updated based on the interest rate, reflecting the dynamics of monetary policy.
3. Visualization:
   • The results for inflation, output gap, and interest rate over time are plotted.

Expected Results

1. Inflation Rate:
   Gradually converges toward the target rate $(\pi^* = 2％)$.
2. Output Gap:
   Oscillates and eventually stabilizes around zero.
3. Interest Rate:
   Adjusts dynamically based on the Taylor Rule to stabilize inflation and output.

The graphs will visually illustrate these dynamics, helping to understand how monetary policy impacts an economy over time.

Scenario

We aim to analyze and visualize the GDP growth of a country over a period of years.
Using a dataset with GDP values, we calculate:

1. Year-over-year growth rates.
2. An average annual growth rate over the period.
3. Visualization of the GDP trend and growth rates.

Python Implementation

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simulated GDP data for analysis
data = {
    "Year": np.arange(2000, 2021),  # Years from 2000 to 2020
    "GDP": [
        1.0, 1.1, 1.2, 1.35, 1.5, 1.65, 1.8, 2.0, 2.3, 2.5, 
        2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.7, 5.0
    ]  # Example GDP values in trillions
}

# Create a DataFrame
df = pd.DataFrame(data)

# Calculate year-over-year GDP growth rates
df["Growth Rate (%)"] = df["GDP"].pct_change() * 100

# Calculate the average annual growth rate
average_growth = df["Growth Rate (%)"].mean()

# Results summary
print("GDP Data with Growth Rates:")
print(df)

print(f"\nAverage Annual Growth Rate: {average_growth:.2f}%")

# Plotting GDP and growth rates
fig, ax1 = plt.subplots(figsize=(12, 6))

# Plot GDP (line chart)
ax1.plot(df["Year"], df["GDP"], color="blue", marker="o", label="GDP (Trillions)")
ax1.set_xlabel("Year")
ax1.set_ylabel("GDP (Trillions)", color="blue")
ax1.tick_params(axis="y", labelcolor="blue")
ax1.legend(loc="upper left")
ax1.grid(True, which="both", linestyle="--", linewidth=0.5)

# Plot Growth Rate (bar chart)
ax2 = ax1.twinx()
ax2.bar(df["Year"], df["Growth Rate (%)"], color="orange", alpha=0.6, label="Growth Rate (%)")
ax2.set_ylabel("Growth Rate (%)", color="orange")
ax2.tick_params(axis="y", labelcolor="orange")
ax2.legend(loc="upper right")

# Add title and show the plot
plt.title("GDP Growth Analysis (2000-2020)")
plt.tight_layout()
plt.show()

Code Explanation

1. Data Setup:

   • Simulated GDP data is provided for years $2000$–$2020$.
   • A DataFrame is created to store the year and GDP values.

2. Calculations:

   • The year-over-year growth rate is calculated using the percentage change formula:
     
   • The average annual growth rate is calculated as the mean of all growth rates.

3. Visualization:

   • GDP is plotted as a line chart to show the trend over time.
   • Growth rates are plotted as a bar chart to emphasize yearly changes.

Results and Graph:

GDP Data with Growth Rates:
    Year   GDP  Growth Rate (%)
0   2000  1.00              NaN
1   2001  1.10        10.000000
2   2002  1.20         9.090909
3   2003  1.35        12.500000
4   2004  1.50        11.111111
5   2005  1.65        10.000000
6   2006  1.80         9.090909
7   2007  2.00        11.111111
8   2008  2.30        15.000000
9   2009  2.50         8.695652
10  2010  2.80        12.000000
11  2011  3.00         7.142857
12  2012  3.20         6.666667
13  2013  3.40         6.250000
14  2014  3.60         5.882353
15  2015  3.80         5.555556
16  2016  4.00         5.263158
17  2017  4.20         5.000000
18  2018  4.40         4.761905
19  2019  4.70         6.818182
20  2020  5.00         6.382979

Average Annual Growth Rate: 8.42%

1. Data Analysis:

   • The printed DataFrame includes the GDP values and their corresponding growth rates.
   • The average annual growth rate is displayed for overall assessment.

2. Visualization:

   • Line Chart: Shows the upward trend of GDP over the years.
   • Bar Chart: Highlights variations in growth rates year-over-year.

This graph effectively communicates both the long-term growth trend and the yearly fluctuations in GDP growth rates, providing a clear picture of economic performance.

Scenario

We want to analyze the equilibrium price and quantity in a market using the following linear supply and demand functions:

1. Demand Function:
   $$
   Q_d = a - b \cdot P
   $$

   • $( Q_d )$: Quantity demanded
   • $( a )$: Demand intercept
   • $( b )$: Slope of demand curve
   • $( P )$: Price

2. Supply Function:
   $$
   Q_s = c + d \cdot P
   $$

   • $( Q_s )$: Quantity supplied
   • $( c )$: Supply intercept
   • $( d )$: Slope of supply curve

The equilibrium occurs where $( Q_d = Q_s )$.

Python Implementation

Below is the $Python$ code to find the equilibrium price and quantity, plot the supply and demand curves, and illustrate the equilibrium graphically.

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

import numpy as np
import matplotlib.pyplot as plt

# Parameters for demand and supply
a = 100  # Demand intercept
b = 2    # Demand slope
c = 20   # Supply intercept
d = 1.5  # Supply slope

# Price range
prices = np.linspace(0, 60, 100)

# Demand and supply functions
Q_d = a - b * prices  # Quantity demanded
Q_s = c + d * prices  # Quantity supplied

# Equilibrium: Solve for P where Q_d = Q_s
equilibrium_price = (a - c) / (b + d)
equilibrium_quantity = a - b * equilibrium_price

# Display equilibrium results
print(f"Equilibrium Price (P): {equilibrium_price:.2f}")
print(f"Equilibrium Quantity (Q): {equilibrium_quantity:.2f}")

# Plotting supply and demand curves
plt.figure(figsize=(10, 6))
plt.plot(prices, Q_d, label="Demand Curve (Qd)", color="blue")
plt.plot(prices, Q_s, label="Supply Curve (Qs)", color="orange")
plt.axvline(x=equilibrium_price, color="red", linestyle="--", label=f"Equilibrium Price: {equilibrium_price:.2f}")
plt.axhline(y=equilibrium_quantity, color="green", linestyle="--", label=f"Equilibrium Quantity: {equilibrium_quantity:.2f}")
plt.scatter([equilibrium_price], [equilibrium_quantity], color="purple", zorder=5, label="Equilibrium Point")

# Labels and legend
plt.title("Supply and Demand Analysis")
plt.xlabel("Price (P)")
plt.ylabel("Quantity (Q)")
plt.grid(True)
plt.legend()
plt.show()

Explanation of the Code

1. Functions:

   • The demand function $( Q_d = a - b \cdot P )$ decreases with price.
   • The supply function $( Q_s = c + d \cdot P )$ increases with price.

2. Equilibrium Calculation:

   • Solving $( Q_d = Q_s )$ algebraically yields $( P = (a - c) / (b + d) )$ and $( Q = a - b \cdot P )$.

3. Visualization:

   • The demand and supply curves are plotted across a price range.
   • The equilibrium price and quantity are highlighted on the graph with lines and a marker.

Results and Graphical Representation

Equilibrium Price (P): 22.86
Equilibrium Quantity (Q): 54.29

1. Equilibrium Values:

   • The equilibrium price and quantity are calculated and displayed.
   • For example, if $( a = 100, b = 2, c = 20, d = 1.5 )$:
     • Equilibrium price $( P \approx 26.67 )$
     • Equilibrium quantity $( Q \approx 46.67 )$

2. Graph:

   • The blue line represents the demand curve, while the orange line represents the supply curve.
   • The intersection point (marked in purple) is the equilibrium.

This example illustrates the core principle of supply and demand analysis, providing insights into how prices and quantities adjust in a competitive market.

Scenario

A firm produces a good using two inputs: labor (L) and capital (K).

The firm seeks to minimize its production costs while meeting a target output $( Q )$.

The cost function is:

$$
C = w \cdot L + r \cdot K
$$

• $( C )$: Total cost
• $( w )$: Wage rate (cost of labor)
• $( r )$: Rental rate of capital (cost of capital)
• $( L )$: Amount of labor used
• $( K )$: Amount of capital used

The production function is:
$$
Q = L^{0.5} \cdot K^{0.5}
$$

The firm must meet a specific output level $( Q = 10 )$.

The goal is to find the optimal values of $( L )$ and $( K )$ to minimize costs.

Python Implementation

Here is the $Python$ code to solve the cost minimization problem using SciPy.

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

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

# Parameters
w = 10  # Wage rate
r = 20  # Rental rate of capital
Q_target = 10  # Target output

# Production function
def production_function(L, K):
    return (L**0.5) * (K**0.5)

# Cost function
def cost_function(x):
    L, K = x
    return w * L + r * K

# Constraint: Output must equal Q_target
def output_constraint(x):
    L, K = x
    return production_function(L, K) - Q_target

# Initial guess
x0 = [1, 1]

# Bounds for L and K (cannot be negative)
bounds = [(0, None), (0, None)]

# Constraints
constraints = {'type': 'eq', 'fun': output_constraint}

# Optimization
result = minimize(cost_function, x0, bounds=bounds, constraints=constraints)

# Optimal values
L_opt, K_opt = result.x
min_cost = result.fun

# Display results
print(f"Optimal Labor (L): {L_opt:.2f}")
print(f"Optimal Capital (K): {K_opt:.2f}")
print(f"Minimum Cost: {min_cost:.2f}")

# Visualization
L_values = np.linspace(0.1, 15, 100)
K_values = Q_target**2 / L_values  # Rearrange Q = L^0.5 * K^0.5
cost_values = w * L_values + r * K_values

plt.figure(figsize=(10, 6))
plt.plot(L_values, cost_values, label="Cost Curve", color="blue", lw=2)
plt.scatter([L_opt], [min_cost], color="red", label="Optimal Point", zorder=5)
plt.title("Cost Minimization with Production Constraint")
plt.xlabel("Labor (L)")
plt.ylabel("Cost (C)")
plt.grid(True)
plt.legend()
plt.show()

Explanation of the Code

1. Cost Function:

   • The function $( C = w \cdot L + r \cdot K )$ calculates total cost for given $( L )$ and $( K )$.

2. Constraint:

   • The output constraint ensures the production level $( Q )$ equals the target $( Q_{target} )$.

3. Optimization:

   • scipy.optimize.minimize is used to minimize the cost function while satisfying the production constraint and non-negativity bounds for $( L )$ and $( K )$.

4. Visualization:

   • The cost curve is plotted for varying labor values $( L )$, with the corresponding $( K )$ values derived from the production constraint.
   • The optimal point is highlighted.

Results and Graphical Representation

Optimal Labor (L): 14.14
Optimal Capital (K): 7.07
Minimum Cost: 282.84

1. Optimal Values:

   • The solution provides the values of $( L )$ and $( K )$ that minimize costs while producing the target output.

2. Graph:

   • The blue curve represents the total cost for different allocations of labor and capital.
   • The red point marks the minimum cost, highlighting the optimal combination of $( L )$ and $( K )$.

This example demonstrates how firms can allocate resources efficiently to minimize production costs while meeting output targets.

Scenario

Consider a scenario where two individuals (Person $A$ and Person $B$) are dividing two resources: Resource X and Resource Y.

Each person has their utility function for the resources.

The goal is to identify the allocations that are Pareto efficient, meaning that no reallocation can make one person better off without making the other worse off.

Utility functions:

• Person $A$: $( U_A = x_A^{0.5} + y_A^{0.5} )$
• Person $B$: $( U_B = x_B^{0.5} + y_B^{0.5} )$

The total resources available are:

• $( X_{total} = 10 )$
• $( Y_{total} = 10 )$

We want to identify Pareto-efficient allocations of $ (x_A, y_A, x_B, y_B) $.

Python Implementation

Here is the $Python$ code to compute and visualize Pareto-efficient allocations:

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

import numpy as np
import matplotlib.pyplot as plt

# Total resources
X_total = 10
Y_total = 10

# Grid for allocations
x_A_values = np.linspace(0, X_total, 100)
y_A_values = np.linspace(0, Y_total, 100)

# Utility functions
def utility_A(x_A, y_A):
    return np.sqrt(x_A) + np.sqrt(y_A)

def utility_B(x_B, y_B):
    return np.sqrt(x_B) + np.sqrt(y_B)

# Compute Pareto-efficient points
pareto_points = []
for x_A in x_A_values:
    for y_A in y_A_values:
        x_B = X_total - x_A
        y_B = Y_total - y_A
        if x_B >= 0 and y_B >= 0:  # Valid allocations
            u_A = utility_A(x_A, y_A)
            u_B = utility_B(x_B, y_B)
            pareto_points.append((u_A, u_B))

# Convert Pareto points to arrays for plotting
pareto_points = np.array(pareto_points)
u_A_values = pareto_points[:, 0]
u_B_values = pareto_points[:, 1]

# Plot results
plt.figure(figsize=(8, 6))
plt.scatter(u_A_values, u_B_values, s=10, c="blue", label="Pareto-Efficient Points", alpha=0.6)
plt.title("Pareto-Efficient Allocations")
plt.xlabel("Utility of Person A (U_A)")
plt.ylabel("Utility of Person B (U_B)")
plt.grid(True)
plt.legend()
plt.show()

Explanation of the Code

1. Utility Functions:
   We defined the utility functions for both individuals based on their allocation of $( X )$ and $( Y )$.
2. Grid Search:
   A grid of possible allocations for $( x_A )$ and $( y_A )$ is generated.
   The remaining allocations for Person $B$ are calculated as $( x_B = X_{total} - x_A )$ and $( y_B = Y_{total} - y_A )$.
3. Pareto-Efficient Points:
   For each valid allocation, we compute the utilities of both individuals and store them as Pareto-efficient points.
4. Visualization:
   • The scatter plot shows the Pareto frontier, which represents the set of allocations where improving one individual’s utility would decrease the other’s.

Results and Graphical Representation

1. Pareto Frontier:
   • The blue points represent the possible Pareto-efficient allocations.
   • Any point on this curve is Pareto efficient because moving away from it would reduce the utility of at least one individual.
2. Utility Trade-off:
   • The graph clearly shows the trade-off between the utilities of Person $A$ and Person $B$.
   • Allocations can favor Person $A$ (higher $( U_A )$) or Person $B$ (higher $( U_B )$) depending on how the resources are distributed.

This approach demonstrates how Pareto efficiency can be visualized and used to evaluate resource allocation in economics or game theory.

Scenario

A company is hosting a sealed-bid auction where bidders independently submit bids without knowing the others’ bids.
The item is awarded to the highest bidder, and they pay their bid amount (a first-price auction).

We will simulate the auction process for $10$ bidders, each with a valuation drawn randomly from a uniform distribution between ＄$50$ and ＄$100$.
Bidders place bids by shading their valuations by a fraction.
The optimal shading strategy depends on the auction type.

Python Implementation

Here is the $Python$ code to simulate and analyze this auction:

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

import numpy as np
import matplotlib.pyplot as plt

# Auction setup
np.random.seed(42)  # For reproducibility
num_bidders = 10
valuations = np.random.uniform(50, 100, num_bidders)  # Valuations of bidders

# Bidding strategy: Shading factor
shading_factor = 0.8  # Bidders bid 80% of their valuation
bids = shading_factor * valuations

# Determine the winner
winning_bid = max(bids)
winner_index = np.argmax(bids)
winner_valuation = valuations[winner_index]

# Display results
print(f"Valuations: {valuations}")
print(f"Bids: {bids}")
print(f"Winning Bid: {winning_bid}")
print(f"Winner's Valuation: {winner_valuation}")
print(f"Winner Index: {winner_index}")

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

# Plot valuations and bids
x = np.arange(num_bidders)
plt.bar(x - 0.2, valuations, width=0.4, label="Valuations", color="blue", alpha=0.7)
plt.bar(x + 0.2, bids, width=0.4, label="Bids", color="orange", alpha=0.7)

# Highlight the winner
plt.scatter(winner_index, winning_bid, color="red", label="Winning Bid", s=100, zorder=5)

# Labels and legend
plt.xlabel("Bidders")
plt.ylabel("Amount")
plt.title("Valuations and Bids in a First-Price Auction")
plt.xticks(x, [f"Bidder {i}" for i in x])
plt.legend()
plt.grid(True)
plt.show()

Explanation of the Code

1. Valuations:
   Each bidder’s valuation is randomly sampled from a uniform distribution.
   This represents their willingness to pay.
2. Bids:
   A simple shading strategy is applied where each bidder bids $80$% of their valuation.
3. Winner Determination:
   The highest bid determines the winner, and their valuation and index are recorded.
4. Visualization:
   The bar chart compares valuations and bids, highlighting the winner.

Explanation of the Results

1. Valuations vs. Bids

• The blue bars represent the valuations of the $10$ bidders.
  These are the true amounts each bidder is willing to pay.
• The orange bars represent the bids submitted by the bidders.
  The bids are shaded down to $80$% of their valuations, as per the strategy.

2. Winning Bid

• The red dot marks the winning bid ($78.03$), submitted by Bidder $1$.
• The winner’s valuation ($97.54$) is the highest among all bidders, which aligns with the competitive nature of a first-price auction.

3. Bidder Strategy

• Since bidders bid only $80$% of their true valuations, no one submits their full willingness to pay.
  This is typical in a first-price auction because bidders want to maximize their surplus (difference between valuation and payment).

4. Observations

• Bidder $1$ wins the auction because their bid of $78.03$ is the highest among all bidders.
• Bidders with higher valuations tend to submit higher bids, as seen in the general trend where valuations (blue) and bids (orange) align proportionally.

Conclusion

This example highlights the effect of strategic bidding in a first-price auction:

• Bidders shade their bids below their valuations to achieve a positive surplus.
• The winner is determined based on the highest bid, but they still pay less than their full valuation.
  This behavior illustrates the trade-off between winning and maximizing surplus in competitive auctions.

Scenario

A company wants to understand how changes in the price of a product affect its demand.

The demand function is estimated from past data as:

$$
Q_d = 500 - 20P
$$

• $( Q_d )$ is the quantity demanded.
• $( P )$ is the price.

We will calculate the price elasticity of demand (PED) for different price levels using the formula:

$$
\text{PED} = \left( \frac{\Delta Q_d}{Q_d} \right) / \left( \frac{\Delta P}{P} \right)
$$

Python Implementation

Here is the $Python$ code to compute and visualize the elasticity of demand for a range of prices:

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

import numpy as np
import matplotlib.pyplot as plt

# Demand function
def demand_function(price):
    return 500 - 20 * price

# Compute price elasticity of demand
def price_elasticity(price, delta_price=0.01):
    # Initial quantity demanded
    initial_quantity = demand_function(price)
    # Quantity demanded after a small price change
    new_quantity = demand_function(price + delta_price)
    # Change in quantity demanded
    delta_quantity = new_quantity - initial_quantity
    # Elasticity calculation
    elasticity = (delta_quantity / initial_quantity) / (delta_price / price)
    return elasticity

# Generate data for analysis
prices = np.linspace(1, 20, 100)
quantities = demand_function(prices)
elasticities = [price_elasticity(p) for p in prices]

# Plot demand curve
plt.figure(figsize=(14, 6))
plt.subplot(1, 2, 1)
plt.plot(prices, quantities, label='Demand Curve', color='blue')
plt.title('Demand Curve')
plt.xlabel('Price')
plt.ylabel('Quantity Demanded')
plt.grid(True)
plt.legend()

# Plot price elasticity of demand
plt.subplot(1, 2, 2)
plt.plot(prices, elasticities, label='Elasticity', color='red')
plt.axhline(0, color='black', linewidth=0.8, linestyle='--')
plt.axhline(-1, color='green', linestyle='--', label='Unit Elasticity (PED=-1)')
plt.title('Price Elasticity of Demand')
plt.xlabel('Price')
plt.ylabel('Elasticity')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

Explanation of the Code

1. Demand Function:
   The function demand_function(price) represents the relationship between price and quantity demanded.
2. Elasticity Calculation:
   The price_elasticity function uses a small change in price $(\Delta P)$ to calculate the corresponding change in quantity $(\Delta Q_d)$, then computes elasticity.
3. Data Generation:
   A range of prices is simulated to compute demand and elasticity values.
4. Visualization:
   • The demand curve shows the relationship between price and quantity demanded.
   • The elasticity plot demonstrates how elastic or inelastic demand is at different price points.
     Unit elasticity is marked for reference.

Graphical Results

1. Demand Curve:
   • Downward-sloping curve indicates that higher prices reduce quantity demanded.
2. Elasticity Curve:
   • Shows the PED at different price points.
   • When PED > -1: Demand is inelastic (less sensitive to price changes).
   • When PED < -1: Demand is elastic (more sensitive to price changes).
   • PED = -1: Unit elastic, where the percentage change in quantity equals the percentage change in price.

This analysis provides insights into pricing strategies.
For example:

• Inelastic demand suggests raising prices to increase revenue.
• Elastic demand suggests lowering prices might increase revenue by boosting demand.