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()

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

def solve_2d_heat_equation(dx=0.1, dy=0.1, dt=0.0001, t_final=0.01):
    """
    Solve 2D heat equation using finite difference method
    ∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²)
    with initial conditions: u(x,y,0) = sin(πx)sin(πy)
    and boundary conditions: u = 0 at boundaries
    """
    # Thermal diffusivity
    alpha = 1.0
    
    # Create grid
    x = np.arange(0, 1+dx, dx)
    y = np.arange(0, 1+dy, dy)
    t = np.arange(0, t_final+dt, dt)
    
    # Initialize temperature grid
    nx, ny = len(x), len(y)
    u = np.zeros((nx, ny))
    
    # Set initial conditions
    X, Y = np.meshgrid(x, y)
    u = np.sin(np.pi*X) * np.sin(np.pi*Y)
    
    # Create arrays for new and old solutions
    u_new = np.zeros((nx, ny))
    
    # Stability condition check
    stability = alpha * dt * (1/(dx**2) + 1/(dy**2))
    print(f"Stability parameter: {stability:.3f} (should be < 0.5)")
    
    # Time stepping
    for n in range(len(t)-1):
        # Interior points
        for i in range(1, nx-1):
            for j in range(1, ny-1):
                u_new[i,j] = u[i,j] + alpha*dt*(
                    (u[i+1,j] - 2*u[i,j] + u[i-1,j])/dx**2 +
                    (u[i,j+1] - 2*u[i,j] + u[i,j-1])/dy**2
                )
        
        # Update solution
        u = u_new.copy()
        
        # Enforce boundary conditions
        u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 0
    
    return X, Y, u

# Solve the equation
X, Y, u_final = solve_2d_heat_equation()

# Create visualization
fig = plt.figure(figsize=(12, 5))

# Surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(X, Y, u_final, cmap=cm.viridis)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('Temperature')
ax1.set_title('2D Heat Distribution (3D View)')
fig.colorbar(surf)

# Contour plot
ax2 = fig.add_subplot(122)
contour = ax2.contourf(X, Y, u_final, levels=20, cmap=cm.viridis)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('2D Heat Distribution (Top View)')
fig.colorbar(contour)

plt.tight_layout()
plt.show()

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

class PipeFlowAnalyzer:
    def __init__(self):
        # Fluid properties (water at 20°C)
        self.density = 998.2      # kg/m³
        self.viscosity = 0.001    # Pa·s
        
        # Pipe properties
        self.length = 10.0        # m
        self.diameter = 0.05      # m
        self.roughness = 0.0015   # mm (typical for commercial steel)
        
    def reynolds_number(self, velocity):
        """Calculate Reynolds number"""
        return self.density * velocity * self.diameter / self.viscosity
    
    def friction_factor(self, reynolds):
        """
        Calculate Darcy friction factor using Colebrook-White equation
        """
        def colebrook(f):
            return (1/np.sqrt(f) + 2.0*np.log10(self.roughness/(3.7*self.diameter*1000) + 
                                               2.51/(reynolds*np.sqrt(f))))
        
        if reynolds < 2300:
            return 64/reynolds
        else:
            # Initial guess using Swamee-Jain equation
            f_0 = 0.25/(np.log10(self.roughness/(3.7*self.diameter*1000) + 5.74/reynolds**0.9)**2)
            return fsolve(colebrook, f_0)[0]
    
    def pressure_drop(self, velocity):
        """Calculate pressure drop using Darcy-Weisbach equation"""
        re = self.reynolds_number(velocity)
        f = self.friction_factor(re)
        return (f * self.length * self.density * velocity**2) / (2 * self.diameter)
    
    def velocity_profile(self, flow_rate, r_points=50):
        """Calculate velocity profile for laminar flow"""
        radius = self.diameter/2
        r = np.linspace(0, radius, r_points)
        
        # Maximum velocity at center
        v_max = 2 * flow_rate / (np.pi * radius**2)
        
        # Parabolic profile
        v = v_max * (1 - (r/radius)**2)
        return r, v
    
    def analyze_flow_range(self, velocities):
        """Analyze flow characteristics over a range of velocities"""
        reynolds = [self.reynolds_number(v) for v in velocities]
        pressure_drops = [self.pressure_drop(v) for v in velocities]
        friction_factors = [self.friction_factor(re) for re in reynolds]
        
        return reynolds, pressure_drops, friction_factors
    
    def plot_analysis(self, velocities):
        """Create visualizations of flow analysis"""
        reynolds, pressure_drops, friction_factors = self.analyze_flow_range(velocities)
        
        fig = plt.figure(figsize=(15, 10))
        
        # Pressure drop vs Velocity
        ax1 = fig.add_subplot(221)
        ax1.plot(velocities, pressure_drops, 'b-', linewidth=2)
        ax1.set_xlabel('Flow Velocity (m/s)')
        ax1.set_ylabel('Pressure Drop (Pa)')
        ax1.set_title('Pressure Drop vs Flow Velocity')
        ax1.grid(True)
        
        # Friction factor vs Reynolds number
        ax2 = fig.add_subplot(222)
        ax2.semilogx(reynolds, friction_factors, 'r-', linewidth=2)
        ax2.set_xlabel('Reynolds Number')
        ax2.set_ylabel('Friction Factor')
        ax2.set_title('Moody Diagram')
        ax2.grid(True)
        
        # Velocity profile for laminar flow
        ax3 = fig.add_subplot(223)
        flow_rate = 0.001  # m³/s
        r, v = self.velocity_profile(flow_rate)
        ax3.plot(r*1000, v, 'g-', linewidth=2)  # Convert radius to mm
        ax3.set_xlabel('Radial Position (mm)')
        ax3.set_ylabel('Velocity (m/s)')
        ax3.set_title('Laminar Flow Velocity Profile')
        ax3.grid(True)
        
        # Head loss vs Flow rate
        ax4 = fig.add_subplot(224)
        head_losses = [dp/(self.density * 9.81) for dp in pressure_drops]
        flow_rates = [v * np.pi * (self.diameter/2)**2 for v in velocities]
        ax4.plot(flow_rates, head_losses, 'm-', linewidth=2)
        ax4.set_xlabel('Flow Rate (m³/s)')
        ax4.set_ylabel('Head Loss (m)')
        ax4.set_title('Head Loss vs Flow Rate')
        ax4.grid(True)
        
        plt.tight_layout()
        return fig

