LAB MANUAL

MATHEMATICS I FOR CSE STREAM (BMATS101)

COMPLIED BY
DR. D. UMADEVI, PROF. SNEHA
Department of Engineering Mathematics,
HKBK College of Engineering, Bengaluru.
 Department of Engineering Mathematics, HKBKCE

Department of Engineering Mathematics

List of Faculties in the Department

1. Dr. C.S NAGABHUSHANA Professor & HOD

2. Prof. UMME SALMA Assistant Professor

3. Dr. D. UMADEVI Associate Professor

4. Prof. SHARMADA.U Assistant Professor

5. Prof. SNEHA.S Assistant Professor

6. Prof. LAKSHMI.S Assistant Professor

7. Prof. AZRA BEGUM Assistant Professor

8. Prof. ROOPASHREE Assistant Professor

9. Prof. ARIF ALI Assistant Professor

10. Prof. ISHRATH Assistant Professor

11. Prof. NARESH Assistant Professor

12. Prof. ARADHANA C.K. Assistant Professor

13. Prof. RASHMI Assistant Professor

14. Prof. JAGADEESH Assistant Professor

1
 Department of Engineering Mathematics, HKBKCE

Programs for Mathematics I CSE Stream Lab

Table of Contents

S.No. Title Page No.

1 2D-Plots of Cartesian and Polar Curves 3

2 Finding Angle Between Two Polar Curves, Curvature and 12
Radius of Curvature
3 Finding Partial Derivatives and Jacobian 16
4 Taylor Series Expansion and L’Hospital’s Rule 19
5 Solution of First Order Differential Equations and Plotting 23
the Solution Curve
6 Finding GCD Using Euclid’s Algorithm 29
7 Solve Linear Congruence of the Form ax ≡ b(modn) 32
8 Numerical Solution of System of Equations, Test for 35
Consistency and Graphical
9 Solution of Linear Equations by Gauss-Seidel Method 39
10 Compute Eigen Value and Corresponding Eigen Vectors. 43
Find the Dominant Eigen Value and Corresponding Eigen
Vector by Rayleigh Power Method.

2
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 1: 2D-Plots of Cartesian and Polar Curves

The above functions are from matplotlib.pyplot library

The above functions are from the numpy library

1.1. Plotting Of Cartesian Curves

1. Write a code for Plotting points (Scattered plot) (1,2), (2,7), (3,9), (4,1), (6,5), (7,10), (8,3).

from matplotlib.pyplot import *

x = [1, 2, 3, 4, 6, 7, 8]
y = [2, 7, 9, 1, 5, 10, 3]
scatter(x, y) # plotting the points
xlabel('x - axis') # naming the x axis
ylabel ('y - axis') # naming the y axis
title ('Scatter points ') # giving a title to my graph
show()

3
 Department of Engineering Mathematics, HKBKCE

Output:

2. Write a code to plot a line (Line plot) passing through the points (1, 2), (2, 7), (3, 9), (4, 1),
(6, 5), (7, 10), (8, 3).

from matplotlib.pyplot import *

x = [1, 2, 3, 4, 6, 7, 8]
y = [2, 7, 9, 1, 5, 10, 3]
plot(x , y , 'r+--') # plotting the points
xlabel('x - axis ') # naming the x axis
ylabel('y - axis ') # naming the y axis
title('My first graph !') # giving a title to my graph
show ()
Output:

3. Write a code for plotting Sine and Cosine curves.

from numpy import *
from matplotlib. pyplot import *
x= arange(-10, 10, 0.001)
y1=sin(x)
y2=cos(x)
plot(x, y1, x, y2)
title("sine curve and cosine curve")
xlabel("Values of x")
ylabel("Values of sin (x) and cos(x)")
grid()
show()

4
 Department of Engineering Mathematics, HKBKCE

Output:

4. Write a code for plotting Exponential curve

from numpy import *

from matplotlib.pyplot import *
x = arange(-10, 10, 0.001)
y = exp(x)
plot(x, y)
title("Exponential curve")
grid()
show()

Output:

5. Write a code for plotting linear, quadratic and cubic curves

from matplotlib.pyplot import *

from numpy import *
x = linspace(0, 2, 100)
plot(x, x, label='linear')
plot(x, x**2, label='quadratic')
plot(x, x**3, label='cubic')
xlabel('x label')
ylabel('y label')
title("Simple Plot")
legend()
show()

5
 Department of Engineering Mathematics, HKBKCE

Output:

1.2. Implicit Functions

The above functions are from the sympy library

1. Write a Code to plot the equation of the circle 𝒙𝟐 + 𝒚𝟐 = 𝟒

from sympy import *

