Mathematics-I for CSE stream Lab Component

Prabodh C P

(Course Code: BMATS101) solutions using Python Programming

Hi everyone! This blog provides solutions for the lab component of the newly introduced subject “Mathematics-I for Computer Science and Engineering stream” (BMATS101) for I 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 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. 2D plots for Cartesian and polar curves
  2. Finding angle between polar curves, curvature and radius of curvature of a given curve
  3. Finding partial derivatives and Jacobian
  4. Applications to Maxima and Minima of two variables
  5. Solution of first-order ordinary differential equation and plotting the solution curves
  6. Finding GCD using Euclid’s Algorithm
  7. Solving linear congruences 𝒂𝒙 ≡ 𝒃(𝒎𝒐𝒅 𝒎)
  8. Numerical solution of system of linear equations, test for consistency and graphical representation
  9. Solution of system of linear equations using Gauss-Seidel iteration
  10. Compute Eigen values and Eigen vectors and find the largest and smallest Eigen value by Rayleigh power method.

2D plots for Cartesian and Polar curves

The following Python programs prompt the user to choose between Cartesian and Polar curves. It then generates those curves and plots the curve. We have taken two curves one in Cartesian and another in Polar form.

The Cartesian curve chosen for this example is $ y = \sin(2*x) * \cos(3*x) $

The Polar curve chosen for this example is $ r = a\cos(n\theta) $ with a = 1 and n = 6

Python Code

import matplotlib.pyplot as plt
import numpy as np

option = int(input("1.Cartesian Curve\n2.Polar Curve\nChoose the option : "))

if option == 1:
    # Define the x values
    x = np.linspace(-5, 5, 100)
    
    # Define the function to plot (y = sin(2x) + cos(3x))
    y = np.sin(2*x) * np.cos(3*x)
    #another example
    # y = x**3 + x**2
    # Create the plot
    plt.plot(x, y)
    
    # Add grid
    plt.grid(True)
    
    # Add labels and title
    plt.xlabel('X Axis')
    plt.ylabel('Y Axis')
    plt.title('Cartesian Curve')
    
    # Show the plot
    plt.show()
    
elif option == 2:
    # Define the theta values for the rose curve
    theta = np.linspace(0, 2*np.pi, 1000)
    
    # Define the parameters for the rose curve
    a = 1
    n = 6
    
    # Calculate the radius values for the rose curve
    r = a*np.cos(n*theta)
    
    # Create a polar plot of the rose curve
    plt.polar(theta, r)
    plt.title("Polar Curve")
    # Show the plot
    plt.show()

Output Sample 1

1.Cartesian Curve
2.Polar Curve
Choose the option : 1
1.Cartesian Curve

Output Sample 2

1.Cartesian Curve
2.Polar Curve
Choose the option : 2

Finding angle between polar curves, curvature and radius of curvature of a given curve

Finding angle between polar curves

Here we show how to find the angle between two polar curves, for this we consider the following two polar curves

  • $ r1 = 2\cos(\theta) $
  • $ r2 = \sin(2\theta) $

Python Code

import numpy as np
import matplotlib.pyplot as plt

# Define the first polar curve
def r1(theta):
    return 2 * np.cos(theta)

# Define the second polar curve
def r2(theta):
    return np.sin(2 * theta)

# Generate arrays of theta and r values for each curve
theta_range = np.linspace(0, 2 * np.pi, 1000)
r1_values = r1(theta_range)
r2_values = r2(theta_range)

# Calculate the angle between the curves
angle = np.arccos(np.abs(np.trapz(r1_values * r2_values, theta_range)) / (np.trapz(r1_values ** 2, theta_range) * np.trapz(r2_values ** 2, theta_range)))

# Print the angle in degrees
print("The angle between the polar curves is:", np.rad2deg(angle), "degrees")

# Plot the curves
plt.polar(theta_range, r1_values, label='r1')
plt.polar(theta_range, r2_values, label='r2')
plt.legend()
plt.show()

Output

The angle between the polar curves is: 90.0 degrees

Finding curvature and radius of curvature of a given curve

Here we will show how to find the curvature and radius of curvature for a given curve lets say $ y = \sin(x) $

Python Program

import numpy as np
import matplotlib.pyplot as plt

# Define the curve
x = np.linspace(-2.5, 2.5, 1000)
y = np.sin(x)

# Calculate the first and second derivatives
dydx = np.gradient(y, x)
d2ydx2 = np.gradient(dydx, x)

# Calculate the curvature
curvature = np.abs(d2ydx2 / (1 + dydx ** 2) ** 1.5)

# Calculate the radius of curvature
radius = 1 / curvature

