Institute of Chemical Technology

Department of Mathematics
Programming Lab I (MAP 2521)
End-Semester, Date: Dec 19, 2022

Timing: 1:30 PM to 4:30 PM Maximum marks: 50

Attempt All 5
1. Consider the sequence of functions defined on Ω = (0, 1) given by
 
k 1 k 1

j 2 /2 for x ∈ − , + where k = 1, 2, 3, · · · , j
fj (x) = j + 1 j3 j + 1 j3
0 otherwise


Plot fj (x) for 1 ≤ j ≤ 5 and label them using an appropriate legend using MatPlotLib
or ggplot2. (6)
2. Consider the metric space X = (R2 , | |). Construct a randomly generated matrix
M50×50 (X). For a fixed v ∈ X of your choice, using NumPy, write a program to find
argminw∈M |w − v|. (6)
3. Using SymPy,
(a) Verify the following integral
Z
cos x π
I= 2
dx =
R x +1 e
Pn
(b) Given Tn = k=1 k, verify the following sum
∞
X 1
=2
k=1
Tk

(6)
4. Define a class by the name ODESolver which is initalized with a matrix M ∈ GL2 (R),
x⃗0 ∈ R2 , and n ∈ R. Using SciPy, equip the class with the following methods
(a) solve(self) which plots the solution the system of ODEs defined by ⃗x′ = M⃗x
with initial condition x⃗0 over time 0 ≤ t ≤ n.
(b) phase_portrait(self,a,b,c,d) which plots the phase portrait of the above
given system when the initial conditions vary discretely over [a, b] × [c, d].
(6)
5. Two real floats x and y are chosen at random from the interval (0, 1). Compute the
probability that the closest integer to x/y is even. Using R, verify that the probability
is precisely (5 − π)/4. (6)
 Attempt Any 4 (Out Of The 5)

1. Let p1 = (0, 0), p2 = (4, 0), and p3 = (2, 3) be the corners of the equilateral triangle
∆. Fix a v0 ∈ R2 of your choice such that v0 lies in the interior of ∆. Using R or
Python, compute and plot the sequence of iterates governed by the rule
1
vn+1 = (vn + prn )
2
where rn is chosen at random from the set {1, 2, 3}. (5)

2. Consider the matrix A15×15 such that the ij-th entry of A is given by

aij = ij−1 (mod j)

Plot A using R or Python to verify Fermat’s little theorem. (5)

3. Let σ(n) denote the sum of the divisors of n. Let S(n, d) denote the sum of the
numbers i ≤ n such that d divides σ(i). Using R or Python verify that S(20, 7) = 49. (5)

4. Let ⃗a and ⃗b(̸= ⃗0) be vectors from R2 . Using R or Python write a function to plot ⃗a, ⃗b,
and the vector projection of ⃗a on ⃗b. (5)

5. Let f : (−1, 1) → R be such that f (x) = cos(x) − x3 . Fix an x0 (̸= 0) ∈ (−1, 1) of

your choice. Using R or Python compute the sequence (xn )n≥1 governed by the rule

f (xn−1 )
xn = xn−1 −
f ′ (xn−1 )

Is (xn )n≥1 convergent? If so, does it converge to the root of f (x)? (5)

Happy Coding!
19/12/2022, 19:55 Endsem Solutions

In [1]: import numpy as np

import matplotlib.pyplot as plt

import math

import sympy as sp

from scipy import integrate

Section A

Question 1

In [2]: def in_intervals(x,I):

def in_interval(x,I):

if I[0]<x<I[1]:

return(True)

else:

return(False)

bools = []

for i in I:

if in_interval(x,i):

bools.append(True)

else:

bools.append(False)

bools = list(set(bools))

if True in bools:

return(True)

else:

return(False)

In [3]: def generate_intervals(j):

intervals = []
for k in range(1,j+1):

intervals.append( (k/(j+1) - j**(-3), k/(j+1) + j**(-3)) )

return(intervals)

In [4]: def f(j,x):

I = generate_intervals(j)

if in_intervals(x,I):

return(j**2/2)

else:

return(0)

In [5]: xvals = np.linspace(0,1,1000)

fig = plt.figure(figsize=(6,6))

ax = plt.axes()

for j in range(1,6):

yvals = [f(j,x) for x in xvals]

ax.plot(xvals,yvals,label = 'j = '+str(j))

plt.legend()

file:///home/kanak/Documents/Endsem Solutions (1).html 1/7

19/12/2022, 19:55 Endsem Solutions

Question 2

In [6]: M = np.zeros((50,50,2))

X = np.random.uniform(0,10,(50,50))

Y = np.random.uniform(0,10,(50,50))

M[:,:,0] = X
M[:,:,1] = Y

In [7]: v = np.array([0,1])

dis = 1000

for i in range(50):

for j in range(50):

w = np.array([X[(i,j)],Y[(i,j)]])

if math.dist(v,w)<dis:

dis = math.dist(v,w)

can = w

array([0.95487155, 7.46257532])
Out[7]:

Question 3

In [8]: x = sp.symbols('x')

f = sp.cos(x)/(x**2 + 1)