x, y = symbols('x y')
p1 = plot_implicit(Eq(x**2+y**2,4), (x,-4,4), (y,-4,4), title = 'Circle:
$x^2+y^2=4$')

Output:

6
 Department of Engineering Mathematics, HKBKCE

2. Write a Code to plot the equation of the Strophoid 𝒚𝟐 (𝒂 − 𝒙) = 𝒙𝟐 (𝒂 + 𝒙) take a = 2

from sympy import *

x, y = symbols('x y')
p2 = plot_implicit(Eq((y**2)*(2-x), (x**2)*(2+x)), (x, -5, 5), (y, -5, 5),
title='Strophoid: $y^2(a-x)=x^2(a+x), a=2$')

Output:

3. Write a code to plot the equation of the Cissiod: 𝒚𝟐 (𝒂 − 𝒙) = 𝒙𝟑 take a = 3

from sympy import *

x, y = symbols('x y')
p3 = plot_implicit(Eq((y**2)*(3-x), x**3), (x,-2,5), (y,-5,5), title
= ‘Cissiod: $y^2(a-x)=x^3$')

Output:

4. Write a code to plot the equation of Lemniscate: 𝒂𝟐 𝒚𝟐 = 𝒙𝟐 (𝒂𝟐 − 𝒙𝟐 ) take a = 2

from sympy import *

x, y = symbols ('x y')
p4 = plot_implicit(Eq(4*(y**2), (x**2)*(4-x**2)), (x,-5,5), (y,-5,5),
title ='Lemniscate: $a^2y^2=x^2(a^2-x^2)$')

7
 Department of Engineering Mathematics, HKBKCE

Output:

5. Write a code to plot the equation of Folium of De-Cartes: 𝒙𝟑 + 𝒚𝟑 = 𝟑𝒂𝒙𝒚, take a=2

from sympy import *

x, y = symbols('x y')
p5 = plot_implicit(Eq(x**3+y**3, 3*2*x*y), (x,-5,5), (y,-5,5), title=
'Folium of De-Cartes:$x^3+y^3=3axy$')

Output:

1.3. Polar Curves

PYLAB

• pylab is a historic interface. The equivalent replacement is matplotlib.pyplot.

• ‘from pylab import *’ imports all the functions from matplotlib.pyplot, numpy, numpy.fft,
numpy.linalg, and numpy.random, and some additional functions into the global namespace.

8
 Department of Engineering Mathematics, HKBKCE

1. Write a code to plot a curve of circle in polar form take 𝒓 = 𝟑

from pylab import *

axes(projection ='polar')
r = 3
rads = arange(0, (2*pi), 0.01)
for i in rads:
polar(i,r,'g.')
show()

Output:

2. Write a code to plot a Cardioid: r = 5(1 + cos θ)

from pylab import *

theta = linspace(0, 2*pi, 1000)
r1=5+5*cos(theta) c
polar(theta, r1,'r')
show()

Output:

3. Write a code to plot a four leaved Rose: 𝒓 = 𝟐|𝐜𝐨𝐬 𝟐𝒙|

from pylab import *

theta = linspace(0, 2*pi, 1000)
r = 2*abs(cos(2*theta))
polar(theta, r, 'r')
show()

9
 Department of Engineering Mathematics, HKBKCE

Output:

4. Write a code to plot a cardioids : 𝒓 = 𝒂 + 𝒂 𝐜𝐨𝐬 𝜽 and 𝒓 = 𝒂 − 𝒂 𝐜𝐨𝐬 𝜽, 𝒂 = 𝟑

from pylab import *

theta = linspace(0, 2*pi, 1000)
a = 3
r1 = a+a*cos(theta)
r2 = a-a*cos(theta)
polar(theta, r1,'r.', theta, r2,'g.')
show()

Output:

1.4. Parametric curves

1. Write a code to plot parametric equation of circle: 𝒙 = 𝒂 𝐜𝐨𝐬 𝜽 , 𝒚 = 𝒂 𝐬𝐢𝐧 𝜽 take 𝒂 = 𝟓
from pylab import *
theta = linspace(-2*pi, 2*pi, 1000)
a = 5
x = (a*cos(theta))
y = (a*sin(theta))
plot(x, y)
show()

10
 Department of Engineering Mathematics, HKBKCE

Output:

2. Write a code to plot parametric Equation of cycloid: 𝒙 = 𝒂(𝜽 − 𝐬𝐢𝐧 𝜽), 𝒚 = 𝒂(𝜽 − 𝒄𝒐𝒔𝜽)
take 𝒂 = 𝟐
from pylab import *
theta = linspace(-2*pi, 2*pi, 100)
a = 2
x = a*(theta-sin(theta))
y = a*(1-cos(theta))
plot(x,y)
show()

Output:

EXERCISE:
Plot the following :
1. Parabola 𝑦 2 = 4𝑎𝑥
𝑥2 𝑦2
2. Hyperbola − 𝑏2 = 1
𝑎2

3. Lower half of the circle 𝑥 2 + 2𝑥 = 4 + 4𝑦 − 𝑦 2

𝜋𝑥
4. cos ( 2 )
𝜋
5. 1+ sin (𝑥 + 4 )

6. Spiral of Archimedes: 𝑟 = 𝑎 + 𝑏𝜃
7. Limacon: 𝑟 = 𝑎 + 𝑏 cos 𝜃

11
 Department of Engineering Mathematics, HKBKCE

LAB Experiment 2: Finding angle between two polar curves,

curvature and radius of curvature

All the above functions are from sympy library

2.1. Angle between two polar curves :

𝑑𝜃 𝑑𝜃
Angle between radius vector and tangent is given by tan𝜙 = 𝑟 𝑑𝑟 ⇒ 𝜙 = tan−1 (𝑟 𝑑𝑟 )

Angle between two polar curves at the point of intersection is ∣ 𝜙1 − 𝜙2 ∣

1. Find the angle between the curves 𝒓 = 𝟒(𝟏 + 𝒄𝒐𝒔𝒕) and 𝒓 = 𝟓(𝟏 − 𝒄𝒐𝒔𝒕).

from sympy import*

r,t=symbols('r,t')
r1=4+4*cos(t)
r2=5-5*cos(t)
dr1=diff(r1,t)
dr2=diff(r2,t)
t1=r1/dr1
t2=r2/dr2
q=solve(r1-r2,t)
w1=t1.subs({t:float(q[0])})
w2=t2.subs({t:float(q[0])})
y1=atan(w1)
y2=atan(w2)
w=abs(y1-y2)
print('Angle between curves in radians is %0.3f'%(w))

Output: Angle between curves in radians is 1.571

12
 Department of Engineering Mathematics, HKBKCE

2. Finding the angle between the curves 𝒓 = 𝟒 𝒄𝒐𝒔𝒕 and 𝒓 = 𝟓 𝒔𝒊𝒏𝒕.

from sympy import*

r,t=symbols('r,t')
r1=4*cos(t)
r2=5*sin(t)
dr1=diff(r1,t)
dr2=diff(r2,t)
t1=r1/dr1
t2=r2/dr2
q=solve(r1-r2,t)
w1=t1.subs({t:float(q[0])})
w2=t2.subs({t:float(q[0])})
y1=atan(w1)
y2=atan(w2)
w=abs(y1-y2)
print('Angle between curves in radians is %0.3f'%(w))

Output: Angle between curves in radians is 1.571

2.2. Radius of curvature

3
(1+𝑦12 )2
Radius of curvature (Cartesian form), 𝜌 =
𝑦2
3
(𝑟 2 +𝑟12 )2
Radius of curvature (polar form), 𝜌 =
𝑟 2 +2𝑟12 −𝑟𝑟2

𝝅
1. Find the radius of curvature, 𝒓 = 𝟒(𝟏 + 𝒄𝒐𝒔𝒕) at 𝒕 = 𝟐 .

from sympy import*

r,t=symbols('r,t')
r=4*(1+cos(t))
r1=Derivative(r,t).doit()
r2=Derivative(r1,t).doit()
rho=(r**2+r1**2)**(1.5)/(r**2+2*r1**2-r*r2);
rho1=rho.subs(t,pi/2)
print('The radius of curvature is',(rho1))

Output: The radius of curvature is 3.77123616632825

13
 Department of Engineering Mathematics, HKBKCE

𝝅
2. Find the radius of curvature for 𝒓 = 𝒂𝒔𝒊𝒏(𝒏𝒕) at 𝒕 = and 𝒏 = 𝟏.
𝟐

from sympy import*

r,t,a,n=symbols('r,t,a,n')
r=a*sin(n*t)
r1=Derivative(r,t).doit()
r2=Derivative(r1,t).doit()
rho=(r**2+r1**2)**(1.5)/(r**2+2*r1**2-r*r2)
rho1=rho.subs([(t,pi/2),(n,1)])
print('The radius of curvature is')
display(simplify(rho1))

(𝑎2 )1.5
Output: The radius of curvature is
2𝑎2

2.3. Parametric curves

3
(1+𝑦12 )2
Radius of curvature (Cartesian form), 𝜌=
𝑦2

𝑑𝑦/𝑑𝑡 𝑑𝑦1 /𝑑𝑡

where 𝑦1 = and 𝑦2 =
𝑑𝑥/𝑑𝑡 𝑑𝑥/𝑑𝑡

1. Find radius of curvature and curvature of 𝒙 = 𝒂 𝐜𝐨𝐬(𝒕), 𝒚 = 𝒂 𝐬𝐢𝐧(𝒕).

from sympy import*

x,a,t,y=symbols('x,a,t,y')
y=a*sin(t)
x=a*cos(t)
y1=simplify(Derivative(y,t).doit())/simplify(Derivative(x,t).doit())
y2=simplify(Derivative(y1,t).doit())/simplify(Derivative(x,t).doit())
rho=simplify(1+y1**2)**(1.5)/y2
display('Radius of curvature is', ratsimp(rho))
rho1= rho.subs([(t,pi/2),(a,5)])
print('The radius of curvature at a=5, t=pi/2 is', rho1)
curvature=1/rho1
print('Curvature at (5,pi/2) is', float(curvature))

Output:

'Radius of curvature is'

1 1.5
−𝑎 ( ) sin (𝑡)3
sin(𝑡)2

The radius of curvature at a=5, t=pi/2 is -5

Curvature at (5, pi/2) is -0.2

14
 Department of Engineering Mathematics, HKBKCE

𝟑 𝟑
2. Find radius of curvature and curvature of 𝒙 = (𝒂 𝐜𝐨𝐬(𝒕)) ; 𝒚 = (𝒂 𝐬𝐢𝐧(𝒕))
𝟐 𝟐

from sympy import*

x,a,t,y=symbols('x,a,t,y')
x=(a*cos(t))**(3/2)
y=(a*sin(t))**(3/2)
y1=simplify(Derivative(y,t).doit())/simplify(Derivative(x,t).doit())
y2=simplify(Derivative(y1,t).doit())/simplify(Derivative(x,t).doit())
rho=simplify(1+y1**2)**(1.5)/y2
display('Radius of curvature is',ratsimp(rho))
rho1=rho.subs([(t,pi/4),(a,1)])
print('the radius of curvature at a=1, t=pi/4 is %0.4f'%rho1)
curvature=1/rho1
print('curvature at (1, pi/4)is %0.3f'%float(curvature))

Output:'Radius of curvature is'

1.5
𝑎 sin(𝑡)3.0
−3.0(a cos(𝑡))3.0 ( + 1) sin(𝑡)3 tan(𝑡)
𝑎 cos(𝑡)3 tan(𝑡)4
(𝑎 sin(𝑡))1.5 cos (𝑡)

the radius of curvature at a=1, t=pi/4 is -2.5227

curvature at (1, pi/4) is -0.396

EXERCISE

1. Find the angle between radius vector and tangent to the following polar curves
𝑎
a) 𝑟 = 𝑎𝜃 and 𝑟 = 𝜃
Ans: Angle between curves in radians is 90.000
b) 𝑟 = 2sin(𝜃) and 𝑟 = 2cos(𝜃)
Ans: Angle between curves in radians is 90.000
𝜋
2. Find the radius of curvature of 𝑟 = 𝑎(1 − cos(𝑡)) at 𝑡 = 2
1.5
0.942809041582063(𝑎2 )
Ans: 𝑎2
3. Find radius of curvature of 𝑥 = 𝑎cos3 (𝑡), 𝑦 = 𝑎sin3 (𝑡) at 𝑡 = 0.
Ans: 𝜌 = 0.75√3 and 𝜅 = 0.769800
𝜋
4. Find the radius of curvature of 𝑟 = 𝑎cos(𝑡) at 𝑡 = 4
1.5
(𝑎2 )
Ans: 2𝑎2
𝜋
5. Find the radius of curvature of 𝑥 = 𝑎(𝑡 − sin(𝑡)) and 𝑦 = 𝑎(1 − cos(𝑡)) at 𝑡 = 2 .
Ans: 𝜌 = 2.82842712 and 𝜅 = 0.353553

