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

# Function to generate a Möbius strip
def generate_mobius_strip(u_samples=100, v_samples=50):
    """
    Generate the x, y, z coordinates of a Möbius strip.

    Parameters:
        u_samples: Number of samples along the strip (angle direction)
        v_samples: Number of samples across the strip (width direction)

    Returns:
        x, y, z: Coordinates of the Möbius strip
    """
    # Create parameter ranges for u (angle) and v (width)
    u = np.linspace(0, 2 * np.pi, u_samples)
    v = np.linspace(-0.5, 0.5, v_samples)
    u, v = np.meshgrid(u, v)

    # Parametric equations for a Möbius strip
    x = (1 + v * np.cos(u / 2)) * np.cos(u)
    y = (1 + v * np.cos(u / 2)) * np.sin(u)
    z = v * np.sin(u / 2)

    return x, y, z

# Generate Möbius strip data
x, y, z = generate_mobius_strip()

# Plot the Möbius strip
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# Surface plot with color gradient
ax.plot_surface(x, y, z, cmap='viridis', edgecolor='k', alpha=0.9)

# Set plot aesthetics
ax.set_title("Möbius Strip", fontsize=16)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Function to calculate the circumcenter and radius
def find_circumcircle(A, B, C):
    # Coordinates of the points
    x1, y1 = A
    x2, y2 = B
    x3, y3 = C

    # Calculate the slopes of the perpendicular bisectors
    D = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
    if D == 0:
        raise ValueError("Points are collinear; circumcircle cannot be determined.")

    # Circumcenter coordinates
    Ux = ((x1**2 + y1**2) * (y2 - y3) + (x2**2 + y2**2) * (y3 - y1) + (x3**2 + y3**2) * (y1 - y2)) / D
    Uy = ((x1**2 + y1**2) * (x3 - x2) + (x2**2 + y2**2) * (x1 - x3) + (x3**2 + y3**2) * (x2 - x1)) / D

    # Radius of the circumcircle
    radius = np.sqrt((x1 - Ux)**2 + (y1 - Uy)**2)

    return (Ux, Uy), radius

# Function to plot the triangle and its circumcircle
def plot_triangle_and_circumcircle(A, B, C):
    # Calculate circumcenter and radius
    circumcenter, radius = find_circumcircle(A, B, C)

    # Create the circle
    circle = plt.Circle(circumcenter, radius, color='blue', fill=False, linestyle='--', label='Circumcircle')

    # Plot the triangle
    x_coords = [A[0], B[0], C[0], A[0]]
    y_coords = [A[1], B[1], C[1], A[1]]

    plt.figure(figsize=(8, 8))
    plt.plot(x_coords, y_coords, '-o', label='Triangle', color='red')
    plt.scatter(*circumcenter, color='green', label='Circumcenter')
    plt.gca().add_artist(circle)

    # Formatting
    plt.title("Triangle and its Circumcircle")
    plt.xlabel("X-axis")
    plt.ylabel("Y-axis")
    plt.axis('equal')
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.legend()
    plt.show()

# Example points
A = (1, 1)
B = (7, 1)
C = (4, 5)

# Plot the triangle and circumcircle
plot_triangle_and_circumcircle(A, B, C)

import time
import random
import matplotlib.pyplot as plt

# Step 1: Implement sorting algorithms

# Bubble Sort
def bubble_sort(arr):
    n = len(arr)
    for i in range(n):
        for j in range(0, n - i - 1):
            if arr[j] > arr[j + 1]:
                arr[j], arr[j + 1] = arr[j + 1], arr[j]

# Merge Sort
def merge_sort(arr):
    if len(arr) > 1:
        mid = len(arr) // 2
        L = arr[:mid]
        R = arr[mid:]

        merge_sort(L)
        merge_sort(R)

        i = j = k = 0

        while i < len(L) and j < len(R):
            if L[i] < R[j]:
                arr[k] = L[i]
                i += 1
            else:
                arr[k] = R[j]
                j += 1
            k += 1

        while i < len(L):
            arr[k] = L[i]
            i += 1
            k += 1

        while j < len(R):
            arr[k] = R[j]
            j += 1
            k += 1

