Mathematics-II for CSE stream Lab Component

import math

def sphere_area(radius):
    return 4 * math.pi * radius ** 2

def sphere_surface_area(radius):
    return 4 * math.pi * radius ** 2

def sphere_volume(radius):
    return (4/3) * math.pi * radius ** 3

def sphere_center_of_gravity(radius):
    return (3/8) * radius

def main():
    radius = float(input("Enter the radius of the sphere: "))

    area = sphere_area(radius)
    surface_area = sphere_surface_area(radius)
    volume = sphere_volume(radius)
    center_of_gravity = sphere_center_of_gravity(radius)

    print(f"Area of the sphere: {area}")
    print(f"Surface area of the sphere: {surface_area}")
    print(f"Volume of the sphere: {volume}")
    print(f"Center of gravity of the sphere: {center_of_gravity}")

main()

Enter the radius of the sphere:  6

Area of the sphere: 452.3893421169302
Surface area of the sphere: 452.3893421169302
Volume of the sphere: 904.7786842338603
Center of gravity of the sphere: 2.25

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

def plot_wireframe_sphere(radius):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)
    x = radius * np.outer(np.cos(u), np.sin(v))
    y = radius * np.outer(np.sin(u), np.sin(v))
    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))

    ax.plot_wireframe(x, y, z, color='b')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Wireframe Sphere')

    plt.show()

def main():
    radius = float(input("Enter the radius of the sphere: "))
    plot_wireframe_sphere(radius)

main()

from scipy import integrate

# Define the function to be integrated
def integrand(x):
    return (x**3 + 2*x + 7) /(x + 2)

# Define the limits of integration
a = 0
b = 1

# Evaluate the improper integral with increased limit
result, error = integrate.quad(integrand, a, b)

print(f"Result of the improper integral: {result}")
print(f"Estimated error: {error}")

Result of the improper integral: 3.3060077927925113
Estimated error: 3.6704059711484397e-14

import numpy as np
import sympy as sp

# Define symbols
x, y, z = sp.symbols('x y z')

# Define vector field
F = sp.Matrix([x**2, y**2, z**2])

# Define variables for differentiation
variables = (x, y, z)

# Define gradient
gradient_F = sp.Matrix([sp.diff(F[i], var) for var in variables for i in range(3)]).reshape(3, 3)

# Calculate divergence
divergence_F = sp.diff(F[0], x) + sp.diff(F[1], y) + sp.diff(F[2], z)

# Calculate curl
curl_F = sp.Matrix([sp.diff(F[(i+1) % 3], variables[(i+2) % 3]) - sp.diff(F[(i+2) % 3], variables[(i+1) % 3])
                    for i in range(3)])

# Print results
print("Vector Field F(x, y, z) = ", F)
print("\nGradient of F:")
print(gradient_F)
print("\nDivergence of F:")
print(divergence_F)
print("\nCurl of F:")
print(curl_F)

# Geometrical interpretations
# Gradient gives the direction of maximum change in the scalar field.
# Divergence represents how much a vector field is spreading out or converging at a point.
# Curl indicates the rotation or "whirl" of the vector field at a point.

Vector Field F(x, y, z) =  Matrix([[x**2], [y**2], [z**2]])

Gradient of F:
Matrix([[2*x, 0, 0], [0, 2*y, 0], [0, 0, 2*z]])

Divergence of F:
2*x + 2*y + 2*z

Curl of F:
Matrix([[0], [0], [0]])

import numpy as np
import matplotlib.pyplot as plt

def basis_and_dimension(matrix):
    # Compute the reduced row echelon form
    rref_matrix = np.linalg.matrix_rank(matrix)
    dimension = np.linalg.matrix_rank(matrix.T)
    
    return rref_matrix, dimension

def plot_linear_transformation(matrix):
    # Create a set of 2D vectors
    vectors = np.array([[1, 0], [0, 1]])
    
    # Apply linear transformation to the vectors
    transformed_vectors = np.dot(matrix, vectors)
    
    # Plot original and transformed vectors
    plt.figure(figsize=(8, 6))
    plt.quiver([0, 0], [0, 0], vectors[:,0], vectors[:,1], angles='xy', scale_units='xy', scale=1, color=['r', 'b'], label='Original Vectors')
    plt.quiver([0, 0], [0, 0], transformed_vectors[:,0], transformed_vectors[:,1], angles='xy', scale_units='xy', scale=1, color=['g', 'y'], label='Transformed Vectors')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.axhline(0, color='black',linewidth=0.5)
    plt.axvline(0, color='black',linewidth=0.5)
    plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
    plt.axis('equal')
    plt.legend()
    plt.title('Linear Transformation')
    plt.show()