15
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 3: Finding partial derivatives and Jacobian

functions of several variables.

All the above functions are from sympy library.

sympy.abc - module exports all latin and greek letters as Symbols.

3.1. Partial derivatives

1. Prove that mixed partial derivatives, 𝒖𝒙𝒚 = 𝒖𝒚𝒙 for 𝒖 = 𝒆𝒙 (𝒙𝐜𝐨𝐬(𝒚) − 𝒚𝐬𝐢𝐧(𝒚)).

from sympy import*

x,y=symbols('x,y')
u=exp(x)*(x*cos(y)-y*sin(y))
uxy=diff(u,x,y)
uyx=diff(u,y,x)
if uxy==uyx:
print('Mixed partial derivatives are equal')
else:
print('Mixed partial deivatives are not equal')

Output: Mixed partial derivatives are equal

2. Prove that if 𝒖 = 𝒆𝒙 (𝒙𝐜𝐨𝐬(𝒚) − 𝒚𝐬𝐢𝐧(𝒚)), then 𝒖𝒙𝒙 + 𝒖𝒚𝒚 = 𝟎.

from sympy import*

x,y=symbols('x,y')
u=exp(x)*(x*cos(y)-y*sin(y))
uxx=diff(u,x,x)
uyy=diff(u,y,y)
w=simplify(uxx+uyy)
print('Answer=',w)

Output: Answer= 0

16
 Department of Engineering Mathematics, HKBKCE

3.2. Jacobians
𝒙𝒚 𝒚𝒛 𝒛𝒙
1. If 𝒖 = ,𝒗= ,𝒘= then prove that 𝑱 = 𝟒.
𝒛 𝒙 𝒚

from sympy import*

x,y,z=symbols('x,y,z')
u=x*y/z
v=y*z/x
w=z*x/y
J=Matrix([[diff(u,x),diff(u,y),diff(u,z)],[diff(v,x),diff(v,y),diff(v,z)],
[diff(w,x),diff(w,y),diff(w,z)]])
display ('The Jacobian matrix is',J)
J1=det(J).doit()
print("The Jacobian value is",J1)

Output: 'The Jacobian matrix is'

𝑦 𝑥 𝑥𝑦
− 2
𝑧 𝑧 𝑧
𝑦𝑧 𝑧 𝑦
− 2
𝑥 𝑥 𝑥
𝑧 𝑥𝑧 𝑥
− 2
[ 𝑦 𝑦 𝑦 ]
The Jacobian value is 4

2. If 𝒖 = 𝒙 + 𝟑𝒚𝟐 − 𝒛𝟑 , 𝒗 = 𝟒𝒙𝟐 𝒚𝒛, 𝒘 = 𝟐𝒛𝟐 − 𝒙𝒚, then prove that at (𝟏, −𝟏, 𝟎), 𝑱 = 𝟐𝟎.

from sympy import*

x,y,z=symbols('x,y,z')
u=x+3*y**2-z**3
v=4*x**2*y*z
w=2*z**2-x*y
J=Matrix([[diff(u,x),diff(u,y),diff(u,z)],[diff(v,x),diff(v,y),diff(v,z)],
[diff(w,x),diff(w,y),diff(w,z)]])
display('The Jacobian Matrix is',J)
J1=det(J).doit()
display(J1)
J2=J1.subs([(x,1),(y,-1),(z,0)])
print("The Jacobian value is",J2)

Output: 'The Jacobian Matrix is'

1 6𝑦 −3𝑧 2
[8𝑥𝑦𝑧 4𝑥 2 𝑧 4𝑥 2 𝑦 ]
−𝑦 −𝑥 4𝑧

4x3y - 24x2y3 + 12x2yz3 + 16x2z2 - 192xy2z2

The Jacobian value is 20
17
 Department of Engineering Mathematics, HKBKCE

𝛛(𝑿,𝒀,𝒁)
3. If 𝑿 = 𝝆𝐜𝐨𝐬(𝝓)𝐬𝐢𝐧(𝜽), 𝒀 = 𝝆𝐜𝐨𝐬(𝝓)𝐜𝐨𝐬(𝜽), 𝒁 = 𝝆𝐬𝐢𝐧(𝜽), then find .
𝛛(𝝆,𝝓,𝜽)

from sympy import *

from sympy.abc import *
X=rho*cos(phi)*sin(theta);
Y=rho*cos(phi)*cos(theta);
Z=rho*sin(phi);
J=Matrix([[diff(X,rho),diff(Y,rho),diff(Z,rho)],[diff(X,phi),diff(Y,phi),
diff(Z,phi)],[diff(X,theta),diff(Y,theta),diff(Z,theta)]])
print('The Jacobian matrix is')
display(J)
print('The Jacobian value is')
display(simplify(Determinant(J).doit()))

Output: The Jacobian matrix is

The Jacobian value is

𝜌2 cos (∅)

Exercise
𝑦 ∂2 𝑢 ∂2 𝑢
1. 𝑢 = 𝑡𝑎𝑛−1 (𝑥 ). Verify that ∂𝑦 ∂𝑥 = ∂𝑥 ∂𝑦

Ans: True
(𝑥 2 +𝑦 2 )
2. If 𝑢 = log (𝑥+𝑦)
, show that 𝑥𝑢𝑥 + 𝑦𝑢𝑦 = 1

Ans: True

3. If 𝑥 = 𝑢 − 𝑣, 𝑦 = 𝑣 − 𝑢𝑣𝑤 and 𝑧 = 𝑢𝑣𝑤, find Jacobian of 𝑥, 𝑦, 𝑧 w.r.t. 𝑢, 𝑣, 𝑤

Ans: 𝑢𝑣
∂(𝑥,𝑦)
4. If 𝑥 = 𝑟cos(𝑡) and 𝑦 = 𝑟sin(𝑡), then find ∂(𝑟,𝑡) .

Ans: 𝐽 = 𝑟
∂(𝑢,𝑣,𝑤)
5. If u=𝑥 + 3𝑦 2 − 𝑧 3 , v=4𝑥 2 𝑦𝑧 and w=2𝑧 2 − 𝑥𝑦 , find at (-2,-1,1) .
∂(𝑥,𝑦,𝑧)

18
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 4: Applications of Maxima and Minima of

functions of two variables, Taylor series expansion, and L' Hospital's
Rule

All the above functions are from sympy library.

4.1. Maxima and Minima problem

1. Find the Maxima and Minima of 𝒇(𝒙, 𝒚) = 𝒙𝟐 + 𝒙𝒚 + 𝒚𝟐 + 𝟑𝒙 − 𝟑𝒚 + 𝟒.