#Plot the curve
fig, ax1 = plt.subplots()
ax1.plot(x, y, 'b-', label='curve')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.tick_params(axis='y')

# plt.plot(x, y)
# plt.grid(True)

# Plot the curvature
ax2 = ax1.twinx()
ax2.plot(x, curvature, 'r-', label='curvature')
ax2.set_ylabel('curvature')
ax2.tick_params(axis='y')

# Plot the radius of curvature
ax3 = ax1.twinx()
ax3.spines["right"].set_position(("axes", 1.2))
ax3.plot(x, radius, 'g-', label='radius')
ax3.set_ylabel('radius of curvature')
ax3.tick_params(axis='y')

# Combine the legends
lines = ax1.get_lines() + ax2.get_lines() + ax3.get_lines()
ax1.legend(lines, [line.get_label() for line in lines])

plt.grid(True)
plt.show()

Output

Finding partial derivatives and Jacobian

In this computation we define an equation of the form $ f(x,y,z) = x^2*y*z + y^2*z*x + z^2*x*y $. Now we first find the partial derivatives of $ f(x,y,z) $ with respect to x, y, and z and then its Jacobian.

import sympy as sp

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

# Define the function f(x, y, z)
f = x**2*y*z + y**2*z*x + z**2*x*y

# Find the partial derivatives of f with respect to x, y, and z
fx = sp.diff(f, x)
fy = sp.diff(f, y)
fz = sp.diff(f, z)

# Print the partial derivatives
print("Partial derivative of f with respect to x:", fx)
print("Partial derivative of f with respect to y:", fy)
print("Partial derivative of f with respect to z:", fz)
print()

# Find the Jacobian of f
J = sp.Matrix([fx, fy, fz])
vars = [x, y, z]
Jacobian = J.jacobian(vars)

# Print the Jacobian
for i in range(len(Jacobian)):
    print("{:25s}".format(str(Jacobian[i])),end = '\t')
    if (i+1)%3 == 0:
        print()

Output

The partial derivatives of $ f(x,y,z) = x^2*y*z + y^2*z*x + z^2*x*y $ and its Jacobian are shown below.

Partial derivative of f with respect to x: 2*x*y*z + y**2*z + y*z**2
Partial derivative of f with respect to y: x**2*z + 2*x*y*z + x*z**2
Partial derivative of f with respect to z: x**2*y + x*y**2 + 2*x*y*z

2*y*z                    	2*x*z + 2*y*z + z**2     	2*x*y + y**2 + 2*y*z     	
2*x*z + 2*y*z + z**2     	2*x*z                    	x**2 + 2*x*y + 2*x*z     	
2*x*y + y**2 + 2*y*z     	x**2 + 2*x*y + 2*x*z     	2*x*y

Solution of first-order ordinary differential equation and plotting the solution curves

The differential equation considered for this program is $ f(t,y) : t * y – y^2 $

Python Code

import numpy as np
import matplotlib.pyplot as plt

# Define the differential equation
def f(t, y):
    return t * y - y**2

# Define the time interval and initial condition
t0 = 0
tf = 5
y0 = 1

# Define the step size
h = 0.01

# Initialize the solution array
t = np.arange(t0, tf+h, h)
y = np.zeros(len(t))
y[0] = y0

# Solve the differential equation using Euler's method
for i in range(len(t)-1):
    y[i+1] = y[i] + h * f(t[i], y[i])

# Plot the solution curve
plt.plot(t, y)
plt.xlabel('t')
plt.ylabel('y')
plt.title('Solution Curve for dy/dt = ty - y^2')
plt.show()

Output

Finding GCD using Euclid’s Algorithm

We write a python program to find the Greatest Common Divisor of two non-negative numbers m and n using the following definition

Euclid’s Algorithm

Python Code

def gcd(a, b):
    """
    Returns the Greatest Common Divisor (GCD) of two integers using Euclid's Algorithm.
    """
    while b != 0:
        a, b = b, a % b
    return a

m = int(input("Enter first number : "))
n = int(input("Enter second number : "))

res = gcd(m,n)

print("GCD of", m, "and", n, "is", res)

Output

Sample 1

Enter first number : 72
Enter second number : 48
GCD of 72 and 48 is 24

Sample 2

Enter first number : 7
Enter second number : 0
GCD of 7 and 0 is 7

Solving linear congruence of the form 𝒂𝒙 ≡ 𝒃(𝒎𝒐𝒅 𝒎)

The following program calculates the near congruence of the form 𝒂𝒙 ≡ 𝒃(𝒎𝒐𝒅 𝒎) if the solution exists. To solve linear congruence we can use the extended Euclidean algorithm. The extended Euclidean algorithm computes the greatest common divisor of two integers a and b, along with coefficients x and y that satisfy the equation ax + by = gcd(a, b).

