BMATS201 LAB PROGRAMMES
Part 1: Double and triple integration
1. Evaluate the integral ∫ ∫
from sympy import *
x,y=symbols('x y ')
w1=integrate(x**2+y**2,(y,0,x),(x,0,1))
print(w1)
Ans: 1/3
2. Evaluate the integral ∫ ∫
from sympy import *
x,y=symbols('x y ')
w1=integrate(x*y,(y,x,x**2),(x,0,1))
print(w1)
Ans: -1/24
3. Evaluate the integral ∫ ∫
from sympy import *
x,y=symbols('x y ')
w1=integrate(x**2+y**2,(y,x, sqrt(x)),(x,0,1))
print(w1)
Ans: 3/35
4. Evaluate the integral ∫ ∫ ∫
from sympy import *
x,y,z=symbols('x y z')
w2=integrate((x*y*z),(z,0,3-x-y),(y,0,3-x),(x,0,3))
print(w2)
Ans: 81/80
5. Prove that ∫ ∫ ∫∫
from sympy import *
x,y,z=symbols('x y z')
w3=integrate(x**2+y**2,y,x)
print(w3)
w4=integrate(x**2+y**2,x,y)
print(w4)
Ans
6. Find the area of an ellipse by double integration
√
∫ ∫
from sympy import *
x=Symbol('x')
y=Symbol('y')
a=Symbol('a')
b=Symbol('b')
a=4
b=6
w3=4*integrate(1,(y,0,(b/a)*sqrt(a**2-x**2)),(x,0,a))
print(w3)
Ans:
7. Find the area of the cardioid by double integration
from sympy import *
r=Symbol('r')
t=Symbol('t')
a=Symbol('a')
a=4
w3=2*integrate(r,(r,0,a*(1+cos(t))),(t,0,pi))
print(w3)
Ans:
Part 2: Gradient. Curl, Divergence
1. To find gradient of
from sympy.vector import *
from sympy import symbols
N=CoordSys3D('N') #Setting the coordinate system
x,y,z=symbols('x y z')
A=N.x**2*N.y+2*N.x*N.z-4 #Variables x,y,z to be used with coordinate system N
delop=Del() #Del operator
print(delop(A)) #Del operator applied to A
gradA=gradient(A) #Gradient function is used
print(f"\n Gradient of {A} is \n")
print(gradA)
Ans:
2. To find divergence of ⃗
from sympy.vector import *
from sympy import symbols
N=CoordSys3D('N')
x,y,z=symbols('x y z')
A=N.x**2*N.y*N.z*N.i+N.y**2*N.z*N.x*N.j+N.z**2*N.x*N.y*N.k
delop=Del()
divA=delop.dot(A)
print(divA)
print(f"\n Divergence of {A} is \n")
print(divergence(A))
Ans: 6xyz
3. To find curl of ⃗
from sympy.vector import *
from sympy import symbols
N=CoordSys3D('N')
x,y,z=symbols('x y z')
A=N.x**2*N.y*N.z*N.i+N.y**2*N.z*N.x*N.j+N.z**2*N.x*N.y*N.k
delop=Del()
curlA=delop.cross(A)
print(curlA)
print(f"\n Curl of {A} is \n ")
print(curl(A))
Ans:
PART 3: Beta Gamma functions and Vector space
1.Evaluate Gamma(5) by using definition
from sympy import *
x=symbols('x')
w1=integrate(exp(-x)*x**4,(x,0,float('inf')))
print(simplify(w1))
Ans: 24
2. Find Beta(3,5), Gamma(5)
from sympy import beta, gamma
m=input('m :');
n=input('n :');
m=float(m);
n=float(n);
s=beta(m,n);
t=gamma(n)
print('gamma (',n,') is %3.3f't)
print('Beta (',m,n,') is %3.3f's)
Ans:
m :3
n :5
gamma ( 5.0 ) is 24.000
Beta ( 3.0 5.0 ) is 0.010
3. Calculate Beta(5/2,7/2) and Gamma(5/2)
from sympy import beta, gamma
m=float(input('m : '));
n=float(input('n :'));
s=beta(m,n);
t=gamma(n)
print('gamma (',n,') is %3.3f'%t)
print('Beta (',m,n,') is %3.3f '%s)
Ans:
m : 2.5
n :3.5
gamma ( 3.5 ) is 3.323
Beta ( 2.5 3.5 ) is 0.037
4. Find the dimesion of subspace spanned by the vectors (1, 2, 3)
(2, 3, 1) and (3, 1, 2)
import numpy as np
import sympy as sp
# Define the vector space V
V = np.array([
[1, 2, 3],
[2, 3, 1],
[3, 1, 2]])
# Find the dimension and basis of V
basis = np.linalg.matrix_rank(V)
dimension = V.shape[0]
print("Basis of the matrix",basis)
print("Dimension of the matrix",dimension)
Ans:
Basis of the matrix 3
Dimension of the matrix 3
5. Verify the rank-nullity theorm for the linear transfromation
defined by . T(x, y, z) = (x + 4y + 7z, 2x + 5y +
8z, 3x + 6y + 9z)
import numpy as np
from scipy.linalg import null_space
# Define a linear transformation interms of matrix
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# Find the rank of the matrix A
rank = np.linalg.matrix_rank(A)
print("Rank of the matrix",rank)
# Find the null space of the matrix A
ns = null_space(A)
print("Null space of the matrix",ns)
# Find the dimension of the null space
nullity = ns.shape[1]
print("Null space of the matrix",nullity)
# Verify the rank-nullity theorem
if rank + nullity == A.shape[1]:
print("Rank-nullity theorem holds.")
else:
print("Rank-nullity theorem does not hold.")
Ans:
Rank of the matrix 2
Null space of the matrix [[-0.40824829]
[ 0.81649658]
[-0.40824829]]
Null space of the matrix 1
Rank-nullity theorem holds.
PART 4: Numerical methods
1. Obtain a root of the equation between 2 and 3 by
regula-falsi method.
from sympy import *
x= Symbol ('x') g = input ('Enter the function ')
f= lambdify (x,g)
a= float ( input ('Enter a valus :'))
b= float ( input ('Enter b valus :'))
N=int( input ('Enter number of iterations :'))
for i in range (1,N+1):
c=(a*f(b)-b*f(a))/(f(b)-f(a))
if ((f(a)*f(c)<0)):
b=c
else :
a=c
print ('itration %d \t the root %0.3f \t function value %0.3f \n'% (i,c,f(c)));
2. Find a root of the equation 3x = cos x+1, near 1, by Newton Raphson
method. Perform 5 iterations
from sympy import *
x= Symbol ('x')
g = input ('Enter the function ')
f= lambdify (x,g)
dg = diff (g);
df= lambdify (x,dg)
x0= float ( input ('Enter the initial approximation '));
n= int( input ('Enter the number of iterations '));
for i in range (1,n+1):
x1 =(x0 - (f(x0)/df(x0)))
print ('itration %d \t the root %0.3f \t function value %0.3f \n'% (i, x1
,f(x1))); x0 = x1
R K F method method:
1. Apply the Runge Kutta method to find the solution of at
y(0.1) taking h = 0.1. Given that y(0) = 1.
def f(x, y):
return x**2 + y**2
def runge_kutta(h, x0, y0, x):
n = int((x - x0) / h)
y = y0
for _ in range(n):
k1 = h * f(x0, y)
k2 = h * f(x0 + h / 2, y + k1 / 2)
k3 = h * f(x0 + h / 2, y + k2 / 2)
k4 = h * f(x0 + h, y + k3)
y += (k1 + 2 * k2 + 2 * k3 + k4) / 6
x0 += h
return y
x0 = 0
y0 = 1
x = 0.1
h = 0.1
result = runge_kutta(h, x0, y0, x)
print("y(0.1) =", result)
Ans: y(0.1) = 1.1114628561787105
2. Apply the Runge Kutta method to find the solution of at y(0.5)
taking h = 0.1. Given that y(0.4) = 1.
def f(x, y):
return 1/(x+y)
def runge_kutta(h, x0, y0, x):
n = int((x - x0) / h)
y = y0
for _ in range(n):
k1 = h * f(x0, y)
k2 = h * f(x0 + h / 2, y + k1 / 2)
k3 = h * f(x0 + h / 2, y + k2 / 2)
k4 = h * f(x0 + h, y + k3)
y += (k1 + 2 * k2 + 2 * k3 + k4) / 6
x0 += h
return y
x0 = 0.4
y0 = 1
x = 0.5
h = 0.1
result = runge_kutta(h, x0, y0, x)
print("y(0.5) =", result)
Ans: y(0.5) = 1
Taylor Series mehod:
def taylor_series(x, y, n):
h = x[1] - x[0]
y_approx = [y[0]]
for i in range(1, n+1):
delta_y = h**i * (x[i-1]**2 + y_approx[i-1]**2) / factorial(i)
y_approx.append(y_approx[i-1] + delta_y)
return y_approx
def factorial(n):
result = 1
for i in range(1, n+1):
result *= i
return result
# Initial condition
x0 = 0
y0 = 1
# Step size and number of terms
h = 0.1
n=3
# Generate x values
x = [x0 + i*h for i in range(n+1)]
# Generate y values using Taylor series method
y_approx = taylor_series(x, [y0], n)
# Print the approximated value of y at x=0.1
print(f"Approximate value of y at x=0.1: {y_approx[-1]}")
Ans: Approximate value of y at x=0.1: 1.1063105762016667