from numpy import *
from sympy import *
x,y=symbols('x,y')
f=x**2+x*y+y**2+3*x-3*y+4
fx=diff(f,x)
fy=diff(f,y)
p=solve([fx,fy],[x,y])
print("The stationary points are",p)
A=diff(fx,x)
B=diff(fy,x)
C=diff(fy,y)
A1=A.subs(p)
B1=B.subs(p)
C1=C.subs(p)
D=A1*C1*-B1**2
print("for",p,"D is",D,"and A is",A1)
if(D>0 and A1<0):
print("The Function attains Maxima")
Max=f.subs(p)
print("Maximum value is",Max)
elif(D>0 and A1>0):
print("The Function attains Minima")
Min=f.subs(p)
print("Minimum value is",Min)
elif(D<0):
print("It is a saddle point ")
elif(D==0):
print("Further test required")

19
 Department of Engineering Mathematics, HKBKCE

Output :

The stationary points are {x: -3, y: 3}

for {x: -3, y: 3} D is -4 and A is 2
It is a saddle point

4.2. Taylor Series and Maclaurin's Series Expansion

𝝅
1. Expand sin(𝒙) as Taylor series about 𝒙 = 𝟐 up to 3rd degree term. Also find 𝒔𝒊𝒏(𝟏𝟎𝟎∘ )

from sympy import *

x=Symbol('x')
y=sin(x)
y1=diff(y,x)
y2=diff(y,x,2)
y3=diff(y,x,3)
yx=lambdify(x,y)
y1x=lambdify(x,y1)
y2x=lambdify(x,y2)
y3x=lambdify(x,y3)
x0=float(pi/2)
TS=yx(x0)+(x-x0)*y1x(x0)+(x-x0)**2*y2x(x0)/2+(x-x0)**3*y3x(x0)/6
print("Taylor Series expansion is")
display(simplify(TS))
t=float(100*pi/180)
print("sin(100)=",yx(t))
Output:

Taylor Series expansion is

-1.02053899928946e-17 𝑥 3 - 0.5 𝑥 2 + 1.5707963267949 𝑥 - 0.23370055013617

sin(100)= 0.984807753012208

2. Find the Maclaurin’s series expansion of 𝒔𝒊𝒏(𝒙) + 𝒄𝒐𝒔(𝒙) upto 𝟑𝒓𝒅 degree term, calculate
𝒔𝒊𝒏(𝟏𝟎∘ ) + 𝒄𝒐𝒔(𝟏𝟎∘ ).
from sympy import *
x=Symbol('x')
y=sin(x)+cos(x)
y1=diff(y,x)
y2=diff(y,x,2)
y3=diff(y,x,3)
yx=lambdify(x,y)
y1x=lambdify(x,y1)
y2x=lambdify(x,y2)
y3x=lambdify(x,y3)
x0=float(pi/2)
MS=yx(0)+x*y1x(0)+x**2*y2x(0)/2+x**3*y3x(0)/6
print("Maclaurin Series expansion is")
display(simplify(MS))
m=float(10*pi/180)
print("sin(10)+cos(10)=",yx(m))

20
 Department of Engineering Mathematics, HKBKCE

Output:
Maclaurin Series expansion is

-0.166666666666667 𝑥 3 - 0.5 𝑥 2 + 1.0 𝑥 + 1.0

sin(10) + cos(10)= 1.1584559306791384

4.3. L’ Hospital Rule

sin 𝒙
1. Evaluate lim𝒙→𝟎
𝒙

from sympy import *

x=Symbol('x')
l=limit((sin(x))/x,x,0)
print("limit value =", l)

Output: limit value= 1

𝟓𝒙𝟒 −𝟒𝒙𝟐 −𝟏
2. Evaluate lim𝒙→𝟏
𝟏𝟎−𝒙−𝟗𝒙𝟑

from sympy import *

x=Symbol('x')
l=limit((5*x**4-4*x**2-1)/(10-x-9*x**3), x, 1)
print("limit value =", l)

Output: limit value= -3/7

𝟏 𝒙
3. Prove that lim𝒙→∞ (𝟏 + ) = 𝒆
𝒙

from sympy import *

x=Symbol('x')
l=limit((1+1/x)**x,x,float('inf'))
print("limit value =", l)

Output: limit value= E

21
 Department of Engineering Mathematics, HKBKCE

EXERCISE

1. Find the Taylor Series expansion of 𝑦 = 𝑒 −2𝑥 at 𝑥 = 0 upto third degree term.

Ans: −0.333333333333333𝑥 3 + 0.666666666666667𝑥 2 − 1.0𝑥 + 1.0

2
2. Expand 𝑦 = 𝑥𝑒 −3𝑥 as Maclaurin’s series upto fifth degree term.

Ans: 𝑥(0.75 ∗ 𝑥 4 − 0.75 ∗ 𝑥 2 + 0.5)

𝜋
3. Find the Taylor Series expansion of 𝑦 = 𝑐𝑜𝑠(𝑥) at 𝑥 = 3 .

Ans: 0.010464𝑥 4 + 0.00544𝑥 3 − 0.155467𝑥 2 − 0.1661389657𝑥 + 0.827151505

−1 (𝑥)
4. Find the Maclaurin’s series expansion of 𝑦 = 𝑒 −𝑠𝑖𝑛 at 𝑥 = 0 upto 𝑥 3 term.

Ans: −0.0833333333333333𝑥 3 + 0.166666666666667𝑥 2 − 0.5𝑥 + 1.0

2sin𝑥−sin2𝑥
5. Evaluate lim𝑥→0 𝑥–sin𝑥

Ans: 6

6. Evaluate lim𝑥→∞ √𝑥 2 + 𝑥 + 1 − √𝑥 2 + 1

Ans: 0.5

22
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 5: Solution of First Order differential

equation and plotting the solution curves.

All the above functions are from sympy library.

𝒅𝑷
1. Solve 𝒅𝒕
= 𝒓, 𝒓 = 𝟓, 𝑷(𝟎) = 𝟏 and plot the solution.
y

from sympy import*

import matplotlib.pyplot as plt
t,r=symbols('t,r')
P=Function('P')
de=Eq(Derivative(P(t),t),r)
display(de)
sol=dsolve(de,P(t),ics={P(0):1})
print("The solution of given linear differential equation is:")
display(sol)
sol=sol.subs(r,5)
print("The graph of solution curve ")
display(sol)
plot(sol.rhs)
plt.show()

Output:

23
 Department of Engineering Mathematics, HKBKCE

𝒅𝒚
2. Solve 𝒅𝒙 = −𝑲𝒚, 𝑲 = 𝟎. 𝟑, 𝒚(𝟎) = 𝟓 and plot the solution.

from sympy import*

import matplotlib.pyplot as plt
x,k=symbols('x, k')
y=Function('y')
de=Eq(Derivative(y(x),x),-k*y(x))
display(de)
sol=dsolve(de,y(x),ics={y(0):5})
print("The solution of given linear differential equation is:")
display(sol)
sol=sol.subs(k,0.3)
print("The graph of solution curve ")
display(sol)
plot(sol.rhs)
plt.show()

Output:

24
 Department of Engineering Mathematics, HKBKCE

𝒅𝒚
3. Solve 𝒙𝟑 𝒅𝒙 − 𝒙𝟐 𝒚 + 𝒚𝟒 𝒄𝒐𝒔𝒙 = 𝟎, 𝒚(𝝅) = 𝟏 and plot the solution.

from sympy import*

import matplotlib.pyplot as plt
x=symbols('x')
y=Function('y')
y1=Derivative(y(x),x)
de=Eq(x**3*y1-(x**2)*y(x)+(y(x)**4)*cos(x),0)
display(de)
sol=dsolve(de,y(x),ics={y(pi):1},hint="Bernoulli")
print("The solution and graph of given Bernoulli's differential equation is")
display(sol)
plot(sol.rhs)
plt.show()

Output:

25
 Department of Engineering Mathematics, HKBKCE

𝒅𝒚
4. Solve 𝒅𝒙 + 𝒚𝒕𝒂𝒏(𝒙) − 𝒚𝟑 𝒔𝒆𝒄(𝒙) = 𝟎, 𝒚(𝟎) = 𝟏 and plot the solution.

from sympy import*

