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

# Constants
H0 = 70  # Hubble constant in km/s/Mpc
H0_SI = H0 * 3.24078e-20  # Convert to 1/s
Omega_m = 1.0  # Matter density parameter

# Friedmann equation as a function of time
def friedmann(t, a):
    return [H0_SI * np.sqrt(Omega_m / a[0]**3)]

# Time span: from the Big Bang (t=0) to 14 billion years
t_start = 0  # Initial time
t_end = 14e9 * 365.25 * 24 * 3600  # 14 billion years in seconds
a_initial = [1e-5]  # Small initial scale factor (close to 0)

# Solve the Friedmann equation
solution = solve_ivp(
    friedmann,
    [t_start, t_end],
    a_initial,
    method="RK45",
    dense_output=True
)

# Extract time (in billion years) and scale factor
time = solution.t / (365.25 * 24 * 3600 * 1e9)  # Convert to billion years
scale_factor = solution.y[0]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(time, scale_factor, label="Scale Factor $a(t)$", color="blue")
plt.axvline(0, color="black", linestyle="--", linewidth=0.8)
plt.axhline(0, color="black", linestyle="--", linewidth=0.8)

# Labels and legend
plt.title("Evolution of the Scale Factor $a(t)$", fontsize=14)
plt.xlabel("Time (Billion Years)", fontsize=12)
plt.ylabel("Scale Factor $a(t)$", fontsize=12)
plt.grid(alpha=0.5)
plt.legend(fontsize=12)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def shannon_entropy(probabilities):
    """
    Calculate the Shannon entropy of a discrete probability distribution.

    Parameters:
        probabilities: List or array of probabilities (must sum to 1).

    Returns:
        entropy: Shannon entropy in bits.
    """
    probabilities = np.array(probabilities)
    # Ensure probabilities are valid
    assert np.isclose(np.sum(probabilities), 1), "Probabilities must sum to 1."
    assert np.all(probabilities >= 0), "Probabilities cannot be negative."

    # Calculate entropy
    entropy = -np.sum(probabilities * np.log2(probabilities + 1e-10))  # Add small value to avoid log(0)
    return entropy

# Define example probability distributions
distributions = {
    "Uniform": [0.25, 0.25, 0.25, 0.25],
    "Biased": [0.7, 0.1, 0.1, 0.1],
    "Highly Skewed": [0.95, 0.05],
    "Equal Binary": [0.5, 0.5],
}

# Calculate entropy for each distribution
entropies = {name: shannon_entropy(dist) for name, dist in distributions.items()}

# Display results
for name, entropy in entropies.items():
    print(f"Entropy of {name} distribution: {entropy:.4f} bits")

# Visualization
fig, ax = plt.subplots(figsize=(10, 6))

# Plot each distribution
for name, dist in distributions.items():
    x = range(len(dist))
    ax.bar(x, dist, alpha=0.7, label=f"{name} (H={entropies[name]:.2f})")
    for i, prob in enumerate(dist):
        ax.text(i, prob + 0.02, f"{prob:.2f}", ha='center', fontsize=10)

# Add labels, legend, and title
ax.set_title("Discrete Probability Distributions and Their Entropy", fontsize=16)
ax.set_xlabel("Events", fontsize=12)
ax.set_ylabel("Probability", fontsize=12)
ax.set_ylim(0, 1.1)
ax.legend()
plt.show()

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}")