If gcd(a, m) does not divide b, then there is no solution to the congruence. Otherwise, we can divide both sides of the congruence by gcd(a, m), giving us a new congruence 𝑎′𝑥≡𝑏′(mod 𝑚′), where 𝑎′=𝑎/gcd(𝑎,𝑚), 𝑏′=𝑏/gcd(𝑎,𝑚), and 𝑚′=𝑚/gcd(𝑎,𝑚).

Python Program

def linear_congruence(a, b, m):
    # check if gcd(a, m) divides b
    gcd, x, y = extended_gcd(a, m)
    if b % gcd != 0:
        return None  # no solution
    else:
        # reduce the congruence
        a //= gcd
        b //= gcd
        m //= gcd
        # find the solution
        x = (x * b) % m
        return x % m

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0
    else:
        gcd, x, y = extended_gcd(b, a % b)
        return gcd, y, x - (a // b) * y


def main():
    a =  int(input("Enter value for a : "))    
    b =  int(input("Enter value for b : "))    
    m =  int(input("Enter value for m : "))  
    
    if linear_congruence(a, b, m) == None:
        print("There is no solution")
    else:
        print("linear_congruence for", (a, b, m), "is", linear_congruence(a, b, m))
        
main()

Output

Sample 1

Enter value for a : 34
Enter value for b : 105
Enter value for m : 55
linear_congruence for (34, 105, 55) is 50

Sample 2

Enter value for a : 35
Enter value for b : 49
Enter value for m : 25
There is no solution

Numerical solution of system of linear equations, test for consistency and graphical representation

The following program takes the following system of linear equations for solving

\[2x+3y=8\]

\[4x+9y=15\]

This can be represented in matrix form as follows

To test for consistency, we can check whether the determinant of the matrix A is nonzero. If the determinant is zero, then the system of equations has no unique solution. Since the determinant of A in our example is 2*9 - 3*4 = 6, which is nonzero, we know that the system of equations is consistent.

Python Code

import numpy as np
import matplotlib.pyplot as plt

"""
Lets solve
2x + 3y = 8
4x + 9y = 15
"""

A = np.array([[2, 3], [4, 9]])
b = np.array([8, 15])

x = np.linalg.solve(A, b)

if np.linalg.det(A) != 0:
    print("The system of equations is consistent.")
else:
    print("The system of equations is inconsistent.")

# Plot the two equations as lines
x_vals = np.linspace(-8, 8, 100)
y1_vals = (8 - 2*x_vals)/3
y2_vals = (15 - 4*x_vals)/9
plt.plot(x_vals, y1_vals, label="2x + 3y = 8")
plt.plot(x_vals, y2_vals, label="4x + 9y = 15")

# Plot the intersection point as a red dot
plt.plot(x[0], x[1], 'ro', label="Intersection")

#show grid
plt.grid(True)
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("System of Linear Equations")
plt.show()

Output

We can graphically represent the system of equations using Matplotlib. We’ll plot the two equations as lines and show the intersection point as a red dot.

The system of equations is consistent.

Solution of system of linear equations using Gauss-Seidel iteration

We will define our system of linear equations. For example, let’s consider the following system of equations:

3x + y + z = 1
x + 4y + z = 4
x + y + 5z = 5

We can represent this system of equations in matrix form as:

[3 1 1] [x]   [1]
[1 4 1] [y] = [4]
[1 1 5] [z]   [5]

We then define the Gauss-Seidel iteration method. This method iteratively computes a new estimate of the solution vector x using the previous estimate, until the error is within a specified tolerance.

Python Code

import numpy as np

def gauss_seidel(A, b, x0, tol=1e-10, max_iter=1000):
    """
    Solves a system of linear equations Ax=b using the Gauss-Seidel iteration method.

    Parameters:
    A (array): the matrix of coefficients
    b (array): the right-hand side vector
    x0 (array): the initial guess for the solution vector
    tol (float): the tolerance for the error (default: 1e-10)
    max_iter (int): the maximum number of iterations (default: 1000)

    Returns:
    x (array): the solution vector
    iter (int): the number of iterations performed
    """
    n = len(b)
    x = np.copy(x0)
    iter = 0
    while iter < max_iter:
        x_new = np.zeros(n)
        for i in range(n):
            s1 = np.dot(A[i, :i], x_new[:i])
            s2 = np.dot(A[i, i + 1:], x[i + 1:])
            x_new[i] = (b[i] - s1 - s2) / A[i, i]
        if np.linalg.norm(x_new - x) < tol:
            return x_new, iter
        x = np.copy(x_new)
        iter += 1
    return x, iter


"""
Lets solve
3x + y + z = 1
x + 4y + z = 4
x + y + 5z = 5
"""

A = np.array([[3, 1, 1], [1, 4, 1], [1, 1, 5]])
b = np.array([1, 4, 5])

x0 = np.zeros(len(b))
x, iter = gauss_seidel(A, b, x0)
print("Solution:", x)
print("Iterations:", iter)

Output

Solution: [-0.24  0.84  0.88]
Iterations: 13

Compute Eigen values and Eigen vectors and find the largest and smallest Eigen value by Rayleigh power method.

First we’ll define our matrix for which we want to compute the eigenvalues and eigenvectors. Let’s consider the following example.

      [3 1 1]
  A = [1 4 1]
      [1 1 5]

We can use the NumPy function linalg.eig to compute the eigenvalues and eigenvectors of A:

w, v = np.linalg.eig(A)

The variable w contains the eigenvalues of A, and the variable v contains the corresponding eigenvectors. The i-th eigenvalue in w corresponds to the i-th eigenvector in v.

To find the largest and smallest eigenvalues of A using the Rayleigh power method, we first need to choose an initial guess for the corresponding eigenvectors. For example, we can choose the first column of A as the initial guess for the eigenvector corresponding to the largest eigenvalue, and the last column of A as the initial guess for the eigenvector corresponding to the smallest eigenvalue:

x = np.array([1, 0, 0])  # initial guess for largest eigenvalue
y = np.array([0, 0, 1])  # initial guess for smallest eigenvalue

Next the Rayleigh power method iteratively computes a new estimate of the eigenvalue and eigenvector using the previous estimate, until the error is within a specified tolerance. Lets now look at the complete program.

Python Code

import numpy as np

def rayleigh_power(A, y, max_iter=100, tol=1e-6):
    """Compute the largest eigenvalue and eigenvector of a matrix A using the Rayleigh power method."""
    iter = 0
    lam = 0
    y = y / np.linalg.norm(y)  # normalize y
    while iter < max_iter:
        y_old = y
        lam_old = lam
        y = np.dot(A, y)
        lam = np.dot(y, y_old)
        y = y / np.linalg.norm(y)  # normalize y
        if np.abs(lam - lam_old) < tol:
            break
        iter += 1
    return lam, y, iter

def inverse_rayleigh_power(A, y, max_iter=100, tol=1e-6):
    """Compute the smallest eigenvalue and eigenvector of a matrix A using the inverse Rayleigh power method."""
    iter = 0
    lam = 0
    y = y / np.linalg.norm(y)  # normalize y
    while iter < max_iter:
        y_old = y
        lam_old = lam
        y = np.linalg.solve(A, y)
        lam = np.dot(y, y_old)
        y = y / np.linalg.norm(y)  # normalize y
        if np.abs(lam - lam_old) < tol:
            break
        iter += 1
    lam = 1/lam
    return lam, y, iter

# Define the matrix A and the initial guess for the eigenvector
A = np.array([[3, -1, 1], [-1, 3, 1], [1, 1, 5]])
y = np.array([1, 1, 1])

# Use the Rayleigh power method to find the largest eigenvalue and eigenvector
lam1, v, iter1 = rayleigh_power(A, y)
print("Largest eigenvalue of A:", lam1)
print()

# Use the inverse Rayleigh power method to find the smallest eigenvalue and eigenvector
lam2, y1, iter2 = inverse_rayleigh_power(A, y)
print("Smallest eigenvalue of A:", lam2)
print()

# Print the eigenvectors corresponding to the eigenvalues
print("Eigenvectors of A:")
print(v)
print()

# Print the number of iterations required for each method
print("Number of iterations for Rayleigh power method:", iter1)
print("Number of iterations for inverse Rayleigh power method:", iter2)

This code first uses the Rayleigh power method to find the largest eigenvalue and eigenvector of the matrix A. It then uses the inverse Rayleigh power method to find the smallest eigenvalue and eigenvector. Finally, it prints the eigenvalues and eigenvectors of the matrix, as well as the number of iterations required for each method.

Output

Largest eigenvalue of A: 5.561552802322129

Smallest eigenvalue of A: 1.4384472022588988

Eigenvectors of A:
[0.26096505 0.26096505 0.92940545]

Number of iterations for Rayleigh power method: 7
Number of iterations for inverse Rayleigh power method: 7

If you are also looking for other subject Lab Manuals, head over to my following blog :

VTU Lab Manuals using FOSS

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

Related Posts

Leave a Reply

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