Mathematics Laboratory Using Python

JSS MAHAVIDYAPEETHA
JSS ACADEMY OF TECHNICAL EDUCATION, BENGALURU
Affiliated to Visvesvaraya Technological University, Belagavi, Karnataka, INDIA
Approved by All India Council for Technical Education, New Delhi
UG programs accredited by NBA: ECE, CSE, ISE, CIVIL, E & IE and MECHANICAL

Mathematics Laboratory Using Python

Laboratory Manual

DEVELOPED BY
Dhanya P
Assistant Professor,
Department of Mathematics,
JSS Academy of Technical Education, Bangalore-60
 Mathematics Laboratory Using Python

List of Laboratory experiments (2 hours/week per batch)

10 lab sessions + 1 repetition class + 1 Lab Assessment

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, Jacobian and plotting the graph
4 Applications to Maxima and Minima of two variables
5 Solution of first-order differential equation and plotting the graphs
6
7
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 eigenvalues and eigenvectors and find the largest and smallest eigenvalue
by
Rayleigh power method.

6th and 7th programs are depend on the stream wise

CSE stream:

6. Finding GCD using Euclid’s Algorithm

7. Applications of Wilson's theorem

E&E stream:

6. Program to compute area, volume and centre of gravity

7. Evaluation of improper integrals

Mechanical and CIV Stream:

6. Solutions of Second-order ordinary differential equations with initial/boundary conditions

7. Solution of a differential equation of oscillations of a spring/deflection of a beam with

different loads
 Mathematics Laboratory Using Python

Prerequisites:
1. Basic mathematics.
2. Basic Programming skills.

Course Objective
The goal of this course is to teach the basics of Python programming in a manner targeted
at mathematics students:
1. To acquire programming skills in core Python.
2. To acquire Object Oriented Skills in Python
3. To develop the skill of designing Graphical user Interfaces in Python

Course outcome

Understand the basic syntaxes, operators and loops of python solve basic
CO1
Mathematical problems.
Make use of python: Identify, formulate and solve the Linear Differential
CO2
Equations.
Apply the python programming to find nth derivative and verify Euler’s
CO3
equation

CO4 Make use of Python programming to plot some standard curves.

CO5 Utilize the python programming, to find the area, volume and surface area
Solve ordinary differential equations of first order pertaining to engineering
CO6
problems by python
 Mathematics Laboratory Using Python

Introduction of Python Programming

Python is a object oriented and it is a widely used general purpose, high
level programming language created by Gudio Rossum in 1989. It is ideally designed for rapid
prototyping of complex applications. It was
mainly developed emphasis on code readability, and its syntax allows programmers to express
concepts in fewer lines of code.
Python is a programming language that lets you work quickly and integrate systems more
efficiently
Python does not allow punctuation characters such as @, $, &, % within identifiers.

Reserve words:
There are 33 reserve words are there in python and these reserve words, we cannot use
them as constant or variable or any other identifiers names. These keywords contain lowercase
letter only.

Value and Types:

A value is one of the basic building component in a program. Each value in python has a
type.
type( ) function is used in python shell to know which class the value belong to
ex: type(“JSSATE”) >>> class ‘str’
type(179) >>> class(int)
type(4.790) >>> class(float)

id( ) function:
id( ) is a built-in function in python, which returns the identity of an object. The identity
is a unique integer for that object during its life time.
Ex: a=2
Print(id(a)) >>>> 140454723286976

Statements and Expressions:

A statement is an instruction that the python interpreter can execute
It indicates some action to be carried out
Ex: while statement, for statement, if statement, import statement.
An expression is a combination of values, variables, operators and call of function.
 Mathematics Laboratory Using Python

Ex: x=1+3, y=x*2-5

Operators in python:
Operators are special symbol in python that carry out arithmetic or logical computation. The
value that the operator operates on is called operand.
1. Arithmetic operators:
The operator which consists the mathematical symbols line +, -, *, ÷, **, etc.,
+ → add 2 objects
- → subtract 2 objects
* → Multiplies 2 number or return the string repeated that many times
** → Power function
/ → Divide 2 number gives results in float
// → Divide 2 number and rounded by nearest whole number
% → Modulo returns the remainder of the division.

2. Relational operators:
The operator which gives relationship between any two object or values.
> → Grater symbol →true if left operand is greater than right
< → lesser symbol → true if left operand is less than the right
== → Equal to → true if the operands are equal
!= → Not equal to →true if the operands are not equal
>= → greater than or equal to → true if left operand is greater than or equal to the right
<= → less than or equal to → true if left operand is less than or equal to the right

3. Logical operator(Boolean operator):

There are only 3 basic logical operators: and, or, not
and → true if both the operators are true
or → true if either of the operator is true
not → true if operand is false