# Example matrix representing a linear transformation
A = np.array([[1, 2], [3, 4]])

# Compute basis and dimension
rref, dimension = basis_and_dimension(A)
print("Reduced Row Echelon Form:\n", rref)
print("Dimension of Column Space:", dimension)

# Plot linear transformation
plot_linear_transformation(A)

Reduced Row Echelon Form:
 2
Dimension of Column Space: 2

import numpy as np

def inner_product(vector1, vector2):
    # Compute the inner product (dot product) between two vectors
    return np.dot(vector1, vector2)

def check_orthogonality(vector1, vector2):
    # Check if two vectors are orthogonal (their inner product is zero)
    return inner_product(vector1, vector2) == 0

# Example vectors
vector1 = np.array([1, 2, 3])
vector2 = np.array([-1, 1, 0])

# Compute the inner product of the vectors
result = inner_product(vector1, vector2)
print("Inner product of vector1 and vector2:", result)

# Check orthogonality
if check_orthogonality(vector1, vector2):
    print("Vectors vector1 and vector2 are orthogonal.")
else:
    print("Vectors vector1 and vector2 are not orthogonal.")
    
# Example vectors
vector3 = np.array([1, 2, 3])
vector4 = np.array([-1, 1, 2])

# Compute the inner product of the vectors
result = inner_product(vector3, vector4)
print("Inner product of vector3 and vector4:", result)

# Check orthogonality
if check_orthogonality(vector3, vector4):
    print("Vectors vector3 and vector4 are orthogonal.")
else:
    print("Vectors vector3 and vector4 are not orthogonal.")

Inner product of vector1 and vector2: 1
Vectors vector1 and vector2 are not orthogonal.

Inner product of vector3 and vector4: 7
Vectors vector3 and vector4 are not orthogonal.

import numpy as np

def ramanujan(f, a, b, tol=1e-6, max_iter=100):
    # Ramanujan's method for solving equations
    for _ in range(max_iter):
        c = (a * f(b) - b * f(a)) / (f(b) - f(a))
        if np.abs(f(c)) < tol:
            return c
        if f(a) * f(c) < 0:
            b = c
        else:
            a = c
    raise ValueError("Ramanujan's method did not converge")

def regula_falsi(f, a, b, tol=1e-6, max_iter=100):
    # Regula-Falsi method for solving equations
    for _ in range(max_iter):
        c = (a * f(b) - b * f(a)) / (f(b) - f(a))
        if np.abs(f(c)) < tol:
            return c
        if f(a) * f(c) < 0:
            b = c
        else:
            a = c
    raise ValueError("Regula-Falsi method did not converge")

def newton_raphson(f, f_prime, x0, tol=1e-6, max_iter=100):
    # Newton-Raphson method for solving equations
    x = x0
    for _ in range(max_iter):
        x_new = x - f(x) / f_prime(x)
        if np.abs(f(x_new)) < tol:
            return x_new
        x = x_new
    raise ValueError("Newton-Raphson method did not converge")

# Example usage:

# Define the function
def f(x):
    return x**3 - 2*x - 5

# Define the derivative of the function
def f_prime(x):
    return 3*x**2 - 2

# Initial guesses
x0_ramanujan = 1
x0_regula_falsi = 1
x0_newton_raphson = 1

# Solve using Ramanujan's method
solution_ramanujan = ramanujan(f, 1, 3)
print("Solution using Ramanujan's method:", solution_ramanujan)

# Solve using Regula-Falsi method
solution_regula_falsi = regula_falsi(f, 1, 3)
print("Solution using Regula-Falsi method:", solution_regula_falsi)

# Solve using Newton-Raphson method
solution_newton_raphson = newton_raphson(f, f_prime, 1)
print("Solution using Newton-Raphson method:", solution_newton_raphson)

Solution using Ramanujan's method: 2.094551400063298
Solution using Regula-Falsi method: 2.094551400063298
Solution using Newton-Raphson method: 2.0945514815642077

import numpy as np
import matplotlib.pyplot as plt

