Scientific Computational Tools

Name – Harsh R Darji

Roll No. – J18INT601
1. Bisection Method
Programme –
import numpy as np
import matplotlib.pyplot as plt
import math
def f(x):
return (x**3)+1
def bisection(a,b,tol):
xl=a
xr=b
while(np.abs(xl-xr)>tol):
c=(xl+xr)/2
prod=f(xl)*f(c)
if prod>tol:
xl=c
else:
xr=c
return c
answer= bisection(1,-2,0.0005)
print("bisection method gives root at x=",answer)

x= np.linspace(-10,10,1000)
y=f(x)
plt.plot(x,y)
plt.grid()
plt.title("Plot of Equation")
plt.xlabel("X")
plt.ylabel("F(x)")
plt.legend(["x"])

Output - Bisection method gives root at x= -0.9925537109375

1
Scientific Computational Tools

2. Secant Method
Programme –

import numpy as np
f = lambda x: x**3-x-1
def secant_method(f, x0, x1, max_iter=100, tolerance =0.05):
steps_taken = 1
while abs(x1-x0) > tolerance:

x2= (x0*f(x1)-x1*f(x0))/(f(x1) - f(x0))

x1, x0 = x2, x1
steps_taken += 1
max_iter +=1
return x2, steps_taken

root, steps = secant_method(f, 1.3, 1.4)

print ("root is:", root)
print ("steps taken:", steps)

Output –
root is: 1.3246060608507002
steps taken: 3

3. Newton Raphson Method

Programme –

import numpy as np
from scipy.misc import derivative
import matplotlib.pyplot as plt

x_n=5.5
y= np.pi
print(y)
x= np.linspace(y,y/2, 500)

2
Scientific Computational Tools

def f(x):

return x**2-5*x+x/2

def x_next(f, x, x_n):

slope= derivative(f, x_n, dx=0.1)

return x_n-f(x_n)/slope

for n in range(10):

#print(x_n)

print("x_iteration_{} = {}".format(n+1,x_n))

x_n = x_next(f,x,x_n)

y=f(x)
plt.plot(x,y,color = "r")
plt.grid()

Output –

3.141592653589793
x_iteration_1 = 5.5
x_iteration_2 = 4.653846153846149
x_iteration_3 = 4.504923076923077
x_iteration_4 = 4.5000053741714385
x_iteration_5 = 4.500000000006418
x_iteration_6 = 4.499999999999999
x_iteration_7 = 4.5
x_iteration_8 = 4.5
x_iteration_9 = 4.5
x_iteration_10 = 4.5

4. LU Decomposition
Programme –
import numpy as np
import pprint
import scipy
import scipy.linalg # SciPy Linear Algebra Library
from scipy.linalg import lu_factor, lu_solve
A = scipy.array([[1,2,3], [2,3,4], [3,4,1]])

3
Scientific Computational Tools

b = np.array([14,20,14])
print(b)
P, L, U = scipy.linalg.lu(A)

print ("A:")
pprint.pprint(A)

print ("P:")
pprint.pprint(P)

print ("L:")
pprint.pprint(L)

print ("U:")
pprint.pprint(U)

lu, piv = lu_factor(A)

x = lu_solve((lu, piv), b)

print(x)
Output –

[14 20 14]
A:
array([[1, 2, 3],
[2, 3, 4],
[3, 4, 1]])
P:
array([[0., 1., 0.],
[0., 0., 1.],
[1., 0., 0.]])
L:
array([[1. , 0. , 0. ],
[0.33333333, 1. , 0. ],
[0.66666667, 0.5 , 1. ]])
U:
array([[3. , 4. , 1. ],
[0. , 0.66666667, 2.66666667],
[0. , 0. , 2. ]])
[1. 2. 3.]

5. Newton Forward - Backward Interpolation

Programme –
import numpy as np
from scipy import interpolate
from matplotlib import pyplot as plt
x=np.array([5,10,15,20])
y=np.array([20,28,40,80])
f=interpolate.interp1d(x,y)
z5=f(12)
z5
print("Value is ",z5)
plt.plot(x,y)

Output –

4
Scientific Computational Tools

Value is 32.8

6. Lagrange’s Interpolation
a. Programme –
import numpy as np
from scipy import interpolate
from scipy.interpolate import lagrange
import pylab as py
from matplotlib import pyplot as plt

x1=np.array([0,10,11])
y1=np.sin(x1)
poly1=lagrange(x,y)
poly1
print(poly1)
plt.plot(x1,y1)

Output –

2
3 x - 2 x

b. Programme –
import numpy as np
from scipy import interpolate
from scipy.interpolate import lagrange

5
Scientific Computational Tools

import pylab as py
from matplotlib import pyplot as plt
x1=np.array([0,1,2])
y1=x**3
poly1=lagrange(x1,y1)
poly1

Output –
poly1d([ 3., -2., 0.])

7. Spline Interpolation
a. Programme –

import numpy as np
from scipy import interpolate
from scipy.interpolate import lagrange
import pylab as py
from matplotlib import pyplot as plt

x=np.linspace(0,10,11)
y=np.sin(x)
xnew=np.linspace(0,10,100)
py.figure(1)
py.plot(x,y,"ro")