# Step 2: Test and time the algorithms
def test_sorting_algorithms():
    sizes = [100, 500, 1000, 5000, 10000]  # List sizes to test
    bubble_times = []
    merge_times = []
    timsort_times = []

    for size in sizes:
        arr = [random.randint(0, 100000) for _ in range(size)]

        # Bubble Sort
        arr_copy = arr[:]
        start_time = time.time()
        bubble_sort(arr_copy)
        bubble_times.append(time.time() - start_time)

        # Merge Sort
        arr_copy = arr[:]
        start_time = time.time()
        merge_sort(arr_copy)
        merge_times.append(time.time() - start_time)

        # Timsort (Python's built-in sorted())
        arr_copy = arr[:]
        start_time = time.time()
        sorted(arr_copy)
        timsort_times.append(time.time() - start_time)

    return sizes, bubble_times, merge_times, timsort_times

# Step 3: Run the test
sizes, bubble_times, merge_times, timsort_times = test_sorting_algorithms()

# Step 4: Visualize the results
plt.figure(figsize=(10, 6))
plt.plot(sizes, bubble_times, label="Bubble Sort", marker='o', color='red')
plt.plot(sizes, merge_times, label="Merge Sort", marker='o', color='blue')
plt.plot(sizes, timsort_times, label="Timsort", marker='o', color='green')
plt.title("Sorting Algorithm Performance")
plt.xlabel("List Size")
plt.ylabel("Time Taken (seconds)")
plt.yscale("log")  # Use logarithmic scale for better comparison
plt.grid(True, linestyle="--", alpha=0.6)
plt.legend()
plt.show()

import matplotlib.pyplot as plt

# Function to compute the Collatz sequence
def collatz_sequence(n):
    steps = 0
    sequence = [n]
    while n != 1:
        if n % 2 == 0:
            n = n // 2
        else:
            n = 3 * n + 1
        sequence.append(n)
        steps += 1
    return sequence, steps

# Simulate Collatz for multiple starting numbers
start_numbers = range(1, 51)  # Starting numbers from 1 to 50
step_counts = []  # To store the number of steps for each starting number

# Compute the sequence for each starting number
for start in start_numbers:
    _, steps = collatz_sequence(start)
    step_counts.append(steps)

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

# Plot the number of steps for each starting number
plt.bar(start_numbers, step_counts, color="skyblue", edgecolor="black")
plt.title("Collatz Conjecture: Number of Steps to Reach 1")
plt.xlabel("Starting Number")
plt.ylabel("Number of Steps")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.show()

# Example sequence for a specific number
n_example = 27
sequence, _ = collatz_sequence(n_example)

# Plot the sequence for the example starting number
plt.figure(figsize=(12, 6))
plt.plot(range(len(sequence)), sequence, marker="o", color="red", label=f"Starting Number: {n_example}")
plt.title(f"Collatz Sequence for Starting Number {n_example}")
plt.xlabel("Step")
plt.ylabel("Value")
plt.grid(True, linestyle="--", alpha=0.7)
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Set random seed for reproducibility
np.random.seed(42)

def simulate_random_walk(n_steps, drift=0.1, volatility=1.0, paths=5):
    """
    Simulate random walk with drift
    Parameters:
    - n_steps: number of time steps
    - drift: trend component
    - volatility: size of random movements
    - paths: number of paths to simulate
    """
    steps = np.random.normal(drift, volatility, (paths, n_steps))
    paths = np.cumsum(steps, axis=1)
    return paths

def simulate_ornstein_uhlenbeck(n_steps, theta=1.0, mu=0.0, sigma=0.5, paths=5):
    """
    Simulate Ornstein-Uhlenbeck process
    Parameters:
    - n_steps: number of time steps
    - theta: mean reversion rate
    - mu: long-term mean
    - sigma: volatility
    - paths: number of paths to simulate
    """
    dt = 0.01
    t = np.linspace(0, n_steps*dt, n_steps)
    paths_ou = np.zeros((paths, n_steps))
    
    for i in range(paths):
        x = np.zeros(n_steps)
        x[0] = mu
        for j in range(1, n_steps):
            dx = theta * (mu - x[j-1]) * dt + sigma * np.sqrt(dt) * np.random.normal()
            x[j] = x[j-1] + dx
        paths_ou[i] = x
    
    return t, paths_ou

