Mathematics-II for CSE stream Lab Component

(Course Code: BMATS201) solutions using Python Programming

Hi everyone! This blog provides solutions for the lab component of the newly introduced subject “Mathematics-II for Computer Science and Engineering stream” (BMATS201) for II semester students of VTU. I have provided solutions for the lab exercises using Python. All the programs shown below have been tested on Python 3.11.7 on Ubuntu 22.04 OS using Anaconda Distribution.

Installation Instructions

A detailed procedure of installing ANACONDA distribution and necessary tools for this is provided in my other blog listed below. https://moodle.sit.ac.in/blog/setting-up-anaconda-python-programming-environment/

Questions

  1. Program to compute area, surface area, volume and centre of gravity
  2. Evaluation of improper integrals
  3. Finding gradient, divergent, curl and their geometrical interpretation
  4. Basis and Dimension for a vector space and Graphical representation of linear transformation
  5. Computing the inner product and orthogonality
  6. Solution of algebraic and transcendental equations by Ramanujan’s, Regula-Falsi and Newton-Raphson method
  7. Interpolation/Extrapolation using Newton’s forward and backward difference formula
  8. Computation of area under the curve using Trapezoidal, Simpson’s (1/3)rd and (3/8)th rule
  9. Solution of ODE of first order and first degree by Taylor’s series and Modified Euler’s method
  10. Solution of ODE of first order and first degree by Runge-Kutta 4th order and Milne’s predictor-corrector method

1) Program to compute area, surface area, volume and centre of gravity

Below is a Python program that computes the area, surface area, volume, and center of gravity of a given 3D shape. This example specifically targets a sphere, but you can adapt it for other shapes by defining appropriate functions for area, surface area, volume, and center of gravity calculations.

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

To run this program online click n the link below

P01 Program

Output Sample

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

Additional code

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

Enter the radius of the sphere: 6



2) Evaluation of improper integrals

Python Code

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

To run this program online click n the link below

P02 Program

Output

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

3) Finding gradient, divergent, curl and their geometrical interpretation

Python Program

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.

To run this program online click n the link below

P03 Program

Output

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

4)Computation of basis and dimension for a vector space and Graphical representation of linear transformation

Python Code

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)

To run this program online click n the link below

P04 Program

Output

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

5) Computing the inner product and orthogonality

Python Code

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

To run this program online click n the link below

P05 Program

Output

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.

6) Solution of algebraic and transcendental equations by Ramanujan’s, Regula-Falsi and Newton-Raphson method

Python Code

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)

To run this program online click n the link below

P06 Program

Output

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

7) Interpolation/Extrapolation using Newton’s forward and backward difference formula

Python Program

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)

To run this program online click n the link below

P07 Program

Output


8) Computation of area under the curve using Trapezoidal, Simpson’s (1/3)rd and (3/8)th rule

Python Code

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)

To run this program online click n the link below

P08 Program

Output



9) Solution of ODE of first order and first degree by Taylor’s series and Modified Euler’s method

Python Code

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

To run this program online click n the link below

P09 Program

Output


10) Solution of ODE of first order and first degree by Runge-Kutta 4th order and Milne’s predictor-corrector method

Python Code

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

To run this program online click n the link below

P10 Program

Output

Prabodh C P is a faculty in the Dept of CSE SIT, Tumkur and also currently a Research Scholar pursuing PhD in IIT Hyderabad. He conducts online classes for C, C++, Python. For more info call +919392302100

Leave a Reply

Your email address will not be published. Required fields are marked *