Example Problem

We will simulate the orbit of a planet around a star using Kepler’s laws of planetary motion and Newtonian mechanics.

The goal is to calculate the orbit of the planet and visualize it.

Problem Description

1. Kepler’s Laws:

   • Planets move in elliptical orbits with the star at one focus.
   • The line connecting a planet to the star sweeps out equal areas in equal times.

2. Newtonian Mechanics:

   • The gravitational force between the star and the planet dictates the motion:
     $$
     F = \frac{G M m}{r^2}
     $$
   • The planet’s position is updated using numerical integration of this force.

Python Code Implementation

Here’s the $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

import numpy as np
import matplotlib.pyplot as plt

# Constants
G = 6.67430e-11  # Gravitational constant, m^3 kg^-1 s^-2
M = 1.989e30  # Mass of the star (e.g., the Sun), kg
m = 5.972e24  # Mass of the planet (e.g., the Earth), kg
AU = 1.496e11  # Astronomical Unit, m (average Earth-Sun distance)

# Initial conditions
r0 = AU  # Initial distance, m
v0 = 30000  # Initial tangential velocity, m/s
dt = 60 * 60  # Time step, seconds
time_steps = 10000  # Number of time steps

# Arrays to store position and velocity
pos = np.zeros((time_steps, 2))
vel = np.zeros((time_steps, 2))
acc = np.zeros((time_steps, 2))

# Initial position and velocity
pos[0] = [r0, 0]
vel[0] = [0, v0]

# Simulation loop
for t in range(1, time_steps):
    # Distance and gravitational force
    r = np.linalg.norm(pos[t - 1])
    F = -G * M * m / r**2
    acc[t - 1] = F / m * pos[t - 1] / r

    # Update velocity and position
    vel[t] = vel[t - 1] + acc[t - 1] * dt
    pos[t] = pos[t - 1] + vel[t] * dt

# Convert to AU for visualization
pos_au = pos / AU

# Plotting the orbit
plt.figure(figsize=(8, 8))
plt.plot(pos_au[:, 0], pos_au[:, 1], label="Planet Orbit")
plt.scatter(0, 0, color="orange", label="Star (e.g., Sun)", s=100)
plt.title("Planetary Orbit Simulation")
plt.xlabel("X Position (AU)")
plt.ylabel("Y Position (AU)")
plt.axis("equal")
plt.grid()
plt.legend()
plt.show()

Code Explanation

1. Initial Conditions:

   • The planet starts at a distance $( r_0 ) (1 AU)$ with an initial velocity perpendicular to the radius.

2. Numerical Integration:

   • At each time step:
     • Calculate the gravitational acceleration $( a = F / m )$.
     • Update the velocity using $( v = v + a \cdot dt )$.
     • Update the position using $( r = r + v \cdot dt )$.

3. Scaling:

   • The positions are scaled to astronomical units $(AU)$ for clear visualization.

4. Plot:

   • The orbit is plotted, showing the elliptical path of the planet around the star.

Results Visualization

1. Orbit Shape:

   • The planet’s trajectory forms an ellipse with the star at one focus.

2. Key Observations:

   • The planet moves faster when it is closer to the star (perihelion).
   • The planet moves slower when it is farther from the star (aphelion), illustrating Kepler’s second law.

This simulation demonstrates how planetary orbits emerge from gravitational interactions.

Example Problem

We will solve an example of both elastic and inelastic collisions in one dimension.

The simulation will calculate the velocities after the collision for two masses and illustrate their behavior.

Problem Description

1. Elastic Collision:

   • Two objects, $(m_1)$ and $(m_2)$, collide with initial velocities $(v_{1i})$ and $(v_{2i})$.
   • Conservation of momentum and kinetic energy holds:
     $$
     m_1 v_{1i} + m_2 v_{2i} = m_1 v_{1f} + m_2 v_{2f}
     $$
     $$
     \frac{1}{2} m_1 v_{1i}^2 + \frac{1}{2} m_2 v_{2i}^2 = \frac{1}{2} m_1 v_{1f}^2 + \frac{1}{2} m_2 v_{2f}^2
     $$

2. Inelastic Collision:

   • Two objects stick together after the collision.
   • Only momentum is conserved:
     $$
     m_1 v_{1i} + m_2 v_{2i} = (m_1 + m_2) v_f
     $$

Python Code Implementation

Here’s the code to calculate and visualize the results:

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

import numpy as np
import matplotlib.pyplot as plt

# Define the masses and initial velocities
m1 = 2.0  # mass of object 1
m2 = 1.0  # mass of object 2
v1i = 4.0  # initial velocity of object 1
v2i = -2.0  # initial velocity of object 2

# Elastic collision calculations
v1f_elastic = ((m1 - m2) * v1i + 2 * m2 * v2i) / (m1 + m2)
v2f_elastic = ((m2 - m1) * v2i + 2 * m1 * v1i) / (m1 + m2)

# Inelastic collision calculations
vf_inelastic = (m1 * v1i + m2 * v2i) / (m1 + m2)

# Time and positions for visualization
time = np.linspace(0, 2, 100)
pos1_elastic = v1i * time
pos2_elastic = v2i * time
pos1_elastic_after = v1f_elastic * time
pos2_elastic_after = v2f_elastic * time
pos_inelastic_after = vf_inelastic * time

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

# Elastic collision
plt.subplot(1, 2, 1)
plt.plot(time, pos1_elastic, label="Object 1 (Before)", linestyle="--")
plt.plot(time, pos2_elastic, label="Object 2 (Before)", linestyle="--")
plt.plot(time, pos1_elastic_after, label="Object 1 (After)")
plt.plot(time, pos2_elastic_after, label="Object 2 (After)")
plt.title("Elastic Collision")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.grid()

# Inelastic collision
plt.subplot(1, 2, 2)
plt.plot(time, pos1_elastic, label="Object 1 (Before)", linestyle="--")
plt.plot(time, pos2_elastic, label="Object 2 (Before)", linestyle="--")
plt.plot(time, pos_inelastic_after, label="Combined (After)")
plt.title("Inelastic Collision")
plt.xlabel("Time (s)")
plt.ylabel("Position (m)")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

Code Explanation

1. Elastic Collision:

   • The final velocities $(v_{1f})$ and $(v_{2f})$ are derived using the conservation equations.
   • These velocities are used to calculate the positions over time for visualization.

2. Inelastic Collision:

   • A single final velocity $(v_f)$ is computed since the objects stick together.
   • This velocity is used to show the combined motion.

3. Visualization:

   • Two subplots illustrate the motion before and after both types of collisions.

Results Visualization

1. Elastic Collision:

   • The two objects bounce off each other, maintaining kinetic energy.
   • Their positions diverge after the collision.

2. Inelastic Collision:

   • The two objects stick together and move with the same final velocity.
   • Their positions are represented by a single line post-collision.

This simulation provides a clear comparison of elastic and inelastic collisions.

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.