# Simulation parameters
n_steps = 1000
n_paths = 5

# Simulate processes
rw_paths = simulate_random_walk(n_steps, drift=0.1, volatility=0.5, paths=n_paths)
t_ou, ou_paths = simulate_ornstein_uhlenbeck(n_steps, theta=2.0, mu=0.0, sigma=1.0, paths=n_paths)

# Create visualizations
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

# Plot Random Walk paths
time = np.arange(n_steps)
for i in range(n_paths):
    ax1.plot(time, rw_paths[i], label=f'Path {i+1}')
ax1.set_title('Random Walk with Drift')
ax1.set_xlabel('Time Steps')
ax1.set_ylabel('Value')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Add theoretical drift line
drift_line = 0.1 * time
ax1.plot(time, drift_line, '--', color='black', label='Theoretical Drift', alpha=0.5)

# Plot Ornstein-Uhlenbeck paths
for i in range(n_paths):
    ax2.plot(t_ou, ou_paths[i], label=f'Path {i+1}')
ax2.axhline(y=0, color='black', linestyle='--', label='Long-term Mean', alpha=0.5)
ax2.set_title('Ornstein-Uhlenbeck Process')
ax2.set_xlabel('Time')
ax2.set_ylabel('Value')
ax2.grid(True, alpha=0.3)
ax2.legend()

plt.tight_layout()
plt.show()

# Calculate and print statistics
print("Random Walk Statistics:")
print(f"Final values: {rw_paths[:, -1]}")
print(f"Mean final value: {np.mean(rw_paths[:, -1]):.3f}")
print(f"Std of final values: {np.std(rw_paths[:, -1]):.3f}")

print("\nOrnstein-Uhlenbeck Statistics:")
print(f"Final values: {ou_paths[:, -1]}")
print(f"Mean final value: {np.mean(ou_paths[:, -1]):.3f}")
print(f"Std of final values: {np.std(ou_paths[:, -1]):.3f}")

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

# Parameters of the normal distribution
mu = 50  # Mean
sigma = 10  # Standard deviation
sample_size = 1000

# Generate a random sample
np.random.seed(42)  # For reproducibility
sample = np.random.normal(mu, sigma, sample_size)

# Compute sample statistics
sample_mean = np.mean(sample)
sample_std = np.std(sample)

print(f"Sample Mean: {sample_mean:.2f}")
print(f"Sample Standard Deviation: {sample_std:.2f}")

# Visualization
x = np.linspace(mu - 4 * sigma, mu + 4 * sigma, 1000)  # x-axis range
pdf = norm.pdf(x, mu, sigma)  # Theoretical PDF

plt.figure(figsize=(10, 6))

# Plot the histogram of the sample
plt.hist(sample, bins=30, density=True, alpha=0.6, color='skyblue', label="Sample Histogram")

# Overlay the theoretical PDF
plt.plot(x, pdf, color='red', linewidth=2, label="Theoretical PDF")

# Add labels and legend
plt.title("Normal Distribution: Histogram and Theoretical PDF")
plt.xlabel("Value")
plt.ylabel("Density")
plt.axvline(sample_mean, color='green', linestyle='--', label=f"Sample Mean: {sample_mean:.2f}")
plt.axvline(mu, color='purple', linestyle=':', label=f"Theoretical Mean: {mu}")
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# Define the function and its derivative
def f(x):
    return x**3 - 6*x**2 + 11*x - 6.1

def f_prime(x):
    return 3*x**2 - 12*x + 11

# Newton-Raphson Method
def newton_raphson(f, f_prime, x0, tol=1e-6, max_iter=100):
    x = x0
    iteration = 0
    errors = []

    for _ in range(max_iter):
        # Calculate next point
        x_new = x - f(x) / f_prime(x)
        error = abs(x_new - x)
        errors.append(error)

        # Check for convergence
        if error < tol:
            return x_new, iteration + 1, errors

        x = x_new
        iteration += 1

    raise ValueError("Newton-Raphson method did not converge.")