def forward_difference_coefficients(y):
    n = len(y)
    coeffs = np.zeros((n, n))
    coeffs[:, 0] = y
    
    for j in range(1, n):
        for i in range(n - j):
            coeffs[i, j] = coeffs[i+1, j-1] - coeffs[i, j-1]
            
    return coeffs

def backward_difference_coefficients(y):
    n = len(y)
    coeffs = np.zeros((n, n))
    coeffs[:, 0] = y
    
    for j in range(1, n):
        for i in range(j, n):
            coeffs[i, j] = coeffs[i, j-1] - coeffs[i-1, j-1]
            
    return coeffs

def newton_interpolation(x, y, coeffs, xi):
    n = len(x)
    h = x[1] - x[0]
    
    fxi = coeffs[0, 0]
    for j in range(1, n):
        term = coeffs[0, j]
        for i in range(j):
            term *= (xi - x[i]) / ((j+1) * h)
        fxi += term
        
    return fxi

# Data points
x = np.array([0, 1, 2, 3, 4])
y = np.array([1, 2, 4, 8, 16])

# Coefficients for forward and backward difference formula
coeffs_forward = forward_difference_coefficients(y)
coeffs_backward = backward_difference_coefficients(y)

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

# Forward Difference Coefficients
plt.subplot(1, 2, 1)
plt.imshow(coeffs_forward, cmap='coolwarm', aspect='auto')
plt.title('Forward Difference Coefficients')
plt.xlabel('Order of Difference')
plt.ylabel('Data Points')
plt.colorbar(label='Coefficient Value')

# Backward Difference Coefficients
plt.subplot(1, 2, 2)
plt.imshow(coeffs_backward, cmap='coolwarm', aspect='auto')
plt.title('Backward Difference Coefficients')
plt.xlabel('Order of Difference')
plt.ylabel('Data Points')
plt.colorbar(label='Coefficient Value')

plt.tight_layout()
plt.show()

# Interpolation
xi = 2.5
forward_interpolation = newton_interpolation(x, y, coeffs_forward, xi)
backward_interpolation = newton_interpolation(x, y, coeffs_backward, xi)

print("Interpolated value at xi using forward difference formula:", forward_interpolation)
print("Interpolated value at xi using backward difference formula:", backward_interpolation)

import numpy as np
import matplotlib.pyplot as plt

def trapezoidal_rule(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    area = h * (np.sum(y) - 0.5 * (y[0] + y[-1]))
    return area, x, y

def simpsons_one_third_rule(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    area = h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-2:2]) + y[-1])
    return area, x, y