# Create analyzer and run analysis
analyzer = PipeFlowAnalyzer()

# Define velocity range for analysis
velocities = np.linspace(0.1, 5, 100)

# Plot results
analyzer.plot_analysis(velocities)

# Print analysis for specific velocity
test_velocity = 1.0  # m/s
re = analyzer.reynolds_number(test_velocity)
dp = analyzer.pressure_drop(test_velocity)
f = analyzer.friction_factor(re)

print("Pipe Flow Analysis Results:")
print(f"\nTest conditions at velocity = {test_velocity} m/s:")
print(f"Reynolds Number: {re:.0f}")
print(f"Flow Regime: {'Laminar' if re < 2300 else 'Turbulent'}")
print(f"Friction Factor: {f:.4f}")
print(f"Pressure Drop: {dp/1000:.2f} kPa")
print(f"Head Loss: {dp/(analyzer.density * 9.81):.2f} m")

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

# Step 1: Create a social network graph
G = nx.erdos_renyi_graph(n=20, p=0.2, seed=42)  # 20 nodes, 20% connection probability

# Step 2: Initialize random opinions
opinions = np.random.rand(len(G.nodes))
initial_opinions = opinions.copy()

# Step 3: Define the adjacency matrix and normalize weights
adj_matrix = nx.to_numpy_array(G)
weights = adj_matrix / adj_matrix.sum(axis=1, keepdims=True)

# Step 4: Simulate opinion updates
iterations = 20
opinion_history = [opinions.copy()]

for _ in range(iterations):
    opinions = weights @ opinions  # Weighted average of neighbors
    opinion_history.append(opinions.copy())

# Step 5: Visualization of the results
plt.figure(figsize=(12, 6))

# Plot the convergence of opinions over time
for i, opinion_track in enumerate(np.array(opinion_history).T):
    plt.plot(opinion_track, label=f'Node {i}', alpha=0.6)

plt.title("Opinion Dynamics in a Social Network (DeGroot Model)")
plt.xlabel("Iteration")
plt.ylabel("Opinion Value")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize='small')
plt.grid(True)
plt.show()

# Plot the initial and final opinions on the graph
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Initial opinions
nx.draw(
    G, 
    with_labels=True, 
    node_color=initial_opinions, 
    cmap=plt.cm.viridis, 
    node_size=500, 
    ax=axes[0]
)
axes[0].set_title("Initial Opinions")

# Final opinions
nx.draw(
    G, 
    with_labels=True, 
    node_color=opinions, 
    cmap=plt.cm.viridis, 
    node_size=500, 
    ax=axes[1]
)
axes[1].set_title("Final Opinions")

