Problem Description

Simulate the trajectory of a projectile launched at an initial speed $( v_0 )$ and angle $( \theta )$ from the horizontal.

Assume:

1. The motion occurs in a vacuum (no air resistance).
2. The acceleration due to gravity is $( g = 9.8 , \text{m/s}^2 )$.

The equations of motion are:
$$
x(t) = v_0 \cos(\theta) t
$$
$$
y(t) = v_0 \sin(\theta) t - \frac{1}{2} g t^2
$$

• $ x(t) $: Horizontal position.
• $ y(t) $: Vertical position.
• $ t $: Time.

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

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.8  # Acceleration due to gravity (m/s^2)
v0 = 50  # Initial speed (m/s)
theta_deg = 45  # Launch angle in degrees
theta_rad = np.radians(theta_deg)  # Convert angle to radians

# Time of flight
t_flight = 2 * v0 * np.sin(theta_rad) / g

# Time array
t = np.linspace(0, t_flight, 500)

# Equations of motion
x = v0 * np.cos(theta_rad) * t
y = v0 * np.sin(theta_rad) * t - 0.5 * g * t**2

# Plot the trajectory
plt.figure(figsize=(10, 6))
plt.plot(x, y, label=f"Launch Angle: {theta_deg}°", color="blue")
plt.title("Projectile Motion Trajectory")
plt.xlabel("Horizontal Distance (m)")
plt.ylabel("Vertical Distance (m)")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants:
   Define gravity $( g )$, initial speed $( v_0 )$, and launch angle $( \theta )$.
2. Time of Flight:
   Calculate the total time the projectile remains in the air using $( t_{\text{flight}} = \frac{2 v_0 \sin(\theta)}{g} )$.
3. Equations of Motion:
   Compute horizontal $( x )$ and vertical $( y )$ positions for a range of time values $( t )$.
4. Plotting:
   Plot the $( x )$ - $( y )$ trajectory to visualize the motion.

Results and Visualization

The graph shows the parabolic trajectory of the projectile:

1. The projectile starts at the origin $( x = 0, y = 0 )$.
2. It reaches a maximum height midway through its flight.
3. The projectile lands back at ground level $( y = 0 )$ after covering a horizontal distance determined by the initial velocity and angle.

This visualization helps understand the fundamental aspects of projectile motion, including its symmetric nature and dependency on the launch angle and speed.

Problem Description

According to Einstein’s theory of special relativity, time dilation occurs for objects moving at significant fractions of the speed of light.

The formula for time dilation is:

$$
t’ = \frac{t}{\sqrt{1 - \frac{v^2}{c^2}}}
$$

• $( t’ )$: Time experienced by the moving observer (proper time).
• $( t )$: Time experienced by the stationary observer.
• $( v )$: Velocity of the moving object.
• $( c )$: Speed of light $(( \approx 3 \times 10^8 , \text{m/s} ))$.

Goal

Simulate how time dilation depends on the velocity $( v )$ of the moving object as a fraction of the speed of light $( c )$.

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

import numpy as np
import matplotlib.pyplot as plt

# Constants
speed_of_light = 3e8  # Speed of light (m/s)
stationary_time = 1.0  # Time experienced by the stationary observer (s)

# Velocity fractions of the speed of light
velocity_fractions = np.linspace(0, 0.99, 500)  # From 0 to 99% of c
velocities = velocity_fractions * speed_of_light

# Calculate time dilation
time_dilation = stationary_time / np.sqrt(1 - (velocities**2 / speed_of_light**2))

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(velocity_fractions, time_dilation, color='red', label="Time Dilation")
plt.title("Time Dilation vs. Velocity (Special Relativity)")
plt.xlabel("Velocity as a Fraction of the Speed of Light (v/c)")
plt.ylabel("Dilated Time (t')")
plt.grid(True)
plt.legend()
plt.show()

Code Explanation

1. Constants:
   The speed of light and stationary observer’s time are defined.
2. Velocity Fractions:
   The object’s velocity is expressed as fractions of $( c )$, ranging from $0$ to $0.99$ ($99$% of $( c )$).
3. Time Dilation Formula:
   The time dilation is computed using the relativistic formula.
4. Plotting:
   The relationship between velocity (as a fraction of $( c )$) and time dilation is plotted.

Results and Visualization

The graph illustrates the effect of velocity on time dilation:

1. At low velocities $(( v/c \approx 0 ))$, time dilation is minimal $(( t’ \approx t ))$.
2. As velocity approaches the speed of light $(( v/c \to 1 ))$, time dilation increases drastically $(( t’ \to \infty ))$.
3. This demonstrates how relativistic effects become significant only at very high velocities, a key feature of special relativity.

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.