ASSIGNMENT 1
ROHIT NILKANTH CHAUDHARI
J18INT640
20/7/2019
BISECTION METHOD
import numpy as np
import [Link] as plt
import math
def f(x):
return (x**3)-(5*x)+1
def bisection(a,b,tol):
xl=a
xr=b
while([Link](xl-xr)>tol):
c=(xl+xr)/2
prod=f(xl)*f(c)
if prod>tol:
xl=c
else:
xr=c
return c
answer= bisection(0,1,0.0005)
print("bisection method gives root at x=",answer)
x= [Link](-10,10,1000)
y=f(x)
[Link](x,y)
[Link]()
[Link]("Plot of Equation")
[Link]("X")
[Link]("F(x)")
[Link](["x"])
OUTPUT
bisection method gives root at x= 0.19775390625
<[Link] at 0xea72c79cf8>
NEWTON RAPHSON METHOD
import numpy as np
from [Link] import derivative
import [Link] as plt
x_n = 5.5
y= [Link]
print(y)
x = [Link](y, y/2, 500)
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)
[Link](x,y,color = "g")
[Link]()
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
CUBIC INTERPOLATION
import numpy as np
from [Link] import CubicSpline
import [Link] as plt
x = [Link](10)
y = [Link](x)
cs = CubicSpline(x,y)
xs = [Link](-0.5,9.6,0.1)
[Link](figsize = (6.5,4))
[Link](x,y,'o',label = 'data')
[Link](xs,[Link](xs),label = 'true')
[Link](xs,cs(xs,1),label = "s")
[Link](xs,cs(xs,2),label = "s'")
[Link](xs,cs(xs,3),label = "s'''")
[Link](-0.5,9.5)
[Link](loc='lower left',ncol=2)
[Link]()
OUTPUT
LU DECOMPOSITION
import numpy as np
import pprint
import [Link]
from [Link] import lu_factor,lu_solve
A = [Link]([[7,5,-1,2],[3,8,1,-4],[-1,1,4,-1],[2,-4,-1,6]])
print(A)
b = [Link]([1,1,1,1])
print(b)
P,L,U = [Link](A)
print("A")
[Link](A)
print("P:")
[Link](P)
lu,piv = lu_factor(A)
x = lu_solve((lu,piv),b)
print(x)
OUTPUT
[[ 7 5 -1 2]
[ 3 8 1 -4]
[-1 1 4 -1]
[ 2 -4 -1 6]]
[1 1 1 1]
A
array([[ 7, 5, -1, 2],
[ 3, 8, 1, -4],
[-1, 1, 4, -1],
[ 2, -4, -1, 6]])
P:
array([[1., 0., 0., 0.],
[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.]])
[-1.09734513 1.20353982 0.00884956 1.33628319]
INTERPOLATE
import numpy as np
from scipy import interpolate
import pylab as py
x = [Link](0, 10, 11)
y = [Link](x)
xnew = [Link](0, 10, 100)
[Link](1)
[Link]('o')
for kind in ['linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic']:
f = interpolate.interp1d(x,y, kind = kind)
ynew = f(xnew)
[Link](xnew, ynew, label = kind)
[Link](loc = 'best')
OUTPUT
TAYLORS SERIES EXPANSION
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
x = [Link]( 0, 0.5, 10 )
f = 1.0 / (1 - x)
f0 = 1 + [Link](10)
f1 = 1 + x
f2 = 1 + x + x**2
[Link]()
[Link](x, f, "k", label = " $f = 1/ (1-x)$")
[Link](x, f, yerr = 0.00005, fmt = "d")
[Link](x, f0, "r", label = " $f - 1$")
[Link](x, f1, "g", label = " $f - 1 + x$")
[Link](x, f2, "--", label = " $f - 1 + x + x**2$")
[Link]("on")
[Link]()
OUTPUT
TRAPAZOIDAL 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.0000253851716194
SIMPSONS 1/3 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)
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
integral1 = h * s
print(integral1)
error1 = (( integral - integral1 )/(integral)) * 100
print ("error1 = ", error1, "%")
OUTPUT
0.999999998453377
0.999999998453377
error1 = 0.0 %
SIMPSONS 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):
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 = 0.0025388007252757916 %
NUMERICAL DIFFERENTIATION
import numpy as np
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)
elif method == "backward":
return (f(a) - f(a - h)) / (h)
else:
return ValueError("method must be central, forward or backward")
derivatives( [Link], 0)
derivatives( [Link], 0, method = "forward")
derivatives( [Link], 0, method = "backward", h = 0.0005)
OUTPUT
0.000249999994705874
import numpy as np
from [Link] import derivative
x = [Link]( 0, 5*[Link], 100 )
dydx = derivative([Link], x)
module = derivative( [Link], x, dx = 0.1 )
dYdx = [Link](x)
dydx
from matplotlib import pyplot as plt
[Link](figsize=(15,5))
[Link](x,dydx,"y",label="central difference", marker="D")
OUTPUT