4. Assignment operator:
Assignment operator assign or updates values from right side operands to left side
operand.
If x is a variable then
= → x=5 ⇒x=5
+= →x+=5 ⇒ x=x+5
-= → x-=5 ⇒ x=x-5
*= → x*=5 ⇒ x=x*5
/= → x/=5 ⇒ x=x/5
%= → x%=5 ⇒ x=x%5
//= → x//=5 ⇒ x=x//5
**= → x**=5⇒ x=x**5

5. Special operators:
Special operators known as identity operator and Membership operators.
Identity operators are used to determine whether a value is of a certain class or type
 Mathematics Laboratory Using Python

is – true if the operands are identical

is not – true if the operands are not identical
Membership operators are used to validates the membership of a value. These operators are used
to find out whether a value is a member of a sequence such as string or list
in- true if value/variable is found in the sequence
not in- true if value/variable is not found in the sequence

Printing information:
In python, to print the message we used a function called print( )
Syntax: print(*object, sep=’,’end=’\n’)
Object- object to be printed
‘*’ Indicates that there may be more than one object
Sep- object are separated by sep and default value: ‘ ’
end- end is printed at last

Conditional execution:

• Simple if statement:
Syntax: if<< test expression>>:
Statement(S)
If the expression is true then only statement is executed