import matplotlib.pyplot as plt
x=symbols('x')
y=Function('y')
y1=Derivative(y(x),x)
de=Eq(y1+y(x)*tan(x)-y(x)**3*sec(x),0)
display(de)
sol=dsolve(de,y(x),ics={y(0):1},hint="Bernoulli")
print("The solution and graph of given Bernoulli's differential equation is")
display(sol)
plot(sol.rhs)
plt.show()

Output:

26
 Department of Engineering Mathematics, HKBKCE

𝒅𝒚
5. Simulate 𝝉 𝒅𝒕 = −𝒚 + 𝑲𝒑; 𝑲𝒑 = 𝟑. 𝟎, 𝝉 = 𝟐. 𝟎 and 𝒚(𝟎) = 𝟏.

from sympy import*

import matplotlib.pyplot as plt
T,t,K=symbols('T, t, K')
y=Function('y')
de=Eq(T*Derivative(y(t),t),-y(t)+K)
display(de)
sol=dsolve(de,y(t),ics={y(0):1})
print("The solution of given linear differential equation is:")
display(sol)
sol=sol.subs([(T,2.0),(K,3.0)])
print("The graph of solution curve ")
display(sol)
plot(sol.rhs)
plt.show()

Output:

27
 Department of Engineering Mathematics, HKBKCE

EXERCISE
1. Solve the following differential equations and plot the solution curves:
a. y sin x dx − (1 + 𝑦 2 + cos2x) dy = 0.
1 3 𝑦3
Ans: 2 𝑦cos2𝑥 + (2) 𝑦 + = 0.
3
𝑑𝑦
b. = 𝑥 + 𝑦 subject to condition 𝑦(0) = 2.
𝑑𝑥

Ans: y = 3ex − x − 1.
𝑑𝑦
c. = 𝑥 2 subject to condition y(0) = 5.
𝑑𝑥

𝑥3
Ans: y = + 5.
3

d. 𝑥 2 𝑦 ′ = y log(y) − 𝑦 ′
−1 𝑥
Ans: y(x) = 𝑒 𝐶1 𝑡𝑎𝑛

e. 𝑦 ′ − 𝑦 − 𝑥𝑒 𝑥 = 0.
𝑥2
Ans: y(x) = (C1+ )𝑒 𝑥
2

28
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 6: Finding GCD using Euclid's algorithm

The above functions are from the sympy library.

Euclidean Algorithm:
1. For 𝑎 > 𝑏, 𝑎 = 𝑏 × 𝑞 + 𝑟 (𝑏𝑦 𝑑𝑖𝑣𝑖𝑠𝑖𝑜𝑛 𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚), where 0 ≤ 𝑟 < 𝑏
2. If r = 0, then GCD = b
3. If r ≠ 0, then assume 𝑎 = 𝑏 & 𝑏 = 𝑟 and repeat steps 1 to 3 until r = 0

6.1. Finding GCD of two numbers using Euclidean algorithm

1. Write a code to find the GCD of 614 and 124 by defining a new function using Euclidean
Algorithm