# Initial guess
x0 = 3.5

# Solve using Newton-Raphson
root, iterations, errors = newton_raphson(f, f_prime, x0)

# Display results
print(f"Root found: {root:.6f}")
print(f"Iterations: {iterations}")

# Plot results
x_vals = np.linspace(2, 4, 100)
y_vals = f(x_vals)

plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals, label="f(x)", color="blue")
plt.axhline(0, color="black", linewidth=0.8, linestyle="--")
plt.scatter(root, f(root), color="red", label=f"Root: x = {root:.6f}", zorder=5)
plt.title("Root Finding Using Newton-Raphson Method")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.grid(True)
plt.show()

# Error Convergence Plot
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(errors) + 1), errors, marker='o', color="green")
plt.title("Error Convergence in Newton-Raphson Method")
plt.xlabel("Iteration")
plt.ylabel("Error")
plt.yscale("log")  # Log scale for better visualization
plt.grid(True, which="both", linestyle="--", linewidth=0.5)
plt.show()

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

class HarmonicOscillator:
    def __init__(self, mass=1.0, spring_constant=10.0, damping=0.5):
        """
        Initialize harmonic oscillator parameters
        mass: mass of the object (kg)
        spring_constant: spring constant k (N/m)
        damping: damping coefficient c (Ns/m)
        """
        self.m = mass
        self.k = spring_constant
        self.c = damping
        
        # Derived parameters
        self.omega0 = np.sqrt(self.k / self.m)  # Natural frequency
        self.gamma = self.c / (2 * self.m)      # Damping ratio
        
    def equations_of_motion(self, state, t):
        """
        Define the system of differential equations
        state[0] = x (position)
        state[1] = v (velocity)
        """
        x, v = state
        dx_dt = v
        dv_dt = (-self.k * x - self.c * v) / self.m
        return [dx_dt, dv_dt]
    
    def simulate(self, t_span, initial_conditions):
        """
        Simulate the oscillator motion
        t_span: time points for simulation
        initial_conditions: [x0, v0] initial position and velocity
        """
        solution = odeint(self.equations_of_motion, initial_conditions, t_span)
        return solution
    
    def calculate_energy(self, x, v):
        """Calculate kinetic and potential energy"""
        KE = 0.5 * self.m * v**2
        PE = 0.5 * self.k * x**2
        return KE, PE

# Create simulation parameters
oscillator = HarmonicOscillator(mass=1.0, spring_constant=10.0, damping=0.5)
t_span = np.linspace(0, 10, 1000)
initial_conditions = [1.0, 0.0]  # x0 = 1m, v0 = 0 m/s

# Run simulation
solution = oscillator.simulate(t_span, initial_conditions)
x = solution[:, 0]  # Position
v = solution[:, 1]  # Velocity

# Calculate energies
KE, PE = oscillator.calculate_energy(x, v)
total_energy = KE + PE

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

# Position and Velocity vs Time
plt.subplot(2, 2, 1)
plt.plot(t_span, x, 'b-', label='Position')
plt.plot(t_span, v, 'r--', label='Velocity')
plt.xlabel('Time (s)')
plt.ylabel('Position (m) / Velocity (m/s)')
plt.title('Position and Velocity vs Time')
plt.grid(True)
plt.legend()

# Phase Space Plot
plt.subplot(2, 2, 2)
plt.plot(x, v, 'g-')
plt.xlabel('Position (m)')
plt.ylabel('Velocity (m/s)')
plt.title('Phase Space Plot')
plt.grid(True)

# Energy vs Time
plt.subplot(2, 2, 3)
plt.plot(t_span, KE, 'r-', label='Kinetic Energy')
plt.plot(t_span, PE, 'b-', label='Potential Energy')
plt.plot(t_span, total_energy, 'g-', label='Total Energy')
plt.xlabel('Time (s)')
plt.ylabel('Energy (J)')
plt.title('Energy vs Time')
plt.grid(True)
plt.legend()