def simpsons_three_eighth_rule(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = f(x)
    area = 3*h/8 * (y[0] + 3*np.sum(y[1:-1:3]) + 3*np.sum(y[2:-2:3]) + 2*np.sum(y[3:-3:3]) + y[-1])
    return area, x, y

# Example function
def f(x):
    return np.sin(x)

# Integration interval
a = 0
b = np.pi

# Number of subintervals
n = 10

# Compute area using different rules
area_trapezoidal, x_trap, y_trap = trapezoidal_rule(f, a, b, n)
area_simpsons_one_third, x_simp13, y_simp13 = simpsons_one_third_rule(f, a, b, n)
area_simpsons_three_eighth, x_simp38, y_simp38 = simpsons_three_eighth_rule(f, a, b, n)

# Plot the function
x = np.linspace(a, b, 100)
y = f(x)

plt.figure(figsize=(10, 6))
plt.plot(x, y, label='Function: sin(x)')
plt.fill_between(x_trap, 0, y_trap, alpha=0.3, label='Trapezoidal Rule Area')
plt.fill_between(x_simp13, 0, y_simp13, alpha=0.3, label="Simpson's 1/3 Rule Area")
plt.fill_between(x_simp38, 0, y_simp38, alpha=0.3, label="Simpson's 3/8 Rule Area")
plt.xlabel('x')
plt.ylabel('y')
plt.title('Area Under the Curve Approximation')
plt.legend()
plt.grid(True)
plt.show()

# Print results
print("Area under the curve using Trapezoidal rule:", area_trapezoidal)
print("Area under the curve using Simpson's 1/3 rule:", area_simpsons_one_third)
print("Area under the curve using Simpson's 3/8 rule:", area_simpsons_three_eighth)

import numpy as np
import matplotlib.pyplot as plt

def taylors_series(f, y0, t0, tn, h):
    t_values = np.arange(t0, tn+h, h)
    y_values = np.zeros_like(t_values)
    y_values[0] = y0
    
    for i in range(1, len(t_values)):
        t = t_values[i-1]
        y = y_values[i-1]
        y_prime = f(t, y)
        y_double_prime = f(t, y_prime)
        y_values[i] = y + h*y_prime + (h**2 / 2)*y_double_prime
    
    return t_values, y_values

def modified_euler(f, y0, t0, tn, h):
    t_values = np.arange(t0, tn+h, h)
    y_values = np.zeros_like(t_values)
    y_values[0] = y0
    
    for i in range(1, len(t_values)):
        t = t_values[i-1]
        y = y_values[i-1]
        y_prime = f(t, y)
        y_modified = y + h*y_prime
        y_prime_modified = f(t+h, y_modified)
        y_values[i] = y + (h/2)*(y_prime + y_prime_modified)
    
    return t_values, y_values

# Example ODE: y' = y - t^2 + 1, y(0) = 0.5
def f(t, y):
    return y - t**2 + 1

# Initial conditions
y0 = 0.5
t0 = 0
tn = 2
h = 0.2

# Solve using Taylor's series method
t_values_taylor, y_values_taylor = taylors_series(f, y0, t0, tn, h)

# Solve using Modified Euler's method
t_values_euler, y_values_euler = modified_euler(f, y0, t0, tn, h)

# Plot the solutions
plt.figure(figsize=(10, 6))
plt.plot(t_values_taylor, y_values_taylor, label="Taylor's Series")
plt.plot(t_values_euler, y_values_euler, label="Modified Euler's Method")
plt.xlabel('t')
plt.ylabel('y')
plt.title('Solution of ODE: y\' = y - t^2 + 1')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def runge_kutta_4th_order(f, y0, t0, tn, h):
    t_values = np.arange(t0, tn+h, h)
    y_values = np.zeros_like(t_values)
    y_values[0] = y0
    
    for i in range(1, len(t_values)):
        t = t_values[i-1]
        y = y_values[i-1]
        
        k1 = h * f(t, y)
        k2 = h * f(t + h/2, y + k1/2)
        k3 = h * f(t + h/2, y + k2/2)
        k4 = h * f(t + h, y + k3)
        
        y_values[i] = y + (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t_values, y_values

def milnes_predictor_corrector(f, y0, t0, tn, h):
    t_values = np.arange(t0, tn+h, h)
    y_values = np.zeros_like(t_values)
    y_values[0] = y0
    
    for i in range(1, 4):
        t = t_values[i-1]
        y = y_values[i-1]
        k1 = h * f(t, y)
        k2 = h * f(t + h/2, y + k1/2)
        k3 = h * f(t + h/2, y + k2/2)
        k4 = h * f(t + h, y + k3)
        
        y_values[i] = y + (k1 + 2*k2 + 2*k3 + k4) / 6
    
    for i in range(3, len(t_values)):
        t_n3 = t_values[i-3]
        t_n2 = t_values[i-2]
        t_n1 = t_values[i-1]
        y_n3 = y_values[i-3]
        y_n2 = y_values[i-2]
        y_n1 = y_values[i-1]
        
        t = t_values[i-1]
        y_p = y_n1 + 4*h/3 * (2*f(t, y_n1) - f(t_n1, y_n2) + 2*f(t_n2, y_n3))
        t = t_values[i]
        y_values[i] = y_n1 + h/3 * (f(t, y_p) + 4*f(t_n1, y_n1) + f(t_n2, y_n2))
    
    return t_values, y_values

# Example ODE: y' = y - t^2 + 1, y(0) = 0.5
def f(t, y):
    return y - t**2 + 1

# Initial conditions
y0 = 0.5
t0 = 0
tn = 2
h = 0.2

# Solve using Runge-Kutta 4th order method
t_values_rk4, y_values_rk4 = runge_kutta_4th_order(f, y0, t0, tn, h)

# Solve using Milne's predictor-corrector method
t_values_milne, y_values_milne = milnes_predictor_corrector(f, y0, t0, tn, h)

# Plot the solutions
plt.figure(figsize=(10, 6))
plt.plot(t_values_rk4, y_values_rk4, label="Runge-Kutta 4th Order")
plt.plot(t_values_milne, y_values_milne, label="Milne's Predictor-Corrector")
plt.xlabel('t')
plt.ylabel('y')
plt.title('Solution of ODE: y\' = y - t^2 + 1')
plt.legend()
plt.grid(True)
plt.show()