def gcdab(a,b):
r=1
if b<a:
a,b=b,a
while(r>0):
r=b%a
print(b,"=",a ,"x", b//a,"+",r)
b=a
a=r
continue
print('GCD =',b)
gcdab(614,124)

Output:

614 = 124 x 4 + 118

124 = 118 x 1 + 6
118 = 6 x 19 + 4
6 = 4 x 1 + 2
4 = 2 x 2 + 0
GCD = 2

29
 Department of Engineering Mathematics, HKBKCE

6.2. Identifying Relatively Prime Numbers

1. Write a code to check whether 163 and 512 are relatively prime

def rp(a,b):
r=1
a1,b1=a,b
if b<a:
a,b=b,a
while(r>0):
r=b%a
b=a
a=r
continue
if b==1:
print(f' {a1} and {b1} are relatively prime')
else:
print(f' {a1} and {b1} are not relatively prime')
rp(163,512)

Output: 163 and 512 are relatively prime

6.3. Checking divisibility

1. Write a code to check the number 8 divides the number 128.

def div(a,b):
r=1
a1,b1=a,b
if b<a:
a,b=b,a
while(r>0):
r=b%a
b=a
a=r
continue
if b==a1:
print(f' {a1} divides {b1} ')
else:
print(f' {a1} doesnot divides {b1}')
div(8,128)

Output: 8 divides 128

30
 Department of Engineering Mathematics, HKBKCE

6.4. Express GCD of a, b as a linear combination of a and b.

1. Calculate GCD of 76, 13 and express GCD as 76x + 13y.

def glin(a, b):

if b == 0:
return a, 1, 0
gcd, x1, y1 = glin(b, a % b)
x = y1
y = x1 - (a // b) * y1
return gcd, x, y
a, b = 76, 13
gcd, x, y = glin (a, b)
print(f"GCD of {a} and {b} is: {gcd}")
print(f"Linear Combination: {gcd} = {a}*({x}) + {b}*({y})")

Output:

GCD of 76 and 13 is: 1

Linear Combination: 1 = 76*(6) + 13*(-35)

EXERCISE

1. Find the GCD of 234 and 672 using Euclidean algorithm.

Ans: 6

2. What is the largest number that divides both 1024 and 1536?

Ans: 512

3. Find the greatest common divisor of 6096 and 5060?

Ans: 4

4. Prove that 1235 and 2311 are relatively prime.

31
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 7: Solving Linear congruence of the form

𝒂𝒙 ≡ 𝒃(𝒎𝒐𝒅𝒎).
Procedure to solve 𝒂𝒙 ≡ 𝒃(𝒎𝒐𝒅 𝒎)
(i) Find gcd(𝑎, 𝑚) = 𝑑
(ii) If 𝑑 does not divide 𝑏, then the linear congruence has no solution
(iii) If 𝑑 divides 𝑏, then the linear congruence has 𝑑 solutions
(𝑚 × 𝑖+𝑏)
(iv) Find an integer 𝑖 from 0 to 𝑚 − 1 such that 𝑥0 = 𝑎
is an integer and 𝑥0 is the initial soln.
𝑚
(v) Other solutions are given by 𝑥 = 𝑥0 + 𝑑
× 𝑡, where 𝑡 = 0,1,2, ⋯ , 𝑑 − 1

1. Show that the linear congruence 𝟔𝒙 ≡ 𝟓(𝒎𝒐𝒅𝟏𝟓) has no solution

from sympy import*
a=int(input('enter integer a '))
b=int(input('enter integer b '))
m=int(input('enter integer m '))
d=gcd(a,m)
if(b%d!=0):
print('The congruence has no integer solution')
else:
for i in range(0,m):
x0=(m*i+b)/a
if(x0//1==x0):
print(f'The {d} solutions are')
for j in range (0,d):
x=int(x0)+(m/d)*j
print(f' x = {x}(mod {m})')
break

Output:
enter integer a 6
enter integer b 5
enter integer m 15
the congruence has no integer solution

2. Find the solution of the congruence 𝟓𝒙 ≡ 𝟑(𝒎𝒐𝒅𝟏𝟑).

from sympy import*
a=int(input('enter integer a '))
b=int(input('enter integer b '))
m=int(input('enter integer m '))
d=gcd(a,m)
if(b%d!=0):
print('The congruence has no integer solution')
else:
for i in range(0,m):
x0=(m*i+b)/a
if(x0//1==x0):
print(f'The {d} solutions are')
for j in range (0,d):
x=int(x0)+(m/d)*j
print(f' x = {x}(mod {m})')
break

32
 Department of Engineering Mathematics, HKBKCE

Output:
enter integer a 5
enter integer b 3
enter integer m 13
The 1 solutions are
x = 11(mod 13)

Finding inverse of 𝑎 mod 𝑚 is equivalent to find 𝑎𝑥 ≡ 1(𝑚𝑜𝑑 𝑚).

3. Find the inverse of 5 mod 13

from sympy import*
a=int(input('enter integer a '))
b=int(input('enter integer b '))
m=int(input('enter integer m '))
d=gcd(a,m)
if(b%d!=0):
print('The congruence has no integer solution')
else:
for i in range(0,m):
x0=(m*i+b)/a
if(x0//1==x0):
print(f'The {d} solutions are')
for j in range (0,d):
x=int(x0)+(m/d)*j
print(f' x = {x}(mod {m})')
break

Output:
enter integer a 5
enter integer b 1
enter integer m 13
The 1 solutions are
x = 8(mod 13)

4. Find the solution of the linear congruence 𝟐𝟖𝒙 ≡ 𝟓𝟔(𝒎𝒐𝒅 𝟒𝟗)

from sympy import*

a=int(input('enter integer a '))
b=int(input('enter integer b '))
m=int(input('enter integer m '))
d=gcd(a,m)
if(b%d!=0):
print('The congruence has no integer solution')
else:
for i in range(0,m):
x0=(m*i+b)/a
if(x0//1==x0):
print(f'The {d} solutions are')
for j in range (0,d):
x=int(x0)+(m/d)*j
print(f' x = {x}(mod {m})')
break

33
 Department of Engineering Mathematics, HKBKCE

Output:
enter integer a 28
enter integer b 56
enter integer m 49
The 7 solutions are
x = 2(mod 49)
x = 9(mod 49)
x = 16(mod 49)
x = 23(mod 49)
x = 30(mod 49)
x = 37(mod 49)
x = 44(mod 49)

EXERCISE

1. Find the solution of the congruence 12𝑥 ≡ 6(𝑚𝑜𝑑23).

Ans: 12

2. Find the multiplicative inverse of 3 𝑚𝑜𝑑 31.

Ans: 21

3. Prove that 12𝑥 ≡ 7(𝑚𝑜𝑑14) has no solution. Give a reason for the answer.

Ans: Because GCD (12,14) = 2 and 2 doesnot divide 7.

34
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 8: Numerical solution of a system of equations,

test for consistency and graphical representation of the solution

8.1. System of homogenous linear equations:

The linear system of equations of the form AX=0 is called the system of homogenous linear
equations.

The n-tuple (0, 0, …, 0) is a trivial solution of the system. The homogeneous system of m
equations AX =0 in n unknowns has a non-trivial solution if and only if the rank of the matrix
A is less than n. Further, if 𝜌(𝐴) = 𝑟 < 𝑛, then the system possesses (n-r) linearly
independent solutions.

Functions with syntax Description

sympy library numpy library

Matrix ([[ row 1], [row 2], …, [row n]]) matrix ([[ row 1], [row 2], …, [row n]]) creates a m x n matrix .

A.rank ( ) linalg.matrix_rank ( A ) gives the rank of matrix A

shape (A) shape (A) gives the dimension of matrix A

A.shape [0], A.shape [1] A.shape [0], A.shape [1] gives the number of rows,
columns in matrix A respectively

A.col_insert (A.shape [1], B) concatenate((A,B), axis=1) creates augmented matrix AB

x,y,z = symbols(‘x, y, z’) linalg.solve (A, B) solves the system of equations

and returns the solutions of x, y, z
solve_linear_system (AB, x, y, z)

1. Check whether the following system of homogenous linear equation has non-trivial
solution. 𝒙𝟏 + 2𝒙𝟐 − 𝒙𝟑 = 𝟎 , 𝟐𝒙𝟏 + 𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎 , 𝟑𝒙𝟏 + 3𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎

from sympy import *

A=Matrix([[1 ,2 ,-1],[2 ,1 , 4],[3 ,3 , 4]])
B=Matrix([0,0,0])
r=A.rank()
n=A.shape[1]
print(f"The rank of the coefficient matrix is {r}")
print(f"The number of unknowns are {n}")
if (r==n):
print("System has trivial solution")
else:
print("System has", n-r, "non-trivial linearly independent solution(s)")

OUTPUT:

The rank of the coefficient matrix is 3

The number of unknowns are 3
System has trivial solution

35
 Department of Engineering Mathematics, HKBKCE

2. Check whether the following system of homogenous linear equation has non-trivial
solution 𝒙𝟏 + 2𝒙𝟐 − 𝒙𝟑 = 𝟎 , 𝟐𝒙𝟏 + 𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎 , 𝒙𝟏 − 𝒙𝟐 + 𝟓𝒙𝟑 = 𝟎

from sympy import *

A=Matrix([[1,2,-1],[2,1,4],[1,-1,5]])
B=Matrix([0,0,0])
r=A.rank()
n=A.shape[1]
print(f"The rank of the coefficient matrix is {r}")
print(f"The number of unknowns are {n}")
if (r==n):
print("System has trivial solution")
else:
print("System has", n-r, "non-trivial linearly independent solution(s)")

OUTPUT:

The rank of the coefficient matrix is 2

The number of unknowns are 3
System has 1 non-trivial linearly independent solution(s)

8.2. System of Non-homogenous Linear Equations

The linear system of equations of the form AX = B is called system of non-homogenous
linear equations if not all elements in B are zeros.
The non-homogeneous system of m equations AX=B in n unknowns is
• Consistent (has a solution) if and only if, 𝜌(𝐴) = 𝜌([𝐴|𝐵])
• has unique solution if 𝜌(𝐴) = 𝜌([𝐴|𝐵]) = 𝑛
• has infinitely many solutions, 𝜌(𝐴) < 𝑛
• inconsistent 𝜌(𝐴) ≠ 𝜌([𝐴|𝐵])
Functions with syntax Description
figure ( ) create a new figure, or activate an existing figure.
fig.add_subplot(111, projection='3d') defines new axes as a typical subplot, with a 3d projection,
to alert Matplotlib that we're about to throw three-
dimensional data at it.
ax.plot_surface (x,y,z,alpha=0.5) # surface plot is a representation of three-dimensional
dataset, where X and Y are 2D array of points of x and y
while Z is 2D array of heights.
# alpha parameter is used to control the transparency of a
plot. It takes a value between 0 (completely transparent) and
1 (completely opaque).
ax.scatter (X, Y, Z, color='red') marks the point x,y,z with the specified color.

The above functions are from matplotlibrary

36
 Department of Engineering Mathematics, HKBKCE

1. Examine the consistency of the following system of equations and solve if consistent, also
plot the graphical solution. 𝒙𝟏 + 2𝒙𝟐 − 𝒙𝟑 = 𝟎 , 𝟐𝒙𝟏 + 𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎 , 𝟑𝒙𝟏 +𝟑𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎

from numpy import *

from matplotlib.pyplot import *
A = matrix([[1, 2, -1], [2, 1, 4], [3,3,4]])
B = matrix([[1],[2],[1]])
AB=concatenate((A, B), axis=1)
if(linalg.matrix_rank(A)==linalg.matrix_rank(AB)):
if(linalg.matrix_rank(A)==A.shape [1]):
print("The system has unique solution")
else:
print("The system has infinitely many solutions")
x0 = linalg.solve(A, B)
print(x0)
X,Y,Z= x0[0],x0[1],x0[2]
x_lim,y_lim = linspace(-10, 10, 5),linspace(-10, 10, 5)
x, y = meshgrid(x_lim, y_lim)
z1,z2,z3 =(x+2*y - 1),(2-2*x-y)/4,(1-3*x-3*y)/4
fig = figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,z1,alpha=0.5),ax.plot_surface(x,y,z2,alpha=0.5),
ax.plot_surface(x,y,z3,alpha=0.5)
ax.scatter(X, Y, Z, color='red')
ax.set_xlabel('X'),ax.set_ylabel('Y'),ax.set_zlabel('Z')
show()
else:
print("The system of equations is inconsistent")

OUTPUT:
The system has unique solution
[[ 7.]
[-4.]
[-2.]]

37
 Department of Engineering Mathematics, HKBKCE

2. Examine the consistency of the following system of equations and solve and plot the
solution if consistent. 𝒙𝟏 + 2𝒙𝟐 − 𝒙𝟑 = 𝟎 , 𝟐𝒙𝟏 + 𝒙𝟐 + 𝟓𝒙𝟑 = 𝟎 , 𝟑𝒙𝟏 + 3𝒙𝟐 + 𝟒𝒙𝟑 = 𝟎

from numpy import *

from matplotlib.pyplot import *
A = matrix([[1, 2, -1], [2, 1, 5], [3,3,4]])
B = matrix([[1],[2],[1]])
AB=concatenate((A, B), axis=1)
if(linalg.matrix_rank(A)==linalg.matrix_rank(AB)):
if(linalg.matrix_rank(A)==A.shape [1]):
print("The system has unique solution")
else:
print("The system has infinitely many solutions")
x0 = linalg.solve(A, B)
print(x0)
X,Y,Z= x0[0],x0[1],x0[2]
x_lim,y_lim = linspace(-10, 10, 5),linspace(-10, 10, 5)
x, y = meshgrid(x_lim, y_lim)
z1,z2,z3 =(x+2*y - 1),(2-2*x-y)/4,(1-3*x-3*y)/4
fig = figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x,y,z1,alpha=0.5),ax.plot_surface(x,y,z2,alpha=0.5),
ax.plot_surface(x, y, z3, alpha=0.5)
ax.scatter(X, Y, Z, color='red')
ax.set_xlabel('X'),ax.set_ylabel('Y'),ax.set_zlabel('Z')
show()
else:
print("The system of equations is inconsistent")

OUTPUT:
The system of equations is inconsistent

EXERCISE

1. Find the solution of the homogenous system of equations

𝑥 + 𝑦 + 𝑧 = 0, 2𝑥 + 𝑦 − 3𝑧 = 0, 4𝑥 − 2𝑦 − 𝑧 = 0
2. Find the solution of the non-homogenous system of equations

25𝑥 + 𝑦 + 𝑧 = 27, 2𝑥 + 10𝑦 − 3𝑧 = 9, 4𝑥 − 2𝑦 − 𝑧 = −10

3. Find the solution of the non-homogenous system of equations
𝑥 + 𝑦 + 𝑧 = 2, 2𝑥 + 2𝑦 − 2𝑧 = 4, 𝑥 − 2𝑦 − 𝑧 = 5

38
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 9: Solution of a System of Linear Equations

by Gauss-Seidel Method

As we already know the def keyword is used to define a normal function in Python. Similarly,
the lambda keyword is used to define an anonymous function in Python
Syntax : lambda arguments : expression
Lambda function can have any number of arguments but only one expression, which is evaluated
and returned. It is efficient whenever one wants to create a function that only contains expressions
in a single line of a statement.

Procedure for Gauss–Seidel Method:

Given a system of linear equations in three unknowns:
𝑎11 𝑥1 + 𝑎12 𝑥2 + 𝑎13 𝑥3 = 𝑏1
𝑎21 𝑥1 + 𝑎22 𝑥2 + 𝑎23 𝑥3 = 𝑏2
𝑎31 𝑥1 + 𝑎32 𝑥2 + 𝑎33 𝑥3 = 𝑏3
• Gauss-Seidel can be applied only if the system of equations is diagonally dominant. That is,
|𝑎11 | > |𝑎12 | + |𝑎13 | , |𝑎22 | > |𝑎21 | + |𝑎23 | , |𝑎33 | > |𝑎31 | + |𝑎32 |
• If not diagonally dominant, then rearrange the system to satisfy the above condition. Write
1
the system of equations as 𝑥1 = [𝑏 − 𝑎12 𝑥2 − 𝑎13 𝑥3 ],
𝑎11 1
1
𝑥2 = [𝑏 − 𝑎21 𝑥1 − 𝑎23 𝑥3 ]
𝑎22 2
1
𝑥3 = [𝑏 − 𝑎31 𝑥1 − 𝑎32 𝑥2 ]
𝑎33 3
• Start the iteration with 𝑥1 = 0, 𝑥2 = 0, 𝑥3 = 0 as the initial approximation values.
• Keep substituting the recent values of 𝑥1 , 𝑥2 , 𝑥3 in the above formula for 𝑥1 , 𝑥2 , 𝑥3 in every
iteration. This process continues until we get the solution to the desired degree of accuracy.

1. Solve the system using Gauss-Seidel method:

𝟐𝟎𝒙 + 𝒚 − 𝟐𝒛 = 𝟏𝟕, 𝟑𝒙 + 𝟐𝟎𝒚 − 𝒛 = −𝟏𝟖, 𝟐𝒙 − 𝟑𝒚 + 𝟐𝟎𝒛 = 𝟐𝟓

f1 = lambda x,y,z:(17-y+2*z)/20
f2 = lambda x,y,z:(-18-3*x+z)/20
f3 = lambda x,y,z:(25-2*x+3*y)/20
x0, y0, z0 = 0, 0, 0
e = float(input('Enter tolerable error : '))
print('\t Iteration \t x \t y \t z\n')
for i in range (0, 25):
x1 = f1(x0, y0, z0)
y1 = f2(x1, y0, z0)
z1 = f3(x1, y1, z0)
print('\t\t %d \t %0.4f \t %0.4f \t %0.4f\n' %(i, x1, y1, z1))
e1, e2, e3 = abs(x0-x1), abs(y0-y1), abs(z0-z1)
x0,y0,z0 = x1,y1,z1
if e1>e and e2>e and e3>e:

39
 Department of Engineering Mathematics, HKBKCE

continue
else:
break
print('\n Solution : x = %0.3f , y = %0.3f and z = %0.3f\n'% (x1, y1, z1))
OUTPUT:
Enter tolerable error: 0.001
Iteration x y z

0 0.8500 -1.0275 1.0109

1 1.0025 -0.9998 0.9998

2 1.0000 -1.0000 1.0000

Solution: x = 1.000, y = -1.000 and z = 1.000

2. Solve the system using Gauss-Seidel method:

𝒙 + 𝟐𝒚 − 𝒛 = 𝟑, 𝟑𝒙 − 𝒚 + 𝟐𝒛 = 𝟏, 𝟐𝒙 − 𝟐𝒚 + 𝟔𝒛 = 𝟐

f1 = lambda x, y, z: (1+y-2*z)/3
f2 = lambda x, y, z: (3-x+z)/2
f3 = lambda x, y, z: (2-2*x+2*y)/6
x0, y0, z0 = 0, 0, 0
e = float(input('Enter tolerable error : '))
print('\t Iteration \t x\t y\t z\n')
for i in range(0, 25):
x1= f1(x0, y0, z0)
y1= f2(x1, y0, z0)
z1= f3(x1, y1, z0)
print('\t\t %d \t %0.4f \t %0.4f \t %0.4f \n' %(i, x1, y1, z1))
e1,e2,e3 = abs(x0-x1),abs(y0-y1),abs(z0-z1)
x0,y0,z0 = x1,y1,z1
if e1>e and e2>e and e3>e:
continue
else :
break
print('\n Solution : x = %0.3f, y = %0.3f and z = %0.3f\n'%(x1, y1, z1))

OUTPUT:

Enter tolerable error: 0.0001

Iteration x y z

0 0.3333 1.3333 0.6667

1 0.3333 1.6667 0.7778

Solution: x = 0.333, y = 1.667 and z = 0.778

40
 Department of Engineering Mathematics, HKBKCE

3. Solve the system using Gauss-Seidel method:

𝟏𝟎𝒙 + 𝒚 + 𝒛 = 𝟏𝟐, 𝒙 + 𝟏𝟎𝒚 + 𝒛 = 𝟏𝟐, 𝒙 + 𝒚 + 𝟏𝟎𝒛 = 𝟏𝟐.

f1 = lambda x, y, z: (12-y-z)/10
f2 = lambda x, y, z: (12-x-z)/10
f3 = lambda x, y, z: (12-x-y)/10
x0, y0, z0 = 0, 0, 0
e = float(input('Enter tolerable error: '))
print('\t Iteration \t x \t y \t z \n')
for i in range (0, 25):
x1 = f1(x0, y0, z0)
y1 = f2(x1, y0, z0)
z1 = f3(x1, y1, z0)
print('\t\t %d \t %0.4f \t %0.4f \t %0.4f\n' %(i, x1, y1, z1))
e1,e2,e3 = abs(x0-x1), abs(y0-y1), abs(z0-z1)
x0,y0,z0 = x1,y1,z1
if e1>e and e2>e and e3>e:
continue
else :
break
print ('\n Solution: x = %0.3f, y = %0.3f and z = %0.3f\n'% (x1, y1, z1))

OUTPUT:

Enter tolerable error: 0.0001

Iteration x y z

0 1.2000 1.0800 0.9720

1 0.9948 1.0033 1.0002

2 0.9996 1.0000 1.0000

3 1.0000 1.0000 1.0000

Solution: x = 1.000, y = 1.000 and z = 1.000

4. Solve the system using Gauss-Seidel method:

𝟓𝒙 − 𝒚 − 𝒛 = −𝟑, 𝒙 − 𝟓𝒚 + 𝒛 = −𝟗, 𝟐𝒙 + 𝒚 − 𝟒𝒛 = −𝟏𝟓

f1 = lambda x, y, z: (-3+y+z)/5
f2 = lambda x, y, z: (9+x+z)/5
f3 = lambda x, y, z: (15+2*x+y)/4
x0, y0, z0 = 0, 0, 0
e = float(input('Enter tolerable error: '))
print('\t Iteration \t x \t y \t z \n')
for i in range (0, 25):
x1 = f1(x0, y0, z0)

41
 Department of Engineering Mathematics, HKBKCE

y1 = f2(x1, y0, z0)

z1 = f3(x1, y1, z0)
print('\t\t %d \t %0.4f \t %0.4f \t %0.4f\n' %(i, x1, y1, z1))
e1,e2,e3 = abs(x0-x1),abs(y0-y1),abs(z0-z1)
x0,y0,z0 = x1,y1,z1
if e1>e and e2>e and e3>e:
continue
else :
break
print('\n Solution: x = %0.3f, y = %0.3f and z = %0.3f\n'%(x1, y1, z1))

OUTPUT:
Enter tolerable error: 0.0001
Iteration x y z

0 -0.6000 1.6800 3.8700

1 0.5100 2.6760 4.6740

2 0.8700 2.9088 4.9122

3 0.9642 2.9753 4.9759

4 0.9902 2.9932 4.9934

5 0.9973 2.9982 4.9982

6 0.9993 2.9995 4.9995

7 0.9998 2.9999 4.9999

8 0.9999 3.0000 5.0000

Solution: x = 1.000, y = 3.000 and z = 5.000

EXERCISE

1. Check whether the following system are diagonally dominant or not

(i) 25𝑥 + 𝑦 + 𝑧 = 27, 2𝑥 + 10𝑦 − 3𝑧 = 9, 4𝑥 − 2𝑦 − 12𝑧 = −10

(ii) 𝑥 + 𝑦 + 𝑧 = 7, 2𝑥 + 𝑦 − 3𝑧 = 3, 4𝑥 − 2𝑦 − 𝑧 = −1

2. Solve the following system of equations using Gauss Seidel method

(i) 4𝑥 + 𝑦 + 𝑧 = 6, 2𝑥 + 5𝑦 − 2𝑧 = 5, 𝑥 − 2𝑦 − 7𝑧 = −8
(ii) 27𝑥 + 6𝑦 − 𝑧 = 85, 6𝑥 + 15𝑦 + 2𝑧 = 72, 𝑥 + 𝑦 + 54𝑧 = 110

42
 Department of Engineering Mathematics, HKBKCE

LAB EXPERIMENT 10: Compute eigenvalues and corresponding

eigenvectors. Find dominant and corresponding eigenvector by
Rayleigh power method

Functions with syntax Description

linalg.eig( A ) computes the eigenvalues and eigenvectors of a square array/matrix.

dot(A, X) returns the dot product of vectors A and X.

The above functions are from numpy library.

10.1. Eigenvalues and Eigenvectors

Let A be a n × n matrix. λ is an eigenvalue of matrix A and x, a non-zero vector, is called an
eigenvector if it satisfies Ax = λx. We say, x is an eigenvector of A corresponding to
eigenvalue λ.
𝟒 𝟑 𝟐
1. Obtain the eigen values and eigen vectors for the given matrix (𝟏 𝟒 𝟏)
𝟑 𝟏𝟎 𝟒

from numpy import *

A = matrix([[4, 3, 2], [1, 4, 1], [3, 10, 4]])
w, v = linalg.eig(A)
print("\n Eigenvalues: \n", w)
print("\n Eigenvectors: \n", v)

OUTPUT:
Eigenvalues:
[8.98205672 2.12891771 0.88902557]

Eigenvectors:
[[-0.49247712 -0.82039552 -0.42973429]
[-0.26523242 0.14250681 -0.14817858]
[-0.82892584 0.55375355 0.89071407]]

𝟏 −𝟑 𝟑
2. Obtain the eigen values and eigen vectors for the given matrix (𝟑 −𝟓 𝟑)
𝟔 −𝟔 𝟒

from numpy import *

A = matrix([[1, -3, 3], [3, -5, 3], [6, -6, 4]])
w, v = linalg.eig(A)
print("\n Eigen values: \n", w)
print("\n Eigen vectors: \n", v)

OUTPUT:
Eigenvalues:
[ 4. -2. -2.]
43
 Department of Engineering Mathematics, HKBKCE

Eigenvectors:
[[ 0.40824829 -0.40824829 -0.30502542]
[ 0.40824829 0.40824829 -0.808424 ]
[ 0.81649658 0.81649658 -0.50339858]]

10.2. Largest eigenvalue and corresponding eigenvector by Rayleigh's power

method
Procedure for Rayleigh’s power method:
Given a matrix A,
• we assume the initial eigenvector X(0) as [1 0 0]T(or) [0 1 0]T (or) [0 0 1]T (or) [1 1
1]T and find the matrix product AX(0)
• take the largest absolute value outside from AX(0) as a common factor to obtain the
form AX(0) =𝜆(1)X(1) (This process is called normalization )
• again, find product AX(1) and normalize it to put in the form AX(1) =𝜆(2)X(2)
• this iterative process is continued till two consecutive iterative values of λ and X are
same upto a desired degree of accuracy.
• the values so obtained are respectively the largest eigenvalue λ and the
corresponding eigen vector X of the given matrix A.

𝟔 −𝟐 𝟐
1. Compute the numerically largest eigenvalue of P = [−𝟐 𝟑 −𝟏] by power method.
𝟐 −𝟏 𝟑

from numpy import *

def normalize(x):
lam = abs(x).max()
x_n = x/lam
return lam, x_n
x = matrix([[1], [1], [1]])
A = matrix([[6, -2, 2], [-2, 3, -1], [2, -1, 3]])
for i in range(10):
x1 = dot(A,x)
l, x = normalize(x1)
print('Eigenvalue:', l)
print('Eigenvector:', x)

OUTPUT:
Eigenvalue: 7.999988555930031
Eigenvector:
[[ 1. ]
[-0.49999785]
[ 0.50000072]]

44
 Department of Engineering Mathematics, HKBKCE

𝟏 𝟏 𝟑
2. Compute the numerically largest eigenvalue of P = [𝟏 𝟓 𝟏] by power method.
𝟑 𝟏 𝟏

from numpy import *

def normalize(x):
lam = abs(x).max()
x_n = x/lam
return lam, x_n
x = matrix([[1], [1], [1]])
A = matrix([[1, 1, 3], [1, 5, 1], [3, 1, 1]])
for i in range(10):
x1 = dot(A,x)
l, x = normalize(x1)
print('Eigenvalue:', l)
print('Eigenvector:', x)

OUTPUT:

Eigenvalue: 6.001465559355154
Eigenvector: [[0.5003663]
[1. ]
[0.5003663]]

EXERCISE
1. Find the eigenvalues and eigenvectors of the following matrices
25 1 2 11 1 2 3 1 1
25 1
a. P = [ ] b. P = [ 1 3 0] c. P = [ 0 10 0] d. P = [1 2 1]
1 3
2 0 −4 0 0 12 1 1 12
25 1 2
2. Find the dominant eigenvalue of the matrix P = [ 1 3 0 ]. Take 𝑋0 = (1, 0, 1)𝑇 .
2 0 −4
6 1 2
3. Find the dominant eigenvalue of the matrix P = [1 10 −1]. Take 𝑋0 = (1, 1, 1)𝑇 .
2 1 −4
5 1 1
4. Find the dominant eigenvalue of the matrix P = [1 3 −1] by power method.
2 −1 −4

45