for kind in ['linear','nearest','zero','slinear','quadratic','cubic']:

f=interpolate.interp1d(x,y,kind=kind)
ynew=f(xnew)
py.plot(xnew,ynew,label=kind)
py.legend(loc='lower right')

Output -

6
Scientific Computational Tools

b. Programme –

import numpy as np
from scipy.interpolate import CubicSpline
from matplotlib import pyplot as plt

x=np.arange(10)
y=np.sin(x)

cs=CubicSpline(x,y)
xs=np.arange(-0.5,9.6,0.1)
plt.figure(figsize=(6.5,4))

plt.plot(x,y,'o',label='data')
plt.plot(xs,np.sin(xs),label='true')
plt.plot(xs,cs(xs),label="s")
plt.plot(xs,cs(xs,1),label="s'")
plt.plot(xs,cs(xs,2),label="s''")
plt.plot(xs,cs(xs,3),label="s'''")
plt.xlim(-0.5,9.5)
plt.legend(loc='lower left',ncol=2)
plt.show

Output –

8. Numerical Differentitation

a. Programme –
from scipy.misc import derivative

7
Scientific Computational Tools

import numpy as np
x=np.linspace(0,5*np.pi,100)
dydx=derivative(np.sin,x)

module=derivative(np.sin,x,dx=0.1)
dydx=np.cos(x)
dydx
module
import matplotlib.pyplot as plt
plt.figure(figsize=(15,5))
plt.plot(x,dydx,"g*",label="Central Difference",marker="o")
plt.plot(x,module,"y",label="by
module",alpha=0.4,marker="o",mec="red",mfc="red",markersize=2)
plt.legend(loc="best")

Output –

b. Programme –
import numpy as np
from scipy.misc import derivative
import matplotlib.pyplot as plt

def derivatives(f,a,method="central",h=0.01):
if method == "central":
return (f(a+h)-f(a-h))/(2*h)
elif method == "forward":
return (f(a+h)-f(a-h))/(h)
elif method == "backward":
return (f(a+h)-f(a-h))(h)
else:
raise ValueError ("method must be central,forward or backward")

derivatives(np.exp,0,method="forward",h=0.000025)

Output –

2.0000000002085017

9. Numerical Integration
Programme –

8
Scientific Computational Tools

#Trapezoidal Rule

from math import sin,pi

f=lambda x:x*sin(x)
a=0
b=pi/2
n=90
h=(b-a)/n

s= 0.5*(f(a) + f(b))

for i in range(1,n):
s+=f(a+i*h)
integral = h*s
print(integral)

Output –

1.000025385171619
Programme –
#Simpson's 1/3rd Rule
from math import sin,pi
f=lambda x:x*sin(x)
a=0
b=pi/2
n=90
h=((b-a)/n)
s=(f(a) + f(b))/3

for i in range(1,n,2):
s+=(4*f(a+i*h))/3

for i in range(2,n,2):
s+=(2*f(a+i*h))/3

integral = h*s
print(integral)

Output –

0.999999998453377

Programme –
#Simpson's 3/8 Rule
from math import sin, pi
f = lambda x : x * sin(x)
a = 0
b = pi/2
n = 90
h = (b - a) / n
s = ( f(a) + f(b) )*(3/8)
for i in range( 1, n, 3 ):
s+= (3*f(a + i*h))*(3/8)
for i in range(3, n, 3):
s+= (2*f( a + i*h))*(3/8)
for i in range(2, n, 3):

9
Scientific Computational Tools

s+= (3*f( a + i*h))*(3/8)

integral2 = h * s
print(integral2)
error2 = 100 - (( integral2 )/(integral) * 100)
print ("error2 = ", error2, "%")

Output –

0.9999999965198877
error2 = 1.9334892442657292e-07 %

10. Taylor’s Series Expansion

Programme –
import numpy as np
from matplotlib import pyplot as plt
% matplotlib inline

x=np.linspace(0,0.5,10)
f=1.0/(1-x)

f0=1+np.zeros(10)
f1=1+x
f2=1+x+x**2
f3=1+x+x**2+x**3
print(f0)
print(f1)
print(f2)
print(f3)
plt.figure()
plt.plot(x,f,"k",label="$f=1/(1-x)$")
plt.errorbar(x,f,yerr=0.00005,fmt="d")
plt.plot(x,f0,"r",label="$f=1$")
plt.plot(x,f1,"g",label="$f=1+x$")
plt.plot(x,f0,"b",label="$f=1+x+x**2$")
plt.plot(x,f0,"y",label="$f=1+x+x**2+x**3$")
plt.grid("on")
plt.legend(loc="best")

Output –
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1.05555556 1.11111111 1.16666667 1.22222222 1.27777778
1.33333333 1.38888889 1.44444444 1.5 ]
[1. 1.05864198 1.12345679 1.19444444 1.27160494 1.35493827
1.44444444 1.54012346 1.64197531 1.75 ]
[1. 1.05881344 1.12482853 1.19907407 1.28257888 1.37637174
1.48148148 1.5989369 1.7297668 1.875 ]
<matplotlib.legend.Legend at 0xcd857093c8>

10
Scientific Computational Tools

11