# Power Spectrum
plt.subplot(2, 2, 4)
freq = np.fft.fftfreq(t_span.size, t_span[1] - t_span[0])
x_fft = np.fft.fft(x)
plt.plot(freq[freq > 0], np.abs(x_fft)[freq > 0])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('Power Spectrum')
plt.grid(True)

plt.tight_layout()
plt.show()

# Print system parameters
print(f"System Parameters:")
print(f"Mass: {oscillator.m} kg")
print(f"Spring Constant: {oscillator.k} N/m")
print(f"Damping Coefficient: {oscillator.c} Ns/m")
print(f"Natural Frequency: {oscillator.omega0:.2f} rad/s")
print(f"Damping Ratio: {oscillator.gamma:.2f}")

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

# Define global parameters
A = 1.0      # Total factor productivity
alpha = 0.3  # Output elasticity of capital
beta = 0.7   # Output elasticity of labor

def cobb_douglas(K, L):
    """
    Cobb-Douglas Production Function
    Y = A * K^α * L^β
    """
    return A * (K**alpha) * (L**beta)

def marginal_products(K, L):
    """Calculate marginal products of capital and labor"""
    MPK = A * alpha * (K**(alpha-1)) * (L**beta)
    MPL = A * beta * (K**alpha) * (L**(beta-1))
    return MPK, MPL

# Create data points
K = np.linspace(0.1, 10, 100)
L = np.linspace(0.1, 10, 100)
K_grid, L_grid = np.meshgrid(K, L)

# Calculate production
Y = cobb_douglas(K_grid, L_grid)

# Calculate marginal products for a specific point
K_point = 5
L_point = 5
MPK, MPL = marginal_products(K_point, L_point)

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

# 3D Production Surface
ax1 = fig.add_subplot(221, projection='3d')
surf = ax1.plot_surface(K_grid, L_grid, Y, cmap='viridis')
ax1.set_xlabel('Capital (K)')
ax1.set_ylabel('Labor (L)')
ax1.set_zlabel('Output (Y)')
ax1.set_title('Cobb-Douglas Production Surface')
fig.colorbar(surf)

# Production Contour Plot
ax2 = fig.add_subplot(222)
contour = ax2.contour(K_grid, L_grid, Y, levels=15)
ax2.clabel(contour, inline=True, fontsize=8)
ax2.set_xlabel('Capital (K)')
ax2.set_ylabel('Labor (L)')
ax2.set_title('Production Isoquants')

# Returns to Scale Analysis
K_scale = np.linspace(0.1, 10, 100)
Y_original = cobb_douglas(K_scale, K_scale)
Y_doubled = cobb_douglas(2*K_scale, 2*K_scale)

ax3 = fig.add_subplot(223)
ax3.plot(K_scale, Y_original, label='Original Scale')
ax3.plot(K_scale, Y_doubled/2, 'r--', label='Doubled Inputs/2')
ax3.set_xlabel('Input Scale')
ax3.set_ylabel('Output/Scale')
ax3.set_title('Returns to Scale Analysis')
ax3.legend()
ax3.grid(True)

# Marginal Products Analysis
MPK_curve = A * alpha * (K_scale**(alpha-1)) * (5**beta)
MPL_curve = A * beta * (K_scale**alpha) * (5**(beta-1))
ax4 = fig.add_subplot(224)
ax4.plot(K_scale, MPK_curve, label='MPK (L=5)')
ax4.plot(K_scale, MPL_curve, label='MPL (L=5)')
ax4.set_xlabel('Capital (K)')
ax4.set_ylabel('Marginal Products')
ax4.set_title('Marginal Products Analysis')
ax4.legend()
ax4.grid(True)

plt.tight_layout()
plt.show()

# Print economic analysis
print(f"\nEconomic Analysis at K={K_point}, L={L_point}:")
print(f"Output: {cobb_douglas(K_point, L_point):.2f}")
print(f"Marginal Product of Capital: {MPK:.2f}")
print(f"Marginal Product of Labor: {MPL:.2f}")
print(f"Returns to Scale: {alpha + beta:.2f}")  # Should be 1.0 for constant returns

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

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

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

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

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

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

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

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

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

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

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

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

plt.tight_layout()
plt.show()