plt.show()

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

# Define the Lotka-Volterra equations
def lotka_volterra(t, z, alpha, beta, delta, gamma):
    x, y = z
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# Parameters
alpha = 0.1   # Prey birth rate
beta = 0.02   # Predation rate
delta = 0.01  # Predator reproduction rate
gamma = 0.1   # Predator death rate

# Initial populations
x0 = 40  # Initial prey population
y0 = 9   # Initial predator population

# Time span for simulation
t_span = (0, 200)
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# Solve the equations
solution = solve_ivp(
    lotka_volterra, t_span, [x0, y0], args=(alpha, beta, delta, gamma), t_eval=t_eval
)

# Extract results
time = solution.t
prey = solution.y[0]
predators = solution.y[1]

# Plot population dynamics over time
plt.figure(figsize=(10, 6))
plt.plot(time, prey, label="Prey Population", color="blue")
plt.plot(time, predators, label="Predator Population", color="red")
plt.title("Lotka-Volterra Predator-Prey Dynamics")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend()
plt.grid()
plt.show()

# Phase plot
plt.figure(figsize=(8, 6))
plt.plot(prey, predators, color="purple")
plt.title("Phase Plot: Predator vs. Prey")
plt.xlabel("Prey Population")
plt.ylabel("Predator Population")
plt.grid()
plt.show()

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

def black_scholes_call(S, K, r, T, sigma):
    """
    Calculate the price of a European call option using the Black-Scholes formula.
    """
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    call_price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    return call_price

def visualize_black_scholes():
    # Parameters
    S = np.linspace(50, 150, 100)  # Stock price range
    K = 100                       # Strike price
    r = 0.05                      # Risk-free interest rate
    T = 1                         # Time to maturity (1 year)
    sigma = 0.2                   # Volatility (20%)
    
    # Calculate call option prices for each stock price
    call_prices = [black_scholes_call(s, K, r, T, sigma) for s in S]
    
    # Plot the results
    plt.figure(figsize=(10, 6))
    plt.plot(S, call_prices, label="Call Option Price")
    plt.axvline(K, color="red", linestyle="--", label="Strike Price (K)")
    plt.title("European Call Option Price (Black-Scholes Model)")
    plt.xlabel("Stock Price (S)")
    plt.ylabel("Call Option Price")
    plt.legend()
    plt.grid()
    plt.show()

if __name__ == "__main__":
    visualize_black_scholes()

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

def analyze_complex_function(z):
    """
    Analyzes the complex function f(z) = z^2 + 2z + 1
    Args:
        z (complex): Input complex number
    Returns:
        complex: Result of f(z)
    """
    return z**2 + 2*z + 1

# Create a grid of points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = X + Y*1j

# Calculate function values
W = analyze_complex_function(Z)

# Extract real and imaginary parts
real_part = W.real
imag_part = W.imag
abs_values = np.abs(W)

def plot_function_analysis():
    """
    Creates three subplots showing different aspects of the complex function
    """
    fig = plt.figure(figsize=(15, 5))
    
    # Real part
    ax1 = fig.add_subplot(131, projection='3d')
    surf1 = ax1.plot_surface(X, Y, real_part, cmap='viridis')
    ax1.set_title('Real Part')
    ax1.set_xlabel('Re(z)')
    ax1.set_ylabel('Im(z)')
    ax1.set_zlabel('Re(f(z))')
    
    # Imaginary part
    ax2 = fig.add_subplot(132, projection='3d')
    surf2 = ax2.plot_surface(X, Y, imag_part, cmap='viridis')
    ax2.set_title('Imaginary Part')
    ax2.set_xlabel('Re(z)')
    ax2.set_ylabel('Im(z)')
    ax2.set_zlabel('Im(f(z))')
    
    # Absolute value
    ax3 = fig.add_subplot(133, projection='3d')
    surf3 = ax3.plot_surface(X, Y, abs_values, cmap='viridis')
    ax3.set_title('Absolute Value')
    ax3.set_xlabel('Re(z)')
    ax3.set_ylabel('Im(z)')
    ax3.set_zlabel('|f(z)|')
    
    plt.tight_layout()
    return fig

# Plot the results
plot_function_analysis()

# Find critical points (where derivative is zero)
# For f(z) = z^2 + 2z + 1, f'(z) = 2z + 2
# Critical point is at z = -1

critical_point = -1
critical_value = analyze_complex_function(critical_point)
print(f"Critical point: z = {critical_point}")
print(f"Value at critical point: f({critical_point}) = {critical_value}")