sp.integrate(f,(x,-sp.oo,sp.oo)).simplify()

π
Out[8]:
e

In [9]: n, k = sp.symbols('n k')

def T(n):

file:///home/kanak/Documents/Endsem Solutions (1).html 2/7

19/12/2022, 19:55 Endsem Solutions

return(sp.summation(k,(k,1,n)))

sp.summation(1/T(k),(k,1,sp.oo))

Out[9]: 2

Question 4

In [10]: class ODESolver:

def __init__(self,M,x0,n):

self.M = M

self.x0 = x0

self.n = n

def solve(self):

def model(X,t):

return([ self.M[(0,0)]*X[0] + self.M[(0,1)]*X[1],

self.M[(1,0)]*X[0] + self.M[(1,1)]*X[1] ])

t = np.linspace(0,self.n,1000)

S = integrate.odeint(model,self.x0,t)

x,y = S.T

fig = plt.figure(figsize=(5,5))

ax = plt.axes()

ax.plot(t,x,label='x')

ax.plot(t,y,label='y')

plt.legend()

plt.show()

def phase_portrait(self,a,b,c,d):

def model(X,t):

return([ self.M[(0,0)]*X[0] + self.M[(0,1)]*X[1],

self.M[(1,0)]*X[0] + self.M[(1,1)]*X[1] ])

t = np.linspace(0,self.n,1000)

x = np.linspace(a,b,100)

y = np.linspace(c,d,100)

X,Y = np.meshgrid(x,y)

U,V = model([X,Y],t)

col = np.hypot(U,V)

fig = plt.figure(figsize=(5,5))

ax = plt.axes()

ax.streamplot(X,Y,U,V,density=1,cmap='Reds_r',color=col)

In [11]: Q = ODESolver(np.array([[1,2],[3,2]]),(0.5,0.5),2)

In [12]: Q.solve()

file:///home/kanak/Documents/Endsem Solutions (1).html 3/7

19/12/2022, 19:55 Endsem Solutions

In [13]: Q.phase_portrait(-2,2,-2,2)

Question 5

In [14]: (5-np.pi)/4

0.4646018366025517
Out[14]:

In [15]: counts = []

for i in range(10000):

x = np.random.uniform()

y = np.random.uniform()

if (np.round(x/y))%2 == 0:

counts.append(1)

len(counts)/10000

0.4591
Out[15]:

file:///home/kanak/Documents/Endsem Solutions (1).html 4/7

19/12/2022, 19:55 Endsem Solutions

Section B

Question 1

In [16]: p = [np.array([0,0]),np.array([4,0]),np.array([2,3])]

iters = [np.array([2,2])]

for i in range(10000):

iters.append(np.add(iters[-1],p[np.random.randint(0,3)])/2)

In [17]: plt.figure(figsize=(5,5))

plt.scatter([i[0] for i in iters],[i[1] for i in iters],s=0.25)

plt.show()

Question 2

In [18]: A = np.zeros((15,15))

for i in range(1,16):

for j in range(1,16):

A[(i-1,j-1)] = (i**(j-1))%j

plt.imshow(A,cmap='Reds')

<matplotlib.image.AxesImage at 0x7f56a8433280>
Out[18]:

file:///home/kanak/Documents/Endsem Solutions (1).html 5/7

19/12/2022, 19:55 Endsem Solutions

Question 3

In [19]: def sigma(n):

divisors = []

for i in range(1,n+1):

if n%i == 0:

divisors.append(i)

return(sum(divisors))

def S(n,d):

vals = []

for i in range(1,n+1):

if sigma(i)%d == 0:

vals.append(i)

return(sum(vals))

S(20,7)

49
Out[19]:

Question 4

In [20]: def proj(u,v):

return((np.dot(u,v)/np.linalg.norm(v)**2)*v)

In [21]: a = np.array([1,1])

b = np.array([2,0.6])

c = proj(a,b)

In [22]: fig = plt.figure(figsize=(5,5))

ax = plt.axes()

v = plt.arrow(0,0,*a,color='Black',linewidth=0.25)

w = plt.arrow(0,0,*b,color='Black',linewidth=0.25)

z = plt.arrow(0,0,*c,color='Red',linestyle='dashed',linewidth=1)

file:///home/kanak/Documents/Endsem Solutions (1).html 6/7

19/12/2022, 19:55 Endsem Solutions

Question 5

In [23]: f = lambda x: np.cos(x)-x**3

fd = lambda x: -np.sin(x)-3*(x**2)

In [24]: iters = [0.5]

for i in range(10):

iters.append(iters[-1] - f(iters[-1])/fd(iters[-1]))

In [25]: iters

[0.5,

Out[25]:
1.1121416370972725,

0.9096726937368068,

0.8672638182088165,

0.8654771352982646,

0.8654740331109566,

0.8654740331016144,

0.8654740331016144,

0.8654740331016144,

0.8654740331016144,

0.8654740331016144]

In [26]: f(iters[-1])

1.1102230246251565e-16
Out[26]:

file:///home/kanak/Documents/Endsem Solutions (1).html 7/7