Ex: a=float(input(‘enter the value:’)

if a<9:
print(a, “is less than 9”)

• If else statement:
Syntax: if<< test expression>>:
body of it
else:
body of else

Ex: a=float(input(‘enter the value:’)

if a<9:
print(a, “is less than 9”)
else:
print(a, “is greater than 9”

• If, elif, else statement:

Syntax: if<< test expression 1>>
body of it
elif<<test expression 2>>:
body of elif

else:
 Mathematics Laboratory Using Python

body of else

Ex: a=float(input(‘enter the value:’)

if a<9:
print(a, “is less than 9”)
elif( a==9):
print(a, “is equal to 9”)
else:
print(a, “is greater than 9”

• For statement:
Syntax: for iterator_var in sequence:
Statement(s)
Ex: n = 4
for i in range(0, n):
print(i)

• while statement:
Syntax: while expression:
Statement(s)
Ex: count = 0
while (count < 3):
count = count + 1
print("Hello Geek")

• while else statement:

Syntax: while condition:
Statement(s)
else:
Statement(s)
Ex: count = 0
while (count < 3):
count = count + 1
print("Hello Geek")
else:
print("In Else Block")

SymPy:
SymPy is a Python library for symbolic mathematics. It aims to become a full-featured
computer algebra system (CAS) while keeping the code as simple as possible in order to be
comprehensible and easily extensible. SymPy is written entirely in Python.
By importing SymPy, simplify expressions, compute derivatives, integrals, and limits,
solve equations, work with matrices.
 Mathematics Laboratory Using Python

NumPy:
NumPy is a Python library used for working with arrays. It also has functions for
working in domain of linear algebra, fourier transform, and matrices. NumPy was created in
2005 by Travis Oliphant. It is an open source project and you can use it freely. NumPy stands
for Numerical Python.

def( ) keyword:
It is a keyword is used to define a function; it is placed before a function name that is
provided by the user to create a user-defined function.

Plotting the graph:

The plotting module allows you to make 2-dimensional and 3-dimensional plots. Presently the
plots are rendered using matplotlib as a backend. It is also possible to plot 2-dimensional plots
using a TextBackend if you don’t have matplotlib.

The plotting module has the following functions:

• plot: Plots 2D line plots.

• plot_parametric: Plots 2D parametric plots.
• plot_implicit: Plots 2D implicit and region plots.
• plot3d: Plots 3D plots of functions in two variables.
• plot3d_parametric_line: Plots 3D line plots, defined by a parameter.
• plot3d_parametric_surface: Plots 3D parametric surface plots.
• [Link]( ):To represent the title for the graph
• [Link]( ):To give the label for x-axis
• [Link]( ):To give the label for y-axis
• [Link]( ): To plot the scattered points
• [Link]( ): To get grid lines in the graph
• [Link]( ): To show the graph or To print the graph
• [Link]( ): To represent the index of the graph

The above functions are only for convenience and ease of use. It is possible to plot any plot by
passing the corresponding Series class to Plot as argument.

Plotting the Implicit function:

Syntax:
plot_implicit(expr,x_var = None, y_var = None, adaptive = True, depth = 0,
points = 300, line_cplor = ’blue’, show = True, **kwargs)

Plotting Polar curves:

 Mathematics Laboratory Using Python

Syntax:
matplotlib . pyplot . polar (theta , r, ** kwargs )

Eq( ): To represent the expression with in a bracket in the form of equation

diff( ):
Syntax: diff(function, independent variable, order of the differentiation).
This function is used to find the differentiation of the function w.r.t the independent variable at
any order.

Symbols ( ): Is the most important class in symPy library, symbolic computations are done with
symbols. SymPy variables are objects of Symbols class This function's argument is a string
containing symbol which can be assigned to a variable.

Matrix( ): Matrix([a11,112,a13],[a21,a22,a23],[a31,a32,a33]) This function is used to define a

matrix

.rank( ): To find the rank of a matrix.

.echelon_form( ) : To rduce the given matrix into row echelon form by Gauss elimination
method.
.rref( ): To reduce the given matrix into normal form bt Gauss Jardons method( also called as
row echelon form by Gauss Jardons Method)

.jacobian( ): This function helps to find the Jacobian’s of a function.

integrate( ): This function is used to integrate the function w.r.t independent variable.

Syntax: integrate(function,independent variable) {without limit}

integrate(function,[variable,lower limit,upper limit]) {with limit}

To write the comments for the program use “#”

 Mathematics Laboratory Using Python

Some simple programs to understand the Relational, Conditional and

Logical operators.
i) Compare two numbers (less than, greater than) using if statement

a=int(input("enter the value of A:")) Output: enter the value of A:8

b=int(input("enter the value of B:")) enter the value of
if (a > b): B:6
print('The number A is greater than B') The number A is
greater than B

a=int(input("enter the value of A:")) Output: 1. enter the value of A:8

b=int(input("enter the value of B:")) enter the value of B:6
if (a > b): The number A is greater than
print('The number A is greater than B') B
else: 2. enter the value of A:7
print('The number B is greater than A') enter the value of B:9
The number B is greater
than A

a=int(input("enter the value of A:")) Output: 1. enter the value of A:8

b=int(input("enter the value of B:")) enter the value of B:6
if (a > b): The number A is greater than
print('The number A is greater than B') B
elif a==b: 2. enter the value of A:7
print('The number A is equal to B') enter the value of B:9
else: The number B is greater than A
print('The number B is greater than A') 3. enter the value of A:7
enter the value of B:7
The number A is equal to B

ii) Sum of natural numbers using while loop

n=int(input(“enter the value of n:”)) Output: enter the value of n: 5
sum=0 The sum of 5 natural number is 15
i=1
while i<=n:
sum=sum+i
i+=1
print(‘The sum of’,n, ‘natural number is’, sum)
 Mathematics Laboratory Using Python

iii) Finding the factors of a number using for loop.

def factors(n): Output: factors(10)
counter=0
for x in range (1,n+1): Output: 1
if n%x==0: 2
print(x) 5
counter+=1 10
print('Number of factors=%d'%counter) Number of factors=4

iv) To check the given number is prime or not (use if... else statement).
def prime(n): Output: prime(10)
for i in range (2,n): The given number 10 is not a prime
if (n%i)==0: number
print("The given number ",n,"is not a prime Output: prime(11)
number") The given number 10 is a prime
break number
else:
print("The given number",n,"is a prime number")

v) Find the factorial of a number (use if...if...else).

n=int(input("Enter the value of n=")) Output:
factorial=1 Enter the value of n=5
if n<0: 5!=120
print("cannot compute factorial for negative
number")
elif (n<2):
print("{}!={}".format(n,factorial))
else:
for num in range (n,1,-1):
factorial=factorial*num
print("{}!={}".format(n,factorial))
 Mathematics Laboratory Using Python

vi) Simple programs to illustrate logical operators (and, or, not)

a = -10 Output:
b = 10 Either of the number is greater than 0
c=0
if a > 0 or b > 0: Atleast one number has boolean
print("Either of the number is greater than 0") value as False
else:
print("No number is greater than 0") 10 is divisible by either 3 or 5
a = 10
b = 12
c=0
if a and b and c:
print("All the numbers have boolean value as
True")
else:
print("Atleast one number has boolean value as
False")
a = 10
if not a:
print("Boolean value of a is True")
if not (a%3 == 0 or a%5 == 0):
print("10 is not divisible by either 3 or 5")
else:
print("10 is divisible by either 3 or 5")
 Mathematics Laboratory Using Python

1. 2D plots for Cartesian and polar curves

#importing the required module

import [Link] as plt

#plotting the scattered points

x=[1,2,3,4,5,6,7,8]
y=[2,7,9,11,13,14,16,17]
[Link](x,y)
[Link]('X-axis')
[Link]('Y-axis')
[Link]('Scatter points')
[Link]()

#importing the required module

import [Link] as plt
#plotting the lines
x=[1,2,3,4,5,6,7,8]
y=[2,7,9,11,13,14,16,17]

#plotting the curve with red color with

points in + and lines in -- for lines

[Link](x,y,'r+--')
[Link]('X-axis')
[Link]('Y-axis')
[Link]('Scatter points')
[Link]()
 Mathematics Laboratory Using Python

#importing the required module

From numpy import*
import [Link] as plt
x=arange(-10,10,0.001)
y=exp(x)
[Link](x,y)
[Link](“Exponential”)
[Link]()
[Link]()

from numpy import*

import [Link] as plt

x=arange(-10,10,0.001)
y1=sin(x)
y2=cos(x)

#plotting sine and cosine function

together with same value of x

[Link](x,y1,'r--',x,y2,'g')

#dotted line represent sine and green line

represent cosine

[Link]("sine and cosine curve")

[Link]("value of x")
[Link]("value of sin(x) and cos(x)")
[Link]()
[Link]()
 Mathematics Laboratory Using Python

# A simple graph
import matplotlib . pyplot as plt
import numpy as np

x = np . linspace (0 , 2 , 100 )
# Plot of y=x a linear curve
plt . plot (x , x , label ='linear ')
# Plot of y=x^2 a quadric curve
plt . plot (x , x ** 2 , label ='quadratic ')
# Plot of y=x^3 a cubic curve
plt . plot (x , x ** 3 , label ='cubic ')
plt . xlabel ('x label ')
plt . ylabel ('y label ')
plt . title (" Simple Plot ")
plt . legend ()
plt . show ()

Plotting Cartesian curve

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$ ') # r= 2
 Mathematics Laboratory Using Python

Strophoid:

#y**2*(a − x) = (x**2)*(a + x), a > 0

p3= 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> 0$ ') # a=2

Lemniscate:

#(a**2)*(y**2)=(a**2)*(a**2-x**2)

p5= plot_implicit (Eq ( 4*( y ** 2 ) ,( x ** 2 )*( 4-x ** 2 ) ) ,(x ,-5 , 5 ) ,(y ,-5 , 5 ) ) # a=2
 Mathematics Laboratory Using Python

Plotting Polar curve

Cardioid:

# Plot cardioid r=5(1+cos theta )

from pylab import *
theta = linspace (0,2*[Link] , 1000 )
r1=5+5*cos ( theta )
polar (theta ,r1 ,'r')
show ()

Four leaved Rose:

# Plot Four Leaved Rose r=2 | cos2x |

from pylab import *
theta = linspace (0,2*pi , 1000 )
r=2*abs(cos(2* theta ))
polar (theta ,r, 'r')
show ()
 Mathematics Laboratory Using Python

Plotting Parametric curve:

Cycloid:

def cycloid (r):

x = [] # create the list of x coordinates
y = [] # create the list of y coordinates
for theta in np. linspace (-2*[Link] , 2*[Link] , 100):
x. append (r*( theta - [Link]( theta )))
y. append (r*(1 - [Link]( theta ))) # same for y
plt. plot (x,y)
plt. show ()
cycloid (2)
 Mathematics Laboratory Using Python

Finding angle between two polar curves, curvature and radius of curvature.

1. Find the angle between the curves r = 4(1 + cos t) and r = 5(1 − cos t).

from sympy import * Output:

r,t = symbols ('r,t') Angle between curves in radians is 1.571
r1=3*(1+cos (t));
r2=6*(1-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[1])})
w2= t2. subs ({t: float (q[1])})
y1= atan (w1)
y2= atan (w2)
w=abs(y1-y2)
print ('Angle between curves in radians is %0.3f
'%(w))

2. Find the radius of curvature, r = 4(1 + cos t) at t=π/2.

from sympy import * Output:
t= Symbol ('t') The radius of curvature is 3.7712 units
r= Symbol ('r')
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) # substitute t in rho
print ('The radius of curvature is %3.4f units '%
rho1 )
 Mathematics Laboratory Using Python

3. Find radius of curvature of x = acos(t), y = asin(t).

from sympy import *
from sympy .abc import rho , x,y,r,K,t,a,b,c, alpha
y=( sqrt (x)-4) ** 2
y=a*sin(t)
x=a*cos(t)

dydx = simplify ( diff(y,t))/ simplify ( diff(x,t))

rho = simplify((1+dydx**2)**1.5/(diff(dydx ,t))/

(diff(x,t)))

print('Radius of curvature is ')

display ( ratsimp (rho ))

t1=pi/2
r1=5
rho1 = rho . subs (t,t1);
rho2 = rho1 . subs (a,r1);

print ('\n Radius of curvature at r=5 and t= pi/2 is ',

simplify ( rho2 ));

curvature =1/ rho2 ;

print ('\n Curvature at (5,pi/2) is ',float ( curvature
))

Finding partial derivatives and Jacobian of functions of several variables.

1. Prove that if 𝐮 = 𝐞𝐱 (𝐱 𝐜𝐨𝐬(𝐲) − 𝐲 𝐬𝐢𝐧(𝐲)) 𝐭𝐡𝐞𝐧 𝐮𝐱𝐱 + 𝐮𝐲𝐲 = 𝟎.

from sympy import *

x,y = symbols ('x y')
u=exp(x)*(x*cos(y)-y*sin(y))
display (u)
uxx = diff(u,x,x)
uyy = diff(u,y,y)
w=uxx+uyy
print('uxx = ',uxx)
print('uyy = ',uyy)
w1= simplify(w)
print ('Ans :',float (w1))
 Mathematics Laboratory Using Python

2. 𝐗 = 𝛒 ∗ 𝐜𝐨𝐬(𝛟) ∗ 𝐬𝐢𝐧(𝛉), 𝐘 = 𝛒 ∗ 𝐜𝐨𝐬(𝛟) ∗ 𝐜𝐨𝐬(𝛉), 𝐙 = 𝛒 ∗ 𝐬𝐢𝐧(𝛟) then

𝛛(𝐗,𝐘,𝐙)
Find 𝛛(𝛒,𝛟,𝛉).

#find the Jacobian of a function

from sympy import *
rho, phi, theta=symbols('rho, phi, theta')
x=rho*cos(phi)*sin(theta)
y=rho*cos(phi)*cos(theta)
z=rho*sin(phi)
A=Matrix([x,y,z])
jacA=[Link]([rho, phi, theta])
jacobian=simplify(det(jacA))
print('jacobian(A)=')
display(jacobian)

Applications of Maxima and Minima of functions of two variables, Taylor

series expansion and L’Hospital’s Rule
1. Find the Maxima and minima of 𝐟(𝐱, 𝐲) = 𝐱 𝟐 + 𝐲 𝟐 + 𝟑𝐱 − 𝟑𝐲 + 𝟒.
from sympy import*
x= Symbol ('x')
y= Symbol ('y')
f=x ** 2+x*y+y ** 2+3*x-3*y+4
d1= diff(f,x)
d2= diff(f,y)
criticalpoints1 = solve (d1)
criticalpoints2 = solve (d2)
print('criticalpoints1',criticalpoints1)
print('criticalpoints2',criticalpoints2)
s1= diff(f,x,2)
s2= diff(f,y,2)
s3= diff( Derivative (f,y),x)
print ('function value is ')
q1=s1. subs ({y: criticalpoints1 ,x: criticalpoints2 }). evalf ()
q2=s2. subs ({y: criticalpoints1 ,x: criticalpoints2 }). evalf ()
q3=s3. subs ({y: criticalpoints1 ,x: criticalpoints2 }). evalf ()
delta =s1*s2-s3 ** 2
print ('delta =',delta ,'\tq1 =', q1)
print ('delta =',delta ,'\tq2 =', q2)
print ('delta =',delta ,'\tq3 =', q3)
if( delta >0 and s1<0):
print (" f takes maximum value ")
elif ( delta >0 and s1>0):
print (" f takes minimum value ")
if ( delta <0):
print ("The point is a saddle point ")
if ( delta ==0):
print (" further tests required ")
 Mathematics Laboratory Using Python

2. Expand sin(x) as Taylor series about x = pi/2 upto 3rd degree term. Also find
sin(1000)

import numpy as np
from matplotlib import pyplot as plt
from sympy import *
x= Symbol ('x')
y=sin(x)
format
x0= float (pi/2)
dy= diff (y,x)
d2y = diff (y,x,2)
d3y = diff (y,x,3)
yat = lambdify (x,y)
dyat = lambdify (x,dy)
d2yat = lambdify (x,d2y)
d3yat = lambdify (x,d3y)

y=yat(x0)+((x-x0)/2)* dyat (x0)+((x-x0) ** 2/6)*

d2yat(x0)+((x-x0) ** 3/24)*d3yat (x0)

print ( simplify (y))

yat = lambdify (x,y)
print ("%.3f" % yat(pi/2+10*(pi/180)))
def f(x):
return [Link] (1*x)
x = np. linspace (-10 , 10)
[Link] (x, yat (x), 'r')
[Link] (x, f(x),'g')
[Link] ([-3, 3])
[Link] ()
[Link] ()

(𝟓𝐱 𝟒 − 𝟒𝐱 𝟐 − 𝟏)
𝟑. 𝐄𝐯𝐚𝐥𝐮𝐚𝐭𝐞 𝐥𝐢𝐦
𝒙→𝟏 𝟏𝟎 − 𝒙 − 𝟗𝒙𝟑

from sympy import * Output:

x= Symbol ('x') -3/7
l=limit((5*x ** 4-4*x ** 2-1)/(10-x-9*x ** 3),x,1)
print (l)
 Mathematics Laboratory Using Python

Solution of First order differential equation and plotting the solution curves
𝑑𝑦 1
1. Solve: = 1+𝑥 2 .
𝑑𝑥
from sympy import*
import numpy as np Output: solution is C1 + atan(x)
import matplotlib . pyplot as plt
from scipy . integrate import odeint
import math
def model (y,x):
dydx=1/(1+x**2)
return dydx
x=[Link](0,20)
y=odeint(model,1,x)
[Link](x,y)
plt . xlabel ('Time ')
plt . ylabel ('Response (y)')
plt . show ()
y=Function('y')
x=symbols('x')
eqn=Eq(y(x).diff(x),(1/(1+x**2)))
sol=simplify(dsolve(eqn,y(x),iex={y(0):1}))
print("solution is ",[Link])
𝑑𝑦
2. Solve 𝑑𝑥 − sin 𝑥 + 𝑘 = 0
from sympy import*
import numpy as np
import matplotlib . pyplot as plt
Output: solution is C1 - k*x - cos(x)
from scipy . integrate import odeint
import math
def model (y,x):
k=1
dydx=-k+sin(x)
return dydx
x=[Link](0,20)
y=odeint(model,0,x)
[Link](x,y)
plt . xlabel ('Time ')
plt . ylabel ('Response (y)')
plt . show ()
y=Function('y')
x,k=symbols('x k')
eqn=Eq(y(x).diff(x),(-k+sin(x)))
sol=simplify(dsolve(eqn,y(x),iex={y(0):0}))
print("solution is ",[Link])
 Mathematics Laboratory Using Python

3. The temperature of a body drops from 100 C to 75 C in 10 minutes where the surrounding
air is at the temperature 20 C . What will be the temperature of the body after half an
hour? Plot the graph of cooling.
import numpy as np Output: k= 0.0374693449441411
from sympy import * T= 20 + 80*exp(-k*t)
from matplotlib import pyplot as plt
t2=20 # surrounding temp
t1=100 # inital temp
# one reading t=1 minute temp is 75 degree
t=10
T=75
k1=(1/t)*log (( t1-t2)/(T-t2))# k calculation
print ('k= ',k1)
k= Symbol ('k')
t= Symbol ('t')
T= Function ('T')(t)
T=t2+(t1-t2)*exp(-k*t) # solution
print ('T=',T)
# ploting the solution curve
T=T. subs (k,k1)
T= lambdify (t,T)
t = np. linspace (0, 70)
plt . plot (t, T(t), 'r')
plt . grid ()
plt . show ()
# When time t=30 minute T is
print ('When time t=30 minute T is ,',T(30),'o
C')
 Mathematics Laboratory Using Python

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

representation of the solution.

1. Examine the consistency of the following system of equations and solve if consistent.
𝑥1 + 2𝑥2 − 𝑥3 = 1, 2𝑥1 + 𝑥2 + 4𝑥3 = 2, 3𝑥1 + 3𝑥2 + 4𝑥3 = 1.

from sympy import* Output:

x,y,z=symbols('x y z') Enter the number of unknowns = 3
n=int(input("Enter the number of unknowns = "))
A=Matrix([[1,2,-1],[2,1,4],[3,3,4]]) Rank of matrix A= 3
B=Matrix([[1],[2],[1]])
AB=Matrix([[1,2,-1],[2,1,4],[3,3,4],[1,2,1]]) Rank of Matrix AB= 3
Agu=AB
r=[Link]() The system as a unique solution
rAB=[Link]() {x: 7, y: -4, z: -2}
print("Rank of matrix A=",r,"\t Rank of Matrix AB=",rAB)
if (r!=rAB):
print("The system of equation is in-consistent")
elif (r==rAB==n):
print("The system as a unique solution")
print(solve([1*x+2*y-1*z-1,2*x+1*y+4*z-2,3*x+3*y+4*z-1],[x,y,z]))
else:
print("The system as a infinite solution")
print(solve([1*x+2*y-1*z-1,2*x+1*y+4*z-2,3*x+3*y+4*z-1],[x,y,z]))
 Mathematics Laboratory Using Python

2. Obtain the solution of 2x + y = 7; 3x − y = 3 graphically.

from sympy import *

import numpy as np
import matplotlib . pyplot as plt
x,y= symbols ('x,y')
sol = solve ([2*x+y-7,3*x-y-3],[x,y])
p=sol[x]
q=sol[y]
print ('Point of intersection is A (', p ,',', q, ')\n' )
x = np. arange (-10 , 10 , 0.001)
y1 = 7-2*x
y2=3*x-3
plt . plot (x,y1 ,'r')
plt . plot (x,y2 ,'g')
plt . plot (p,q, marker = 'o')
plt . annotate ('A', xy=(p,q), xytext =(p+0.5, q))
plt . xlim (-5,7)
plt . ylim (-7,7)
plt . axhline (y=0)
plt . axvline (x=0)
plt . title ("$2x+y=7; 3x -y=3$")
plt . xlabel (" Values of x")
plt . ylabel (" Values of y ")
plt . legend (['$2x+y=7$ ', '$3x -y-3$ '])
plt . grid ()
plt . show ()
 Mathematics Laboratory Using Python

Solution of system of linear equations by Gauss-Seidel method.

1. Solve the system of equations using Gauss-Seidel method:

20x + y − 2z = 17; 3x + 20y − z = −18; 2x − 3y + 20z = 25.

f1 = lambda x,y,z: (17-y+2*z)/20 Output:

f2 = lambda x,y,z: (-18-3*x+z)/20
f3 = lambda x,y,z: (25-2*x+3*y)/20 Enter tolerable error : 0.01
x0 = 0
y0 = 0 Count x y z
z0 = 0
count = 1 1 0.8500 -1.0275 1.0109
e = float ( input ('Enter tolerable error : '))
print ('\n Count \tx\ty\tz\n') 2 1.0025 -0.9998 0.9998
condition = True
while condition : 3 1.0000 -1.0000 1.0000
x1 = f1(x0 ,y0 ,z0)
y1 = f2(x1 ,y0 ,z0)
z1 = f3(x1 ,y1 ,z0) Solution : x=1.000 , y=-1.000 and z = 1.000
print ('%d\t%0.4f\t%0.4f\t%0.4f\n' %(count , x1 ,y1,z1))
e1 = abs (x0-x1);
e2 = abs (y0-y1);
e3 = abs (z0-z1);
count += 1
x0 = x1
y0 = y1
z0 = z1
condition = e1>e and e2>e and e3>e
print ('\n Solution : x=%0.3f , y=%0.3f and z = %0.3f\n'%
(x1 ,y1 ,z1))
 Mathematics Laboratory Using Python

Compute Eigen-values and corresponding eigenvectors. Find dominant and

corresponding eigenvector by Rayleigh power method.

4 3 2
1. Obtain the Eigen values and Eigen vectors for the given matrix.[1 4 1]
3 10 4
import numpy as np Output:Given matrix :
I=np. array ([[4,3,2],[1,4,1],[3,10 ,4]]) [[ 4 3 2]
print ("\n Given matrix : \n", I) [ 1 4 1]
#x=np. linalg . eigvals (I) [ 3 10 4]]
w,v = np. linalg .eig (I)
print ("\n Eigen values : \n", w) Eigen values :
print ("\n Eigen vectors : \n", v) [8.98205672 2.12891771 0.88902557]
## To display one eigen value and correspondingeigen
vector Eigen vectors :
print (" Eigen value :\n ", w[0]) [[-0.49247712 -0.82039552 -0.42973429]
print ("\n Corresponding Eigen vector :", v[:,0]) [-0.26523242 0.14250681 -0.14817858]
[-0.82892584 0.55375355 0.89071407]]
Eigen value :
8.982056720677651

Corresponding Eigen vector : [-0.49247712 -0.2

6523242 -0.82892584]

1 1 3
2. Compute the numerically largest Eigen-value of [1 5 1] by power method.
3 1 1

import numpy as np Output:

def normalize (x): Eigenvalue : 6.001465559355154
fac = abs(x).max () Eigenvector : [0.5003663 1. 0.5003663]
x_n = x / [Link] ()
return fac , x_n
x = np. array ([1, 1,1])
a = np. array ([[1,1,3 ],[1,5,1],[3,1,1]])
for i in range (10):
x = [Link](a, x)
lambda_1 , x = normalize (x)
print ('Eigenvalue :', lambda_1 )
print (' Eigenvector :', x)
 Mathematics Laboratory Using Python

Computer Science and Engineering Stream

6. Finding GCD using Euclid’s algorithm.

def gcd1 (a,b): Output:
c=1 gcd1 (620 ,123 )
if b < a:
t=b 123 5
b=a 53
a=t 32
while (c > 0): 21
c = b%a 10
print (a,c) GCD = 1
b=a
a=c
continue
print ('GCD =',b)

7. Solving linear congruence of the form ax ≡b( mod m).

from sympy import * Output:
from math import * enter integer a 6
a=int( input ('enter integer a ')); enter integer b 3
b=int( input ('enter integer b ')); enter integer m 3
m=int( input ('enter integer m ')); the solution of the congruence is 1.0
d=gcd(a,m)
if (b%d!=0):# Reminder calculation
print ('the congruence has no integer solution ');
else :
for i in range (1,m-1):
x=(m/a)*i+(b/a)
if(x // 1==x):# check whether x is an integer
print ('the solution of the congruence is ', x)
break
 Mathematics Laboratory Using Python

Electrical & Electronics Engineering Stream

6. a. Find the area of an ellipse by double integration
from sympy import * Output:
x= Symbol ('x')
y= Symbol ('y') 24.0*pi
#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)
b. Find the volume of the tetrahedron bounded by the planes
from sympy import * Output:
x= Symbol ('x') a*b*c/6
y= Symbol ('y')
z= Symbol ('z')
a= Symbol ('a')
b= Symbol ('b')
c= Symbol ('c')
w2= integrate (1 ,(z,0,c*(1-x/a-y/b)) ,(y,0,b*(1-x/a)) ,(x,0,a))
print (w2)
c. Find the center of gravity of cardioids. Plot the graph of cardioids and mark the center of
gravity.
import numpy as np Output:
import [Link] as plt 5*a/6
import math
from sympy import *
r= Symbol ('r')
t= Symbol ('t')
a= Symbol ('a')
I1= integrate (cos (t)*r ** 2 ,(r,0,a*(1+cos (t))) ,(t,-pi ,pi))
I2= integrate (r ,(r,0,a*(1+cos (t))) ,(t,-pi ,pi))
I=I1/I2
print (I)
I=[Link] (a,5)
[Link](projection ='polar')
a=5
rad = [Link] (0, (2 * [Link]), 0.01)
for i in rad :
r = a + (a*[Link](i))
plt. polar (i,r,'g.')
[Link](0,I,'r.')
[Link]()
 Mathematics Laboratory Using Python

Evaluation of improper integrals, Beta and Gamma functions

5 7 5
1. Calculate Beta(2 , 2) and Gamma(2).

from sympy import beta , gamma Output:

m= float ( input ('m : ')); m : 3.5
n= float ( input ('n :')); n :5.5
s= beta (m,n); gamma ( 5.5 ) is 52.343
t= gamma (n) Beta ( 3.5 5.5 ) is 0.004
print ('gamma (',n,') is %3.3f '%t)
print ('Beta (',m,n,') is %3.3f '%s)

2. Verify that Beta(m, n) = Gamma(m)Gamma(n)/Gamma(m + n) for m=5 and n=7

from sympy import beta , gamma Output:

m=5; 0.000432900432900433 0.000432900432900433
n=7; beta and gamma are related
m= float (m);
n= float (n);
s= beta (m,n);
t=( gamma (m)* gamma (n))/ gamma (m+n);
print (s,t)
if (abs (s-t)<=0.00001 ):
print ('beta and gamma are related ')
else :
print ('given values are wrong ')
 Mathematics Laboratory Using Python

Mechanical & Civil Engineering Stream

Solution of second order ordinary differential equation and plotting the solution
curve
1. Plot the solution curves of 𝑦 ′′ + 2𝑦 ′ + 2𝑦 = cos(2𝑥) , 𝑦(0) = 0, 𝑦 ′ (0) = 0
import numpy as np
from scipy . integrate import odeint
import matplotlib . pyplot as plt
def dU_dx (U, x):
return [U[1], -2*U[1] - 2*U[0] + [Link] (2*x)]
U0 = [0, 0]
xs = np. linspace (0, 10 , 200 )
Us = odeint (dU_dx , U0 , xs)
ys = Us[:,0]
ys1 =Us[:,1]
plt . xlabel ("x")
plt . ylabel ("y")
plt . title (" Solution curves ")
plt . plot (xs ,ys , label ='y');
plt . plot (xs ,ys1 , label ='z');
plt . legend ()
plt . show ()
 Mathematics Laboratory Using Python

Solution of differential equation of oscillations of a spring with various load

𝑑2𝑥 dx
1. Solve 9 𝑑𝑡 2 + 2 dt + 1.2 x = 0, x(0) = 1.5, x ′ (0) = 2.5 and plot the solution curve .
import numpy as np Output:
from scipy . integrate import odeint [[ 2.5 1.5 ]
import matplotlib . pyplot as plt [ 2.546773 1.47613782]
def f(u,x): [ 2.59279255 1.45205195]
return (u[1],-(1/9)*(1.2*u[1]+2*u[0])) ...
y0=[2.5,1.5] [-0.06627678 0.00620788]
xs=np. linspace (0,20*[Link] , 2000 ) [-0.06607481 0.00664319]
us= odeint (f,y0 ,xs) [-0.0658592 0.00707522]]
print (us)
ys=us[:,0]
plt . plot (xs ,ys ,'r-')
plt . xlabel ('Time ')
plt . ylabel ('Amplitude ')
plt . title ('Solution of free and damped case ')
plt . grid ( True )
plt . show ()