Contents

PART A

1 Least Square Curve Fitting 2

2 Matrix Operations 15

3 Bisection Method 26

4 Runge-Kutta Method 30

5 Numerical Integration 38

6 Newton Raphson Method 44

7 Monte Carlo Method 49

PART B

8 Logistic map 52

9 Two Dimensional Projectile Motion 57

10 Schr̈odinger Equation 65

11 Rutherford Scattering 69

12 Harmonic Oscillator 72

1
Experiment 1

Least Square Curve Fitting

Aim 1
To fit a straight line from the following data and to print the slope and y-intercept correct to
three decimal places.

x 0 0.2 0.4 0.6 0.8 1.0

y 1.357 3.501 5.805 8.012 10.335 12.689

Principle
When the given data is accurate, Interpolation by polynomial fitting gives results. If the given
data has errors as in the case of experimentally generated data, a curve be fitted approximately.
Usually a mathematical equation is fitted to experimental data, by plotting the data on a graph
paper and then passing a straight line through the data [Link] method has obvious drawback
is that the straight line drawn may not be unique. The method of least square is probably the
most systematic procedure to fit a unique curve through the given data point and is widely used
in practical computations. It can be shown that the best fit is obtained if the sum of squares of
deviations from the approximate curve is minimum. This is called least square approximation.
If the curve fitted is a straight line, it is called linear least square approximation.
Let the set of data points be (xi , yi ); i = 1,2,...,m and let the curve given by y = f(x) be fitted
to this data. At x =xi , the experimental (or observed) value of the ordinate is yi and the
corresponding value on the fitting curve is f(xi ). If ei is the error of approximation at x = xi ,
then we have

ei = yi − f (xi ) (1.1)
If we write

S = [y1 − f (x1 )]2 + [y2 − f (x2 )]2 + [y3 − f (x3 )]2 ......[ym − f (xm )]2 = e21 + e22 + e23 .... + e2m (1.2)

Then the method of least squares consists of minimizing S, i.e. the sum of squares of the errors.
Let, y = b + ax be the straight line fitted to the given data. Then corresponding to Eqn 1.2, we
have
S = [y1 − (b + ax1 )]2 + [y2 − (b + ax2 )]2 + ...... + [ym − (b + axm )]2

2
For S to be minimum, we have
∂S
= 0 = −2[y1 − (b + ax1 )] − 2[y2 − (b + ax2 )] − ..... − 2[ym − (b + axm )]
∂b
∂S
= 0 = −2x1 [y1 − (b + ax1 )] − 2x2 [y2 − (b + ax2 )] − ..... − 2xm [ym − (b + axm )]
∂a
The above equation simplifies to

mb + a(x1 + x2 + ..... + xm ) = y1 + y2 + ..... + ym

b(x1 + x2 + .... + xm ) + a(x21 + x22 + .... + x2m ) = x1 y1 + x2 y2 + .... + xm ym

or more compactly to
m
X m
X
mb + a xi = yi
i=i i=i

and m m
X X X
b xi + a x2i = xi y i
i=i i=i

Solving above equations, 2 unknowns a and b can be obtained as follows:

P P P
n xi y i − xi y i
a =
n x2i − xi xi
P P P

n x2i
P P P P
yi − xi xi yi
b =
n x2i − xi xi
P P P

Algorithm
Step 1: Start.

Step 2: Input x values.

Step 3: Input y values.

Step 4: Create a blank list, xx

Step 5: Create a blank list, xy

Step 6: n=len(x)

Step 7: For i in range(n), do steps 8 and 9

Step 8: [Link](x[i]**2)

Step 9: [Link](x[i]*y[i])

Step 10: sumx=sum(x)

3
Step 11: sumy=sum(y)
Step 12: sumxx=sum(xx)
Step 13: sumxy=sum(xy)
Step 14: denom=(n*sumxx)-(sumx*sumx)
Step 15: a=(n*sumxy-sumx*sumy)/denom
Step 16: b=(n*sumxx-sumx*sumxy)/denom
Step 17: Print a and b
Step 18: Stop

Program
from pylab import *
x= l i s t ( e v a l ( i n p u t ( ’ Enter x=v a l u e s ’ ) ) )
y= l i s t ( e v a l ( i n p u t ( ’ Enter y=v a l u e s ’ ) ) )
xx = [ ]
xy = [ ]
n=l e n ( x )
f o r i i n range ( n ) :
xx . append ( x [ i ] * * 2 )
xy . append ( x [ i ] * y [ i ] )
sumx=sum ( x )
sumy=sum ( y )
sumxx=sum ( xx )
sumxy=sum ( xy )
denom=(n * sumxx) = (sumx * sumx )
a=(n * sumxy=sumx * sumy ) / denom
b=(sumy * sumxx=sumx * sumxy ) / denom
p r i n t ( ’ S l o p e i s : %0.3 f and y= i n t e r c e p t i s :%0.3 f ’ % ( a , b ) )
scatter (x , y)
x f i t=l i n s p a c e ( x [ 0 ] , x [ = 1 ] , 1 0 0 )
y f i t=a * x f i t + b
t i t l e ( ’ Least square f i t ’ )
xlabel ( ’x ’)
ylabel ( ’y ’)
plot ( xfit , y f i t )
show ( )

Output
Enter x= v a l u e s 0 . 0 , 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 , 1 . 0
Enter y= v a l u e s 1 . 3 5 7 , 3 . 5 0 1 , 5 . 8 0 5 , 8 . 0 1 2 , 1 0 . 3 3 5 , 1 2 . 6 8 9
S l o p e i s : 1 1 . 3 3 8 and y= i n t e r c e p t i s : 1 . 2 8 1

4
5
Aim 2
The force F required to pull load L are tabulated [Link] F and L obey a linear relationship,find
the force needed to pull a load L=150kg.
F(N) 12 15 21 25 29
L(kg) 50 70 100 120 140

Algorithm
Step 1: Start.

Step 2: Input x values.

Step 3: Input y values.

Step 4: Create a blank list, xx

Step 5: Create a blank list, xy

Step 6: n=len(x)

Step 7: For i in range(n), do steps 8 and 9

Step 8: [Link](x[i]**2)

Step 9: [Link](x[i]*y[i])

Step 10: sumx=sum(x)

Step 11: sumy=sum(y)

Step 12: sumxx=sum(xx)

Step 13: sumxy=sum(xy)

Step 14: denom = (n ∗ sumxx) − (sumx ∗ sumx)

Step 15: a = (n ∗ sumxy − sumx ∗ sumy)/denom

Step 16: b = (n ∗ sumxx − sumx ∗ sumxy)/denom

Step 17: Define x0 = (y1 − b)/a

Step 18: Print values of y at x=x0

Step 19: Plot xfit versus yfit curve with Force on x-axis and Load on y-axis

Step 20: Stop

6
Program
from pylab import *
x= l i s t ( e v a l ( i n p u t ( ’ Enter F=v a l u e s ’ ) ) )
y= l i s t ( e v a l ( i n p u t ( ’ Enter L=v a l u e s ’ ) ) )
y1=e v a l ( i n p u t ( ’ Enter th e v a l u e o f ”L” a t which ”F” i s c a l c u l a t e d : ’))
xx = [ ]
xy = [ ]
n=l e n ( x )
f o r i i n range ( n ) :
xx . append ( x [ i ] * * 2 )
xy . append ( x [ i ] * y [ i ] )
sumx=sum ( x )
sumy=sum ( y )
sumxx=sum ( xx )
sumxy=sum ( xy )
denom=(n * sumxx) = (sumx * sumx )
a=(n * sumxy=sumx * sumy ) / denom
b=(sumy * sumxx=sumx * sumxy ) / denom
x0=(y1=b ) / a
p r i n t ( ’ t he v a l u e s o f y a t F= ’ , y1 , ’ i s ’ , x0 )
scatter (x , y)
x f i t=l i n s p a c e ( x [ 0 ] , x [ = 1 ] , 1 0 0 )
y f i t=a * x f i t + b
t i t l e ( ’ Least square f i t ’ )
x l a b e l ( ’ Force ’ )
y l a b e l ( ’ Load ’ )
plot ( xfit , y f i t )
show ( )

Output
Enter F=v a l u e s 1 2 ,15 ,21 ,25 ,29
Enter L=v a l u e s 5 0 ,70 ,100 ,120 ,140
Enter t he v a l u e o f ”L” a t which ”F” i s c a l c u l a t e d : 150
the values of y a t F= 150 i s 3 0 . 7 5 4 4 2 0 4 3 2 2 2 0 0 4

7
8
Aim 3
The voltage across a capacitor at time t sec. is given by the following table. Use the principle of
least squares to fit an exponential curve to the data.

t(sec) 0 2 4 6 8
V (volts) 100 60 22 13 8

Principle
As the voltage across the capacitor decreases with time, it represents discharge of a capacitor
through a resistor. The governing equation connecting v(t) and t is non-linear. However the
equation can be bought to linear form.

q(t) = q0 ∗ e−t/τ

q(t)/C = q0 /C ∗ e−t/τ
v(t) = v0 ∗ e−t/τ
1
ln(v) = ln(v0 ) −
∗t
τ
which is the equation of a straight line with slope = -1/τ and y-intercept = ln(v0 )
Obtain the slope and y-intercept by fitting a straight line to ln(v) versus t curve. Use these
values to plot the curve using the above formula.

Algorithm
Step 1: Start.

Step 2: import numpy as np

Step 3: import pylab.

Step 4: Define t and v values.

Step 5: Define v1 as log(v).

Step 6: Define n as length of v.

Step 7: Define p for finding sum of t and q for finding sum of v1.

Step 8: Define r for finding sum of squares of t and s for finding sum of product of t and v1.
p∗q−n∗s
Step 9: slope =
p ∗ ∗2 − n ∗ r
p∗s−r∗q
Step 10: intercept =
p ∗ ∗2 − n ∗ r
Step 11: Define a new array t1 for new t values.

9
Step 12: Define :
vcurve = eintercept+slope∗t1

Step 13: Print Time constant.

Step 14: Print voltage at t = 0.

Step 15: plot t1 versus vcurve with time along x-axis and voltage along y-axis.

Step 16: Stop.

Program
import numpy as np
from pylab import *
t=np . a ra nge ( 0 , 9 , 2 )
v=np . a r r a y ( [ 1 0 0 , 6 0 , 2 2 , 1 3 , 8 ] )
v1=np . l o g ( v )
n=l e n ( v )
p , q=sum ( t ) , sum ( v1 )
r , s=sum ( t * * 2 ) , sum ( t * v1 )
s l o p e =(p * q=n * s ) / ( p **2 = n * r )
i n t e r c e p t=np . a v e r a g e ( v1)= s l o p e * np . a v e r a g e ( t )
t 1=np . l i n s p a c e ( 0 , 8 , 1 0 0 0 )
vc ur ve=np . exp ( i n t e r c e p t+s l o p e * t 1 )
p r i n t ( ’ Time c o n s t a n t =%0.6 f s e c ’% abs ( 1 . 0 / s l o p e ) )
p r i n t ( ’ v o l t a g e a t t=0 i s :%0.6 f v o l t s ’%np . exp ( i n t e r c e p t ) )
p l o t ( t1 , vcurve )
p l o t ( t , v , ’ ro ’ )
x l a b e l ( ’ time ’ )
ylabel ( ’ voltage ’ )
show ( )

Output
Time c o n s t a n t =3.039120 s e c
v o l t a g e a t t=0 i s : 9 9 . 7 9 9 1 9 5 v o l t s

10
11
Aim 4
Fit a parabola, y=a*x2 +b*x+c to the following set of observations. Print values of a, b and c.

x 0 1 2 3 4
y 1 6 11 21 39

Principle
The given data is to be fitted to a parabola given by the equation,y=a*x2 +b*x+c. So the error
function,
E = Σ[yi − (a ∗ x2i + b ∗ xi + c)]2
For E to be minimum,
∂E
=0
∂a
∂E
=0
∂b
∂E
=0
∂c
∂E
= 0 implies
∂a
a ∗ Σ(x4i ) + b ∗ Σ(x3i ) + c ∗ Σ(x2i ) = Σ(x2i ∗ yi )
∂E
= 0 implies
∂b
a ∗ Σ(x3i ) + b ∗ Σ(x2i ) + c ∗ Σ(xi ) = Σ(xi ∗ yi )
∂E
=0
∂c
a ∗ Σ(x2i ) + b ∗ Σ(xi ) + n ∗ c = Σ(yi )
Solving the above equations by forming arrays with the coefficients of a,b and c we could find
the values of a,b and c. Let p = Σ(x4i ), q = Σ(x3i ), r = Σ(x2i ), s = Σ(xi ), t = Σ(x), u = Σ(xi ∗ yi ),
v = Σ(yi )      
p q r a t
q r s  ∗  b  = u
r s n c v
Solving the above matrix a, b and c can be found.

Algorithm
Step 1: Start.

Step 2: Input x and y values.

Step 3: Define p for finding sum of x values to the power of 4.

Step 4: Define q for finding sum of cube of x values.

12
 Step 5: Define r for finding sum of square of x values.

Step 6: Define s for finding sum of x values.

Step 7: Define t for finding sum of product of squares of x and y values.

Step 8: Define u for finding sum of product of x and y.

Step 9: Define v for finding sum of y values.

Step 10: Define m as array of coefficient of a,b and c and n the array of t,u and v.

Step 11: Print the fit parameters a,b and c.

Step 12: Assign xp and yp as the new values of x and y.

Step 13: Plot xp and yp with x on x-axis and y on y-axis.

Step 14: Stop.

Program
# python pgm o f p a r a b o l i c f i t u s i n g l e a s t s q u a r e s method
# y = ax ˆ2 + bx + c
from numpy import *
from numpy . l i n a l g import *
from pylab import *

x=ar an ge ( 0 , 5 , 1 ) # Given data p o i n t s

y=a r r a y ( [ 1 , 6 , 1 1 , 2 1 , 3 9 ] ) # Given data p o i n t s

p=sum ( x * * 4)
q=sum ( x * * 3)
r=sum ( x * * 2)
s=sum ( x )

t=sum ( x ** 2 * y )
u=sum ( x * y )
v=sum ( y )

m=a r r a y ( [ [ p , q , r ] , [ q , r , s ] , [ r , s , 1 ] ] ) # C o e f f i c i e n t matrix
n=a r r a y ( [ [ t ] , [ u ] , [ v ] ] )
o=s o l v e (m, n )
a=o [ 0 ]
b=o [ 1 ]
c=o [ 2 ]

p r i n t ( ’ F i t p a r a m e t e r s a r e : a = ’ , a , ’ b= ’ ,b , ’ c = ’ , c )

13
xp=ara ng e ( x [ 0 ] , x [ = 1 ] , 0 . 0 1 ) # x array for fitting
yp=a * xp ** 2 + b * xp +c # y array for fitting
s c a t t e r ( x , y , l a b e l =’ Given data ’ ) # To p l o t g i v e n data
p l o t ( xp , yp , l a b e l =’ P a r a b o l i c f i t ’ ) # To p l o t th e f i t t e d c u r v e
xlabel ( ’x ’)
ylabel ( ’y ’)
legend ()
show ( )

# End o f pgm

Output
F i t p a r a m e t e r s a r e : a= [ 1 . 8 0 3 3 7 0 7 9 ] b= [ 2 . 4 6 1 7 9 7 7 5 ] c= [ = 0 . 7 1 9 1 0 1 1 2 ]

14
Experiment 2

Matrix Operations

Aim 1
Write a program to read two matrices from the keyboard and find their sum and product using
core python.

Principle
Matrix is an array of numbers arranged into rows and columns. Matrix of m rows and n columns
is said to be of the order, (mxn). Numpy module in python supports the operation on compound
data types like arrays and matrices.

ADDITION OF TWO MATRICES

Two matrices A and B of the orders,(mxn) and (p,q) can be added if and only if m=p and
n=q. Adding two matrices gives another matrix of the same order.
If A+B=C, then A,B and C are of the same order.

MULTIPLICATION OF TWO MATRICES

Two matrices can be multiplied only if the number of columns of the first matrix is equal to the
number of rows of the second matrix.
Consider two matrices A and B of the orders (mxn) and (pxq) respectively. A and can e mul-
tiplied only if n=p. The resulting matrix will be of the order (mxq).We can find the matrix
multiplication by the dot() python code. dot(A,B) is used to find the product of two matrices A
and B.

TRANSPOSE OF A MATRIX

Transpose of a matrix is obtained by interchanging the rows and columns of the original matrix.
The element at ith row and jth column in A becomes in jth row and ith column in AT . transpose()
command is used to find the transpose of a matrix in pyhton.

TRACE OF A MATRIX

15
The trace of a matrix is the sum of the diagonal elements. It is possible to find the trace of
a matrix only if the number of rows and columns are equal. Thus, trace is only valid for square
matrices. trace() command is used to find trace in python.

INVERSE OF A MATRIX

If A.B = B.A = I, then the matrix B is called the inverse of matrix A, denoted by A−1 . A
matrix whose determinant is zero is called a singular matrix.
For finding the inverse of a matrix,
(i) The matrix should be non-singular.
(ii)It should be a square matrix.
The inverse of a matrix A is found by,

adjA
A−1 = ; —A— ̸= 0
|A|
We use det() to find the determinant of the matrix and linalg inv() to find the inverse, where
both are found in the [Link] sub package.

Algorithm
Step 1: Start.

Step 2: Import numpy module.

Step 3: Input the number of rows as m1 and coloumns as n1 of the first matrix.

Step 4: Input the matrix elements as x.

Step 5: y=array(x).reshape(m1,n1)

Step 6: Print y.

Step 7: Input the number of rows as m2 and coloumns as n2 of second matrix.

Step 8: Input the matrix elements as x.

Step 9: z=array(x).reshape(m2,n2)

Step 10: Print z.

Step 11: if m1==m2, do step 12 to step 15, else do step 16

Step 12: if n1==n2, do step 13 and 14, else do step 15

Step 13: c=y+z

Step 14: Print c

Step 15: print(’Summation not possible’)

16
Step 16: print(’the summation is not possible’)

Step 17: if n1==m2, do step 18 and 19, else do step 20

Step 18: d=dot(y,z)

Step 19: Print d

Step 20: print(’multiplication is not possible’)

Step 21: Stop.

Program
from numpy import *
m1, n1=e v a l ( i n p u t ( ’ e n t e r t h e number o f rows and columns o f f i r s t
matrix s e p a r a t e d by comma ’ ) )
x=e v a l ( i n p u t ( ’ e n t e r th e e l e m e n t s o f f i r s t matrix ’ ) )
y=a r r a y ( x ) . r e s h a p e (m1, n1 )
print (y)
m2, n2=e v a l ( i n p u t ( ’ e n t e r t h e number o f rows and columns o f second
matrix s e p a r a t e d by comma ’ ) )
x=e v a l ( i n p u t ( ’ e n t e r th e e l e m e n t s o f second matrix ’ ) )
z=a r r a y ( x ) . r e s h a p e (m2, n2 )
print ( z )
i f (m1==m2) :
i f ( n1==n2 ) :
c=y+z
p r i n t ( ’ th e sum o f two matrix i s ’ , c )
else :
p r i n t ( ’ th e summation i s not p o s s i b l e ’ )
else :
p r i n t ( ’ th e summation i s not p o s s i b l e ’ )
i f ( n1==m2) :
d=dot ( y , z )
p r i n t ( ’ th e product o f two m a t r i c e s i s ’ , d )
else :
p r i n t ( ’ m u l t i p l i c a t i o n i s not p o s s i b l e ’ )

Output
e n t e r t h e number o f rows and columns o f f i r s t matrix s e p a r a t e d by
comma2 , 2
e n t e r t h e e l e m e n t s o f f i r s t matrix1 , 2 , 3 , 4
[ [ 1 2]
[3 4 ] ]
e n t e r t h e number o f rows and columns o f second matrix s e p a r a t e d by
comma2 , 2

17
e n t e r t h e e l e m e n t s o f second matrix3 , 4 , 2 , 6
[ [ 3 4]
[2 6 ] ]
t h e sum o f two matrix i s [ [ 4 6 ]
[ 5 10]]
t h e product o f two m a t r i c e s i s [ [ 7 16]
[17 36]]

18
Aim 2
Write a program to read a matrix from the keyboard and print its trace and transpose using core
python.

Algorithm
Step 1: Start

Step 2: import numpy module.

Step 3: Input the number of rows and coloumns of the matrix as m1 and n1

Step 4: Input the matrix elements as x.

Step 5: y=array(x).reshape(m1,n1)

Step 6: Print y.

Step 7: if m1==n1, do step 7 to step 11, else do step 12 and 13.

Step 8: p=[Link]()

Step 9: Print p.

Step 10: t=trace(y)

Step 11: Print t

Step 12: print(’Transpose of the matrix cannot be be found as the given matrix is not a square
matrix.’)

Step 13: print(’Trace cannot be found as the given matrix is not a square matrix.’)

Step 14: Stop

Program
from numpy import *
m1, n1=e v a l ( i n p u t ( ’ e n t e r t h e number o f rows and columns o f t h e matrix
s e p a r a t e d by comma ’ ) )
x=e v a l ( i n p u t ( ’ e n t e r th e e l e m e n t s o f t h e matrix ’ ) )
y=a r r a y ( x ) . r e s h a p e (m1, n1 )
print (y)
i f (m1==n1 ) :
p=y . t r a n s p o s e ( )
p r i n t ( ’ Transpose o f th e g i v e n matrix i s ’ , p )
t=t r a c e ( y )
p r i n t ( ’ The t r a c e o f t h e g i v e n matrix i s ’ , t )
else :

19
 p r i n t ( ’ Transpose o f th e matrix cannot be be found as th e g i v e n
matrix i s not a s q u a r e matrix . ’ )
p r i n t ( ’ Trace cannot be found as t he g i v e n matrix i s not a s q u a r e
matrix . ’ )

Output
e n t e r t h e number o f rows and columns o f th e matrix s e p a r a t e d by
comma2 , 2
e n t e r t h e e l e m e n t s o f t h e matrix1 , 4 , 1 2 , 3
[ [ 1 4]
[12 3 ] ]
Transpose o f t he g i v e n matrix i s [ [ 1 1 2 ]
[ 4 3]]
The t r a c e o f t he g i v e n matrix i s 4

20
Aim 3
Write a program to read a matrix from the keyboard and print its inverse using core python.

Algorithm
Step 1: Start

Step 2: import numpy module.

Step 3: import [Link]

Step 4: Input the number of rows and coloumns of the matrix as m1 and n1

Step 5: Input the matrix elements as x.

Step 6: y=array(x).reshape(m1,n1)

Step 7: Print y.

Step 8: if m1==n1, do step 9 to 12, else do step 13.

Step 9: if det(y)==0, do step 10, else do step 11 and 12.

Step 10: print(’The given matrix is singular and hence inverse cannot be obtained’)

Step 11: yinv=inv(y)

Step 12: Print yinv

Step 13: print(’The inverse of the given matrix cannot be obtained as it is not a square matrix.’)

Step 14: Stop

Program
from numpy import *
from numpy . l i n a l g import *
m1, n1=e v a l ( i n p u t ( ’ e n t e r t h e number o f rows and columns o f t h e matrix
s e p a r a t e d by comma ’ ) )
x=e v a l ( i n p u t ( ’ e n t e r th e e l e m e n t s o f t h e matrix ’ ) )
y=a r r a y ( x ) . r e s h a p e (m1, n1 )
print (y)
i f m1==n1 :
i f de t ( y ) ==0:
p r i n t ( ’ The g i v e n matrix i s s i n g u l a r and hence i n v e r s e cannot
be o b t a i n e d . ’ )
else :
yinv=i n v ( y )
p r i n t ( ’ The i n v e r s e o f th e g i v e n matrix i s ’ , yi nv )

21
else :
p r i n t ( ’ The i n v e r s e o f th e g i v e n matrix cannot be found as i t i s
not a s q u a r e matrix . ’ )

Output
e n t e r t h e number o f rows and columns o f th e matrix s e p a r a t e d by
comma2 , 2
e n t e r t h e e l e m e n t s o f t h e matrix6 , 5 , 3 , 4
[ [ 6 5]
[3 4 ] ]
The i n v e r s e o f th e g i v e n matrix i s [ [ 0 . 4 4 4 4 4 4 4 4 = 0.55555556]
[ = 0.33333333 0 . 6 6 6 6 6 6 6 7 ] ]

22
Aim 4
Write a program to read N×N matrix A from the keyboard and print its trace t, transpose AT
and inverse A−1 using relevant methods defined in NUMPY. Also print AAT , AA−1 and A−1 AT .

Algorithm
Step 1: Start

Step 2: Import numpy as np.

Step 3: Input the number of rows and columns as m,n

Step 4: Read the row values as a list.

Step 5: print matrix

Step 6: Initialise a for loop

Step 7: print A

Step 8: A←[Link](A)

Step 9: print trace of matrix.

Step 10: Atrans←[Link](A)

Step 11: print transpose.

Step 12: Initialise a for loop.

Step 13: print transpose of matrix.

Step 14: Ainv←[Link](A)

Step 15: print inverse.

Step 16: Initialise a for loop.

Step 17: print inverse of matrix.

Step 18: AAtrans←[Link]([Link])

Step 19: print A × AT

Step 20: Initialise a for loop.

Step 21: print AAT

Step 22: AAinv←[Link]([Link])

Step 23: Initialise a for loop.

Step 24: print A × A−1

23
Step 25: Initialise a for loop.

Step 26: print AAinv[i]

Step 27: AinvAtrans←[Link](Ainv,Atrans)

Step 28: print A−1 AT .

Step 29: Initialise a for loop.

Step 30: print AinvAtrans[i]

Step 31: Stop

Program
import numpy as np
m, n=e v a l ( i n p u t ( ’ Enter t he number o f rows and columns : ’ ) )
A=[ l i s t ( e v a l ( i n p u t ( ’ Enter row ’+ s t r ( i +1)+ ’ : ’ ) ) ) f o r i i n range (m) ]
p r i n t ( ’ The matrix i s : ’ )
f o r i i n range (m) :
p r i n t (A[ i ] )
A=np . a r r a y (A)
p r i n t ( ’ Trace i s : ’ , np . t r a c e (A) )
Atrans=np . t r a n s p o s e (A)
p r i n t ( ’ Transpose i s : ’ )
f o r i i n range (m) :
p r i n t ( Atrans [ i ] )
Ainv=np . l i n a l g . i n v (A)
print ( ’ Inverse i s : ’)
f o r i i n range (m) :
p r i n t ( Ainv [ i ] )
AAtrans=np . dot (A, Atrans )
p r i n t ( ’A*AˆT i s : ’ )
f o r i i n range (m) :
p r i n t ( AAtrans [ i ] )
AAinv=np . dot (A, Ainv )
p r i n t ( ’A*Aˆ=1 i s : ’ )
f o r i i n range (m) :
p r i n t ( AAinv [ i ] )
AinvAtrans=np . dot ( Ainv , Atrans )
p r i n t ( ’Aˆ=1 AˆT i s : ’ )
f o r i i n range (m) :
p r i n t ( AinvAtrans [ i ] )

Output

24
Enter t he number o f rows and columns : 3 , 3
Enter row1 : 1 , 0 , 2
Enter row2 : 0 , 3 , 0
Enter row3 : 4 , 0 , 5
The matrix i s :
[1 , 0 , 2]
[0 , 3 , 0]
[4 , 0 , 5]
Trace i s : 9
Transpose i s :
[1 0 4]
[0 3 0]
[2 0 5]
Inverse i s :
[ = 1.66666667 0 . 0.66666667]
[0. 0.33333333 0. ]
[ 1.33333333 0. = 0.33333333]
A*AˆT i s :
[ 5 0 14]
[0 9 0]
[14 0 41]
A*Aˆ=1 i s :
[ 1 . 0. 0 . ]
[ 0 . 1. 0 . ]
[ 0 . 0. 1 . ]
Aˆ=1 AˆT i s :
[ = 0.33333333 0 . = 3.33333333]
[ 0 . 1. 0 . ]
[0.66666667 0. 3.66666667]

25
Experiment 3

Bisection Method

Aim 1
Using bisection method, find the real root of the equation x3 − 2x − 5 = 0.

Principle
Bisection method is the simplest among all the numerical schemes to solve the transcendental
equations. This scheme is based on intermediate value theorem for continuous functions. The
theorem states that if f (a) is continuous between a and b and f (a) and f (b) are of opposite
signs, then there exists at least one root between a and b. For definiteness, let f (a) be negative
and f (b) be positive, then the root lies between a and b and let the approximate value be given
(a + b)
by x0 =
2
If f (x0 ) = 0, then x0 is the root of equation. If f (x0 ) is either positive or negative, the root lies
between either x0 and a or between x0 and b respectively. The new interval is designated as[a1,
|b − a|
b1], whose length is
2
As before this is bisected at x1 and the new interval will be exactly half the length of previous
one. The process is repeated until the latest interval containing the root is as small as desired,
say r. The interval width is reduced by a factor one half at each step and at the end of nth step,
|b − a| b−a
the new interval will be [an,bn] of length n
.Then we have | n |, which on simplification
2 2
(|b − a|)
loge
gives, n≥ ϵ .
loge 2
The method can be represented graphically as shown in figure 3.1.
If there are more roots than one in the interval bisection method finds only one of the roots. It
can be easily programmed using the following computational steps:
1. Choose 2 real numbers a and b such thatf (a) ∗ f (b) < 0.
(a + b)
2. Set c = .
2
3. Iff (a) ∗ f (c) < 0 the roots lies in the interval (a,c). Then set b = c and go to step 2.
4. Iff (a) ∗ f (c) are greater than zero , then the root lies in the interval (c,b) . Then seta = cand
go to step 2.

26
 Figure 3.1: Bisection Method

5. Iff (a) ∗ f (c) = 0 it means that c is the root of the equation f (x) = 0and computation maybe
terminated.

Algorithm
Step 1: Start.

Step 2: Define a function and return its value.

Step 3: Define average of a function x and return its value.

Step 4: Read the upper and lower limits.

Step 5: While abs(a − b) >= allowed error value

Step 6: if f (x) ∗ f (a) > 0

Step 7: set a = x

Step 8: else b = x

Step 9: Print the solution, x.

Step 10: else print , no root exist in between this limit.

Step 11: Stop.

Program
def f (x ) :
r e t u r n x ** 3 = 2 * x = 5.0
def average (a , b ) :
r e t u r n ( a+b ) / 2 . 0
a=f l o a t ( i n p u t ( ’ e n t e r t h e upper l i m i t a : ’))
b=f l o a t ( i n p u t ( ’ e n t e r t he upper l i m i t b : ’))
i =0

27
w h i l e abs ( a=b) >1.0 e = 8:
i f f ( a v e r a g e ( a , b ) ) * f ( a ) <0:
b=a v e r a g e ( a , b )
else :
a=a v e r a g e ( a , b )
i=i +1
p r i n t ( ’ The r o o t a f t e r %d i t e r a t i o n s i s %10.7 f ’%( i , a ) )
p r i n t ( ’ f (%0.6 f ) = %0.3 e ’ % ( a , f ( a ) ) )

Output
e n t e r t h e upper l i m i t a : 2 . 0
e n t e r t h e upper l i m i t b : 3 . 0
The r o o t a f t e r 1 i t e r a t i o n s i s 2.0625000
The r o o t a f t e r 2 i t e r a t i o n s i s 2.0937500
The r o o t a f t e r 3 i t e r a t i o n s i s 2.0942383
The r o o t a f t e r 4 i t e r a t i o n s i s 2.0944824
The r o o t a f t e r 5 i t e r a t i o n s i s 2.0945435
The r o o t a f t e r 6 i t e r a t i o n s i s 2.0945511
The r o o t a f t e r 7 i t e r a t i o n s i s 2.0945513
The r o o t a f t e r 8 i t e r a t i o n s i s 2.0945514
The r o o t a f t e r 9 i t e r a t i o n s i s 2.0945515
The r o o t a f t e r 10 i t e r a t i o n s i s 2.0945515
f ( 2 . 0 9 4 5 5 1 ) = = 2.632 e =09

28
Aim 2
Write a Python Program to find the positive root of the equation x ∗ ex = 1 using bisection
method.

Program
from math import exp
def f (x ) :
r e t u r n x * exp ( x) =1
def average (a , b ) :
r e t u r n ( a+b ) / 2 . 0
a=f l o a t ( i n p u t ( ’ Enter t h e upper l i m i t a : ’ ) )
b=f l o a t ( i n p u t ( ’ Enter t he l o w e r l i m i t b : ’ ) )
i =0
w h i l e abs ( a=b) >1.0 e = 8:
i f f ( a v e r a g e ( a , b ) ) * f ( a ) <0:
b=a v e r a g e ( a , b )
else :
a=a v e r a g e ( a , b )
i=i +1
p r i n t ( ’ The r o o t a f t e r %d i t e r a t i o n s i s %10.7 f ’%( i , a ) )
p r i n t ( ’ f (%0.6 f ) = %0.3 e ’ % ( a , f ( a ) ) )

Output
Enter t he upper l i m i t a : 0 . 0
Enter t he l o w e r l i m i t b : 0 . 7
The r o o t a f t e r 1 i t e r a t i o n s i s 0.3500000
The r o o t a f t e r 2 i t e r a t i o n s i s 0.5250000
The r o o t a f t e r 3 i t e r a t i o n s i s 0.5468750
The r o o t a f t e r 4 i t e r a t i o n s i s 0.5578125
The r o o t a f t e r 5 i t e r a t i o n s i s 0.5632812
The r o o t a f t e r 6 i t e r a t i o n s i s 0.5660156
The r o o t a f t e r 7 i t e r a t i o n s i s 0.5666992
The r o o t a f t e r 8 i t e r a t i o n s i s 0.5670410
The r o o t a f t e r 9 i t e r a t i o n s i s 0.5671265
The r o o t a f t e r 10 i t e r a t i o n s i s 0.5671371
The r o o t a f t e r 11 i t e r a t i o n s i s 0.5671425
The r o o t a f t e r 12 i t e r a t i o n s i s 0.5671432
The r o o t a f t e r 13 i t e r a t i o n s i s 0.5671432
The r o o t a f t e r 14 i t e r a t i o n s i s 0.5671433
The r o o t a f t e r 15 i t e r a t i o n s i s 0.5671433
f ( 0 . 5 6 7 1 4 3 ) = = 1.837 e =09

29
Experiment 4

Runge-Kutta Method

Aim 1
To find the solution of first order differential equations using Runge-Kutta method.

Principle
Euler’s method is less efficient in practical problems since it requires h to be small for obtaining
reasonable accuracy. The Runge-Kutta methods are designed to give greater accuracy and they
possess the advantages of requiring only the function values at some selected points on the
interval.
Runge - Kutta methods refers to a family one step methods for numerical solutions of differential
equations. Runge-Kutta methods are known by their order. Euler’s method is a first order R-K
method. Similarly we have second order R-k, fourth order R-K etc. In R-K family, as order
increases accuracy also increases.
In the Euler’s method,

Yi+1 = Yi + slope × stepsize

Yi+1 = Yi + mh

Here m is the slope of a function at a point.

In Euler’s method, we are considering the slope of a single point. That is the reason for deviation
in the answer. Instead of taking the slope at one point, we can measure the slope at different
points. The average of this can be used for further calculation. Order of R-K method is exactly
, the number of slopes which we are using for the calculation of next step. So Euler’s method is
a first order method.
In second order method we have to consider the slope at two points. One slope is from the
previous point and the other from the forth-coming point by extrapolation. The average of these
two can be used to find the final solutions. In the fourth order R-K there will be four slopes.
The fourth order Runge- Kutta formula is,
1
Y1 = Y0 + × (K1 + 2K2 + 2K3 + K4 )
6

30
where,

K1 = hf (X0 , Y0 )
h K1
K2 = hf (X0 + , Y0 + )
2 2
h k2
K3 = hf (X0 + , Y0 + )
2 2
K4 = hf (X0 + h, Y0 + K3 )
1
Y1 = Y0 + (K1 + 2K2 + 2K3 + K4 )
6

Algorithm
Step 1: Start.
Y
Step 2: Define the function −1.0 ×
(1.0 + X)
Step 3: Read X,Y,h,X0 .
(X0 − X)
Step 4: p ← (1 + )
h
Step 5: Assign an interval range for X

Step 6: Initialize for loop

Step 7: print X,Y

Step 8: K1 = hf (X, Y )
K1
Step 9: K2 = hf (X + h2 , Y + 2
)
K2
Step 10: k3 = hf (X + h2 , Y + 2
)

Step 11: K4 = hf (X + h, Y + K3 )

Step 12: Y1 = Y + 16 (K1 + 2K2 + 2K3 + K4 )

Step 13: End for loop

Step 14: Stop.

31
Program
import numpy as np
d e f f ( x , y ) : r e t u r n = 1.0 * ( y /(1.0+ x ) )
x , y=e v a l ( i n p u t ( ’ Enter th e i n i t i a l c o n d i t i o n ( x0 , y0 ) : ’ ) )
h=e v a l ( i n p u t ( ’ Enter s t e p s i z e : ’ ) )
x0=e v a l ( i n p u t ( ’ Value a t which th e s o l u t i o n i s needed ( x0 ) : ’))
p=(1+( f l o a t ( x0=x ) / h ) )
x=np . l i n s p a c e ( x , x0 , i n t ( p ) )
for x in x :
p r i n t ( ( ’ x=%0.3 f y=%0.6 f ’ %(x , y ) ) )
k1=h * f ( x , y )
k2=h * f ( x+(h / 2 ) , y+(k1 / 2 ) )
k3=h * f ( x+(h / 2 ) , y+(k2 / 2 ) )
k4=h * f ( x+h , y+k3 )
y=y+(( k1 +(2 * k2 )+(2 * k3)+k4 ) / 6 . 0 )

Output
Enter t he i n i t i a l c o n d i t i o n ( x0 , y0 ) : 0 . 0 , 2 . 0
Enter s t e p s i z e : 0 . 2 5
Value a t which t h e s o l u t i o n i s needed ( x0 ) : 2 . 5
x =0.000 y =2.000000
x =0.250 y =1.600000
x =0.500 y =1.333333
x =0.750 y =1.142857
x =1.000 y =1.000000
x =1.250 y =0.888889
x =1.500 y =0.800000
x =1.750 y =0.727273
x =2.000 y =0.666667
x =2.250 y =0.615385
x =2.500 y =0.571429

32
Aim 2
Solve the differential equation using RK4 method; dy/dx=x+y2 given that y(0)=[Link] y(0.45)

Algorithm
Step 1: Start.

Step 2: Import numpy.

Step 3: Define function.

Step 4: Enter the initial condition as (x0,y0)

Step 5: Enter the step size.

Step 6: Enter the value at which solution is needed as x0.

Step 7: calculate k1 =h*f(x,y).

Step 8: Calculate k2 =h*f(x+h/2,y+k1 /2)

Step 9: Calculate k3 =h*f(x+h/2,y+k2 /2).

Step 10: Calculate k4 =h*f(x+h,y+k3 ).

Step 11: Calculate y=y+((k1 +2*k2 +2*k3 +k4 )/6.0).

Step 12: x=x+h.

Step 13: Print the solution.

Step 14: Stop.

Program
import numpy as np
d e f f ( x , y ) : r e t u r n x+y * * 2 . 0
x , y=i n p u t ( ’ e n t e r th e i n i t i a l c o n d i t i o n ( x0 , y0 ) : ’ )
h=i n p u t ( ’ e n t e r s t e p s i z e : ’ )
x0=i n p u t ( ’ v a l u e a t which th e s o l u t i o n i s needed ( x0 ) : ’)
p r i n t ’ x=%0.3 f y=%0.6 f ’ % ( x , y )
w h i l e x<=x0 :
k1=h * f ( x , y )
k2=h * f ( x+h / 2 , y+k1 / 2)
k3=h * f ( x+h / 2 , y+k2 / 2)
k4=h * f ( x+h , y+k3 )
y=y+(( k1+2* k2+2* k3+k4 ) / 6 . 0 )
x=x+h
p r i n t ’ x=%0.3 f y=%0.6 f ’ %(x , y )

33
Output
e n t e r t h e i n i t i a l c o n d i t i o n ( x0 , y0 ) : 0 . 0 , 1 . 0
enter step s i z e : 0.03
v a l u e a t which t h e s o l u t i o n i s needed ( x0 ) : . 4 5
x =0.000 y =1.000000
x =0.030 y =1.031387
x =0.060 y =1.065708
x =0.090 y =1.103225
x =0.120 y =1.144242
x =0.150 y =1.189107
x =0.180 y =1.238226
x =0.210 y =1.292070
x =0.240 y =1.351193
x =0.270 y =1.416250
x =0.300 y =1.488022
x =0.330 y =1.567444
x =0.360 y =1.655651
x =0.390 y =1.754025
x =0.420 y =1.864281
x =0.450 y =1.988555

34
Aim 3
d2 y dy 2
Solve the second order differential equation dx2
- x( dx ) + y 2 =0 using RK4 method with y(0)=1
and y′(0)=0 with h=0.2

Principle
The Runge-kutta methods are designed to give greater accuracy and they possess the advantages
of requiring only the function values at some selected points on the interval.

kx1 , kv1 = h ∗ pos(t, x, v), h ∗ velo(t, x, v)

kx2 , kv2 = h ∗ pos(t + h/2, x + kx1 /2, v + kv1 /2), h ∗ velo(t + h/2, x + kx1 /2, v + kv1 /2)

kx3 , kv3 = h ∗ pos(t + h/2, x + kx2 /2, v + kv2 /2), h ∗ velo(t + h/2, x + kx2 /2, v + kv2 /2)

kx4 , kv4 = h ∗ pos(t + h, x + kx3 , v + kv3 ), h ∗ velo(t + h, x + kx3 , v + kv3 )

x = x + (kx1 + 2 ∗ kx2 + 2 ∗ kx3 + kx4 )/6.0

v = v + (kv1 + 2 ∗ kv2 + 2 ∗ kv3 + kv4 )/6.0

Algorithm
Step 1: Import numpy

Step 2: Import pi,sin

Step 3: Define position

Step 4: Define velocity

Step 5: Enter the initial condition

Step 6: Enter the value at which solution is needed

Step 7: Enter the step size

Step 8: calculate kx1 , kv1 = h ∗ pos(t, x, v), h ∗ velo(t, x, v)

Step 9: calculate kx2 , kv2 = h ∗ pos(t + h/2, x + kx1 /2, v + kv1 /2), h ∗ velo(t + h/2, x + kx1 /2, v +
kv1 /2)

35
Step 10: calculate kx3 , kv3 = h ∗ pos(t + h/2, x + kx2 /2, v + kv2 /2), h ∗ velo(t + h/2, x + kx2 /2, v +
kv2 /2)

Step 11: calculate kx4 , kv4 = h ∗ pos(t + h, x + kx3 , v + kv3 ), h ∗ velo(t + h, x + kx3 , v + kv3 )

Step 12: calculate x = x + (kx1 + 2 ∗ kx2 + 2 ∗ kx3 + kx4 )/6.0

Step 13: calculate v = v + (kv1 + 2 ∗ kv2 + 2 ∗ kv3 + kv4 )/6.0

Program
from pylab import *
import numpy as np
from math import pi , s i n
d e f pos ( t , x , v ) : r e t u r n v
d e f v e l o ( t , x , v ) : r e t u r n 1 . 0 * t * v **2 = x ** 2
x , v=i n p u t ( ’ e n t e r x [ t 0 ] and x prime [ t 0 ] : ’ )
x0=f l o a t ( i n p u t ( ’ v a l u e a t which th e s o l u t i o n i s needed ( t 0 ) : ’ ) )
t=f l o a t ( i n p u t ( ’ th e v a l u e o f ”x” a t which t h e i n i t i a l c o n d i t i o n i s g i v e n : ’))
h=i n p u t ( ’ s t e p s i z e : ’ )
p=1+( f l o a t ( x0=t ) / h )
g=np . l i n s p a c e ( 0 , x0 , p )
for t in g :
p r i n t ’ x=%0.3 f y=%0.6 f ’ % ( t , x )
kx1 , kv1=h * pos ( t , x , v ) , h * v e l o ( t , x , v )
kx2 , kv2=h * pos ( t+h / 2 , x+kx1 / 2 , v+kv1 / 2 ) , h * v e l o ( t+h / 2 , x+kx1 / 2 , v+kv1 /2 )
kx3 , kv3=h * pos ( t+h / 2 , x+kx2 / 2 , v+kv2 / 2 ) , h * v e l o ( t+h / 2 , x+kx2 / 2 , v+kv2 /2 )
kx4 , kv4=h * pos ( t+h , x+kx3 , v+kv3 ) , h * v e l o ( t+h , x+kx3 , v+kv3 )
x=x+(kx1+2* kx2+2* kx3+kx4 ) / 6 . 0
v=v+(kv1+2* kv2+2* kv3+kv4 ) / 6 . 0

Output
e n t e r x [ t 0 ] and x prime [ t 0 ] : 1 . 0 , 0 . 0
v a l u e a t which t h e s o l u t i o n i s needed ( t 0 ) : 1 . 0
t h e v a l u e o f ”x” a t which th e i n i t i a l c o n d i t i o n i s g i v e n : 0 . 0
step s i z e : 0.05
x =0.000 y =1.000000
x =0.050 y =0.990417
x =0.100 y =0.978939
x =0.150 y =0.965652
x =0.200 y =0.950657
x =0.250 y =0.934063
x =0.300 y =0.915995
x =0.350 y =0.896583
x =0.400 y =0.875969
x =0.450 y =0.854297

36
x =0.500 y =0.831717
x =0.550 y =0.808379
x =0.600 y =0.784431
x =0.650 y =0.760018
x =0.700 y =0.735278
x =0.750 y =0.710344
x =0.800 y =0.685336
x =0.850 y =0.660366
x =0.900 y =0.635534
x =0.950 y =0.610930
x =1.000 y =0.586631

37
Experiment 5

Numerical Integration

Aim 1
R1 1
Integrate 0
dx using Simpson’s 1/3 rule and Trapezoidal rule and compare the results..
1+x

Principle
Let we have some data points (x1 , y1 ), (x2 , y2 ), ..., (xn , yn ) connected by a function y = f(x), where
f(x)is not known explicitly. Numerical integration is a method to compute the definite integral
Rb
of the function by using the observed data points, a f (x)dx. Divide the interval [a,b] into n
subintervals of equal width, such that, x0 = a and xn = b
Z b
I= f (x)dx
a

Approximating y by Newton’s forward difference formula, we obtain,

Z xn
p(p − 1) 2 p(p − 1)(p − 2) 3
I= [y0 + p∆y0 + ∆ y0 + ∆ y0 + ......]dx (5.1)
x0 2 6

since x = x0 + ph, dx = hdp and the above integral becomes,

Z n
p(p − 1) 2 p(p − 1)(p − 2) 3
I=h [y0 + p∆y0 + ∆ y0 + ∆ y0 + ......]dp (5.2)
0 2 6
which gives on simplification,
Z xn Z xn
n n(2n − 3) 2 n(n − 2)2 3
ydx = nh [y0 + ∆y0 + ∆ y0 + ∆ y0 + ......]dx (5.3)
x0 x0 2 12 24

From this general formula, we can obtain different integration formulae by putting n = 1, 2,
3,....etc. We derive here a few of these formulae but it should be remarked that the Trapezoidal
and Simpson’s 1/3 rules are found to give sufficient accuracy for use in practical problems.

TRAPEZOIDAL RULE

38
For n=1 the first order approximation is the Trapezoidal Rule. We can neglect the higher
order terms, so the equation becomes,
Z xn
n
I= ydx = h[y0 + ∆y0 ]
x0 2

For
first interval [x0 , x1 ] the equation becomes,
Z x1
I= ydx = h[y0 + y1 ]
x0

For the next interval [x1 , x2 ] we deduce similarly,

Z x2
I= ydx = h[y1 + y2 ]
x1

And so on...... Combining all these expression, we obtain the rule,

Z xn
I= ydx = h[y0 + 2(y1 + y2 + .... + yn−1 ) + yn ]
x0

This is known as Trapezoidal Rule.

SIMPSONS RULE:

In Simpson’s 1/3rd rule, the interval [a,b] is divided into n sub elements. Since two elements are
coupled for analysis, n = 2. Then,
Z x2
1 h
ydx = 2h(y0 + ∆y0 + ∆2 y0 ) = (y0 + 4y1 + y2 )
x0 6 3

Similarly, Z x4
h
ydx = (y2 + 4y3 + y4 )
x2 3
And finally, Z xn−4
h
ydx = (yn−2 + 4yn−1 + yn )
xn−1 3
Summing we obtain,
Z xn
h
ydx = [y0 + 4(y1 + y3 + ... + yn−1 ) + 2(y2 + y4 .... + yn−2 ) + yn ] (5.4)
x0 3

Algorithm
Step 1: Start.
Step 2: import math
Step 3: define the empty data sets

39
 Step 4: initialise the sum (tot)

Step 5: enter the number of datasets

Step 6: initialize the ‘for‘ loop

Step 7: x=a+i*h

Step 8: append the x data value to the empty x dataset

Step 9: append y data values to the empty y data sets

Step 10: end the for loop.

Step 11: step size h=xdata[1]-xdata[0]

Step 12: total=ydata[0]+ydata[n-1]

Step 13: initialize for loop

Step 14: total=total+4*ydata end of for loop

Step 15: initialize for loop

Step 16: total=total+2*ydata

Step 17: integral=total*h/3

Step 18: print integral

Step 19: Stop.

Program
from math import s q r t , cos , p i
def f (x ) : return sqrt ( cos (x ))
p r i n t (’======== T r a p e z o i d a l r u l e ========’)
a , b , n=e v a l ( i n p u t ( ’ I n t e g r a t i o n l i m i t s ( a , b ) and no . o f i n t e r v a l s ”n” : ’))
h=f l o a t ( b=a ) / n
y=[ f ( a+i * h ) f o r i i n range ( n +1)]
I n t e g r a l =(y [ 0 ] + y[ = 1]+2 * sum ( y [ 1 : = 1 ] ) ) * h / 2 . 0
p r i n t ( ’ Value o f i n t e g r a l by T r a p e z o i d a l r u l e = %0.7 f ’ % I n t e g r a l )
p r i n t (’======== Simpson 1/3 rd r u l e ========’)
a , b , n=e v a l ( i n p u t ( ’ I n t e g r a t i o n l i m i t s ( a , b ) and no . o f i n t e r v a l s ”n” : ’))
h=f l o a t ( b=a ) / n
y=[ f ( a+i * h ) f o r i i n range ( n +1)]
I n t e g r a l =(y [ 0 ] + y[ = 1]+4 * sum ( y [ 1 : = 1 : 2 ] ) + 2 * sum ( y [ 2 : = 1 : 2 ] ) ) * h / 3 . 0
p r i n t ( ’ Value o f i n t e g r a l by Simpsons r u l e = %0.7 f ’ % I n t e g r a l )

Output

40
======== T r a p e z o i d a l r u l e ========
I n t e g r a t i o n l i m i t s ( a , b ) and no . o f i n t e r v a l s ”n” : 0 . 0 , p i /2 ,1 000
Value o f i n t e g r a l by T r a p e z o i d a l r u l e = 1 . 1 9 8 1 2 7 3
======== Simpson 1/3 rd r u l e ========
I n t e g r a t i o n l i m i t s ( a , b ) and no . o f i n t e r v a l s ”n” : 0 . 0 , p i /2 ,1 000
Value o f i n t e g r a l by Simpsons r u l e = 1 . 1 9 8 1 3 5 2

41
Aim 2
A particle starts from rest and attains instantaneous velocities as given in the table. Find the
approximate distance covered in 20 seconds using Simpson’s 1/3 rule.

t(s) 2 4 6 8 10 12 14 16 18 20
v(m/s) 10 18 25 29 32 20 11 5 2 0

Algorithm
Step 1: Start

Step 2: Input the integration limits a, b and the step size as h.

Step 3: Input the given v values

Step 4: Integral=(v[0]+v[-1]+4*sum(v[1:-1:2])+2*sum(v[2:-1:2])*h/3.0

Step 5: print Integral

Step 6: Stop

Program
a , b , h=e v a l ( i n p u t ( ’ I n t e g r a t i o n l i m i t s ( a , b ) and s t e p s i z e ( h ) : ’))
v=[10.0 ,18.0 ,25.0 ,29.0 ,32.0 ,20.0 ,11.0 ,5.0 ,2.0 ,0.0]
I n t e g r a l =(v [ 0 ] + v[ = 1]+4 * sum ( v [ 1 : = 1 : 2 ] ) + 2 * sum ( v [ 2 : = 1 : 2 ] ) ) * h / 3 . 0
p r i n t ( ’ D i s t a n c e t r a v e l l e d i n 2 0 . 0 s = %0.3 f ’ % I n t e g r a l )

Output
I n t e g r a t i o n l i m i t s ( a , b ) and s t e p s i z e ( h ) : 2 . 0 , 2 0 . 0 , 2 . 0
Distance t r a v e l l e d in 20.0 s = 292.000

42
Aim 3
A curve is drawn to pass through the following points. Find the volume of solid generated by
revolving the area bounded by the curve, x-axis and ordinates at x=0 and x=1 about x-axis.

x 0.00 0.25 0.50 0.75 1.00 1.25

y 1.0 0.9896 0.9589 0.9089 0.8415 0.8850

Algorithm
Step 1: Start

Step 2: Import numpy.

Step 3: Import pylab.

Step 4: Input x values as an array.

Step 5: Input y values as an array.

Step 6: y=y*y

Step 7: h=x[1]-x[0]

Step 8: V=[Link]*(h/3.0)*(y[0]+y[-1]+4*sum(y[1:-1:2])+2*sum(y[0:-1:2 ] ) )

Step 9: Print V

Step 10: Stop

Program
import numpy as np
from pylab import *
x=np . a r r a y ( [ 0 . 0 , 0 . 2 5 , 0 . 5 , 0 . 7 5 , 1 . 0 ] )
y=np . a r r a y ( [ 1 . 0 , 0 . 9 8 9 6 , 0 . 9 5 8 9 , 0 . 9 0 8 9 , 0 . 8 4 1 5 ] )
y=y * y
h=x [1] = x [ 0 ]
V=np . p i * ( h / 3 . 0 ) * ( y [ 0 ] + y[ = 1]+4 * sum ( y [ 1 : = 1 : 2 ] ) + 2 * sum ( y [ 0 : = 1 : 2 ] ) )
p r i n t ( ’ Volume g e n e r a t e d : %0.5 f ’ % V)

Output
Volume g e n e r a t e d : 3 . 3 4 2 8 5

43
Experiment 6

Newton Raphson Method

Aim
x
To obtain the solution of the given equation sin x = using Newton-Raphson method.
2

Principle
Newton-Raphson method is a more accurate method to improve the approximate root of alge-
braic or transcendental equations, obtained by some other methods which are not so accurate.
Let x0 be an approximate or guessed root of the equation f (x) = 0. Then the correct root
x1 = x0 + h, so that f (x1 ) = 0. We can assume that h is the correction. Substituting for x,

f (x0 + h) = 0.
Expanding with Taylor series,

h ′ h2 ′′ h3 ′′′
f (x0 ) + f (x0 ) + f (x0 ) + f (x0 ) + ... = 0.
1! 2! 3!
Neglecting the higher order terms, we have

f (x0 ) + hf ′ (x0 ) = 0
f (x0 )
h=− ′
f (x0 )

Now the approximate root can be converted into a better root x1 by,

x1 = x0 + h,
f (x0 )
x1 = x0 – ′ .
f (x0 )

44
 if x1 is not an root with desired accuracy level, improve it again by

f (x1 )
x2 = x1 – ,
f ′ (x1 )
f (x2 )
x3 = x2 – ′ ,
f (x2 )
...

In general,
f (xn )
xn+1 = xn – . (6.1)
f ′ (xn )
This can be continued till we get the root with the desired accuracy. The procedure given by
equation 6.1 is known as Newton-Raphson Method.

Algorithm
Step 1: Start.
Step 2: import sine and cosine from math module
Step 3: Define the function: f (x) = sin x − 0.5 ∗ x (as per the given function).
Step 4: Define derivative of the function and return its value: g(x) = cos x − 0.5
Step 5: Input the initial guess as assigning x
Step 6: introduce while loop with the condition f (x) > 1.0e−8
f (x)
Step 7: Calculate x = x − .
g(x)
Step 8: Print the solution, x and f(x).
Step 9: Stop.

Program
from math import s i n , c o s
d e f f ( x ) : r e t u r n s i n ( x ) = 0.5 * x
d e f g ( x ) : r e t u r n c o s ( x ) = 0.5
x=e v a l ( i n p u t ( ’ Enter i n i t i a l g u e s s : ’ ) )
w h i l e abs ( f ( x )) >1.0 e = 8:
x=x= f ( x ) / g ( x )
p r i n t ( ’ The r o o t i s : %0.6 f and f (%0.6 f ) = %0.3 e ’ % ( x , x , f ( x ) ) )

Output
Enter i n i t i a l g u e s s : 2
The r o o t i s : 1 . 8 9 5 4 9 4 and f ( 1 . 8 9 5 4 9 4 ) = = 1.431 e =10

45
Aim 2
Using Newton-Raphson method, find any two minima of the function f (x) = sin(x). ln(x) be-
tween x = 0 and x = 20

Algorithm
Step 1: Start.

Step 2: Import pylab

Step 3: Define the function: original function.

Step 4: Define function f(z).

Step 5: Define function g(z)

Step 6: xp ← linspace(0,20,500)

Step 7: yp ← foriginal(xp)

Step 8: Plot graph with xp in x axis and yp in y axis.

Step 9: i←1

Step 10: Initialise while loop

Step 11: Flag ← False

Step 12: Initialise while loop.

Step 13: Read x

Step 14: Initialise while loop.

Step 15: abs(f(x)) > 1.0e−7

f (x)
Step 16: x←x-
g(x)
Step 17: Print x, f(x), g(x).

Step 18: Print maxima or minima.

Step 19: Conditional statement: if g(x) > 0 do step 20 and 21 otherwise, do step 22.

Step 20: Print g(x) and minima.

Step 21: Flag ← True

Step 22: Print not a minima.

Step 23: increment i

Step 24: Stop.

46
Program
from pylab import *
def f o r i g i n a l ( z ) : return sin ( z )* log ( z )
d e f f ( z ) : r e t u r n ( s i n ( z ) / z ) +( l o g ( z ) * c o s ( z ) )
d e f g ( z ) : r e t u r n =( s i n ( z ) / z * * 2) +(2 * c o s ( z ) / z )= s i n ( z ) * l o g ( z )

xp=l i n s p a c e ( 0 , 2 0 , 5 0 0 )
yp=f o r i g i n a l ( xp )
p l o t ( xp , yp )
show ( )

i =1
w h i l e ( i <=2) :
f l a g =’ F a l s e ’
w h i l e f l a g ==’F a l s e ’ :
x=f l o a t ( i n p u t ( ’ I n i t i a l g u e s s o f th e r o o t : ’) )
w h i l e ( abs ( f ( x ) ) >1.0 e = 7) :
x=x= f ( x ) /g ( x )

p r i n t ( ”To f i n d minima %d , \ t x= %5.3 f , \ t f ( x )= %10.3 e \ t g (

x )= %10.3 e”%( i , x , f ( x ) , g ( x ) ) )
p r i n t ( ’ One o f t he minima/maxima i s a t x = %0.5 f ’ % x )
i f g ( x ) >0:
p r i n t ( ’ S i n c e g ( x ) = %0.5 f i s p o s i t i v e , i t i s i n d e e d a
minima o f t he g i v e n f u n c t i o n \n’% g ( x ) )
f l a g =’True ’
else :
p r i n t ( ’ I t i s not a minima , r e p e a t t h e run with a
d i f f e r e n t i n i t i a l g u e s s \n ’ )
i=i +1

Output
Output 1:
I n i t i a l guess of the root : 4
I t e r a t i o n s= 1 , x= 4 . 8 4 3 , f ( x )= 9 . 7 5 3 e =09
One o f t he minima i s a t x = 4 . 8 4 2 5 6
S i n c e g ( x ) = 1 . 6 5 9 9 9 i s p o s i t i v e , i t i s i n d e e d a minima o f t h e g i v e n
function
Output 2:
I n i t i a l g u e s s o f t h e r o o t : 10
I t e r a t i o n s= 1 , x= 1 1 . 0 3 3 , f ( x )= 1 . 7 8 9 e =14
One o f t he minima i s a t x = 1 1 . 0 3 3 3 1
S i n c e g ( x ) = 2 . 4 1 4 2 6 i s p o s i t i v e , i t i s i n d e e d a minima o f t h e g i v e n
function

47
48
Experiment 7

Monte Carlo Method

Aim
To find the value of π using Monte Carlo algorithm. Tabulate the value for every 20000 interval
of iteration steps from 100000 to 200000.

Principle
In Monte Carlo method we use repeat random sampling to obtain numerical results. One of the
basic examples of getting started with the Monte Carlo algorithm is the estimation of Pi. The
idea is to simulate random (x,y) points in a 2-D plane with domain as a square of side 1 unit.
Imagine a circle inside the same domain with same diameter and inscribed into the square. We
then calculate the ratio of number points that lied inside the circle and total number of generated
points.
We know that area of the square is 1 unit sq while that of circle is π ∗ ( 21 )2 = π4 . Now for a very
large number of generated points,

Figure 7.1: montecarlo fig

49
 area of the circle no. of points generated inside the circle
=
area of the square total no. of points generated or no. of points generated inside the square
that is,
no. of points generated inside the circle
π=4
no. of points generated inside the square
Generate random (x, y) pairs and check if x2 + y 2 ≤ 1 . If yes, we increment the number of
points that appears inside the [Link] more the number of iterations, the more accurate the
result is. Thus, value of π can be calculated.

Algorithm
Step 1: Start.
Step 2: Import random and math modules
Step 3: Assign initial and final value
Step 4: Assign value of interval in between iteration
Step 5: Assigning initial count as zero
Step 6: For loop start for the range of values that to be generated till final value
Step 7: Conditional statement if to check the generated randoms are in inside circle
Step 8: Conditional statement if, for every iteration steps starting from initial value having a
interval
cnt
Step 9: Print value of π using 4.0 ∗ for every value i that generated
i
Step 10: Stop

Program
from random import *
from math import *
i n i v a l =100000
f i n v a l =200000
i n t =20000
c n t=0
f o r i i n range ( 1 , f i n v a l +1):
i f ( random () ** 2+ random () ** 2) <=1:
c nt+=1
i f ( i>=i n i v a l and ( i )% i n t ==0):
p r i n t ( ’ Value o f p i= ’ , 4 . 0 * c nt / i , ’ f o r v a l u e o f i= ’ , i )

50
Output
Value of p i= 3 . 1 4 3 4 4 f o r v a l u e o f i= 100000
Value of p i= 3 . 1 4 0 3 3 3 3 3 3 3 3 3 3 3 3 4 f o r v a l u e o f i= 120000
Value of p i= 3 . 1 3 9 0 8 5 7 1 4 2 8 5 7 1 4 3 f o r v a l u e o f i= 140000
Value of p i= 3 . 1 4 0 1 5 f o r v a l u e o f i= 160000
Value of p i= 3 . 1 3 8 8 4 4 4 4 4 4 4 4 4 4 4 5 f o r v a l u e o f i= 180000
Value of p i= 3 . 1 4 0 6 f o r v a l u e o f i= 200000

51
Experiment 8

Logistic map

Aim
To plot the set of fixed points x against the control parameter k of the Logistic map function
xn+1 = kxn (1 − xn ), 1 ≤ k ≤ 4, 0 ≤ xn ≤ 1.

Principle
Logistic map is a function which relates the co-ordinates of a point Pn+1 in terms of those of the
previous points Pn . The map is always discrete as it uses the previous values of the dependent
variable as the present value of independent [Link] is thus no question of differentiability
for a [Link] logistic map was developed by Robert May in 1876 as a mathematical model of
population growth whose generations do not overlap with a fixed environment. It is given by,

xn+1 = kxn (1 − xn ).

This is called Logistic map where 0 ≤ x ≤ 1 and k ≥ 1.

A continuous form of this map is the logistic equation,

f (x) = kx(1 − x).

Characteristics of Logistic equation and map:

1. The roots of the logistic equation are obtained by setting f (x) = 0.

ie., kx(1 − x) = 0 =⇒ x = 0, 1.

df (x)
2. Extremum occurs when = 0.
dx
f (x) = kx(1 − x) = k(x − x2 ),
df (x)
= k(1 − 2x),
dx
df (x) 1
= 0 at x = .
dx 2

52
 d2 f (x)
This is a maximum because = −2k < 0. The point x = 1/2 is called the critical
dx2
point of the function possessing only single maximum at a given interval, 0 < x < 1 in this
case.

3. After some iterations of the map, it often converges to some fixed value called an ‘attractor’.
Any further iteration of it will yield the same value. If x∗n is such a value, then

x∗n = kx∗n (1 − x∗n ),

1 = k(1 − x∗n ),
1 1
= (1 − x∗n ) =⇒ x∗n = 1 − .
k k

Condition for stability of the attractor:

If xn < 0, then iterations will lead xn+1 → −∞. If xn = 0, then xn+1 = 0 always.
For x = k1 , xn+1 = 1 − k1 = x∗n . The range 0 < xn < k1 is called the ‘Basin of attraction’ of x∗n .
A value xn approaches x∗n if successive iterations bring it closer to x∗n .

xn+1 − x∗n
< 1,
xn − x∗n
f (xn ) − x∗n
< 1.
f (xn−1 ) − x∗n

When f (xn−1 ) − x∗n → 0,

df (xn )
< 1,
dxn
xn =x∗n
1 2
|k(1 − 2x∗n )| = k(1 − 2[1 − ]) = k(1 − 2 + ) = |k − 2k + 2| = |2 − k| < 1
k k

This is possible only if 1 < k < 3. When k = 3, the attractor bifurcates to two fixed points x∗1
and x∗2 in such a way that,
x∗1 = f (x∗2 ) and x∗2 = f (x∗1 )
x∗2 = f (f (x∗2 ))
= f (kx∗2 (1 − x∗2 ))
= k(kx∗2 (1 − x∗2 ))[1 − kx∗2 (1 − x∗2 )]
= k 2 x∗2 (1 − x∗2 )[1 − kx∗2 (1 − x∗2 )]
Each x2 is said to be a fixed point of period 2. In general, if xp is a fixed point of period p, xp
repeats after a set of p iteration of f. That is,

f (p) (xp ) = f (f (f (. . . . . . . . . p times(xp )))) = xp

53
This bifurcation of the attractor at k = 3 is called pitchfork bifurcation due to its shape.

|f (f (xn ))| ≤ 1
√
This requires k ≥ 1 + 6 = 3.449489743. Then each branch of fixed points bifurcates into two
separate branches. The points on these branches will be of period 4.
If k is further increased, further branching occurs. Fixed points of period p give rise to 2p
branches. It is found that for k = 3.5699 . . . , an infinite number of bifurcations occur. In logistic
map, fixed points never repeat. The band of fixed points forms a continuum. Complete chaos
begins from this point. Thus bifurcation is the route to chaos for logistic equation.

Algorithm
Step 1: Start.

Step 2: Input cstart ,cend , stability point and [Link] iteration , n.

Step 3: Step size is (cend- cstart) /n

Step 4: Initialise c = cstart

Step 5: Star iteration while c < cend

Step 6: For i from 0 to n , compute x=c x(1-x)

Step 7: Append x value and c value

Step 8: Increment c value

Step 9: Continue till c < cend

Step 10: Plot commands are given

Step 11: End

Program
from pylab import *
c s t a r t=f l o a t ( i n p u t ( ’ e n t e r upper range o f c o n t r o l parameter ’ ) )

cend=f l o a t ( i n p u t ( ’ e n t e r l o w e r range o f c o n t r o l parameter ’ ) )

s t b=i n p u t ( ’ s t a b i l i t y p o i n t ’ )
s t b=i n t ( s t b )
n=f l o a t ( i n p u t ( ’ e n t e r t he number o f i t e r a t i o n ’ ) )
h=(cend = c s t a r t ) / n
c l i s t =[]
x l i s t =[]
c=c s t a r t

54
w h i l e c<cend :
x =0.1
f o r i i n range ( 0 , s t b ) :
x=c * x * (1 = x )
f o r i i n range ( 0 , i n t ( n ) ) :
x=c * x * (1 = x )
c l i s t . append ( c )
x l i s t . append ( x )
c=c+h
x l a b e l ( ’ c o n t r o l parameter ’ )
y l a b e l ( ’ population ’ )
t i t l e ( ’ L o g i s t i c Map ’ )
plot ( clist , xlist , ’ . ’ )

show ( )

Output
e n t e r upper range o f c o n t r o l parameter0
e n t e r l o w e r range o f c o n t r o l parameter4
s t a b i l i t y point300
e n t e r t h e number o f i t e r a t i o n 1 0 0 0

To plot the function for a single value of k.

Algorithm
Step 1: Start.

Step 2: Import pylab.

55
 Step 3: Define f (k, x) as kx(1 − x).

Step 4: Initialize variables x0 = 0.1, x = x0 and the number of iterations, n.

Step 5: Set initial values of the lists Li and Lx as 0 and x0 respectively.

Step 6: In a for loop repeat steps 8 to 10 n times.

Step 7: Evaluate x = f (k, x).

Step 8: Append i to list Li ans x to list Lx.

Step 9: Plot Li versus Lx.

Step 10: Label x and y axes .

Step 11: Show.

Step 12: Stop.

Program

#L o g i s t i c map= e v o l u t i o n o f x f o r g i v e n k
from pylab import *
def f (k , x ) :
r e t u r n k * x * (1 = x )
x0 =0.1
x=x0
k =3.45
n=1000
Li , Lx = [ 0 ] , [ x0 ]
f o r i i n range ( 1 , n ) :
x=f ( k , x )
Li . append ( i )
Lx . append ( x )
xlabel ( ’n ’)
y l a b e l ( ’ xn ’ )
p l o t ( Li , Lx , ’ . ’ )
show ( )
#end program

56
Experiment 9

Two Dimensional Projectile Motion

Aim
To plot the trajectory of a projectile that moving near earth’s surface and compare the trajectories
with and without considering air resistance. Also evaluate the time of flight, maximum horizontal
range and vertical height in both cases.

Principle
Two dimensional projectile motion is a freely falling body that is projected near earth surface
and follows a paraboloic path.
Here we consider the trajectories with and without considering air drag

PROJECTILE MOTION - WITHOUT CONSIDERING AIR DRAG

Consider a body thrown at an angle θ with an initial velocity v0 . The body will go up first
and eventually fall back on the ground.
Consider the influence of gravity only, then

Fx = 0
Fy = −mg
ax = 0
ay = −g

Since the projectile at an angle,the initial velocity can be split into two components. They are
v0 cosθ along x axis and v0 sinθ along y axis. Then by Euler method, along x-direction,

Acceleration = 0
Initial velocity = v0 cos θ
vx(i+1) = vxi + haxi
x(i+1) = xi + hvxi

57
and along y-direction,

Acceleration = − g
Initial velocity = v0 sinθ
vy(i+1) = vyi + hayi
y(i+1) = yi + hvyi

Using these formulas, the position and velocity at any stage can be [Link] the case of
a projectile, the maximum value of displacement in x-direction is called maximum range. The
maximum value of displacement in y direction is called maximum height.

PROJECTILE MOTION - WITH CONSIDERING AIR DRAG

Consider the motion of a projectile when moving through air. All the parameters will be changed
due to the viscous drag force of air.
The airdrag,
1
Fd = Cπρr2 v 2 (9.1)
2
The effect of air drag can be split into two components Fd cosϕ along x axis and Fd sinϕ along y
axis,where ϕ is the angle made by velocity component with x direction at any instant.
Along x-direction

Force due to Earth’s gravity = 0

1
Air drag = − Fd cosϕ = − Cπρr2 v 2 cosϕ
2

Cπρr2 v 2 cosϕ
Net acceleration at any instant = − = −kv 2 cosϕ
2m
where,

1
k= Cπρr2
2m
vx
At any instant, cosϕ =
v
where,
q
v = vx2 + vy2 (9.2)
vx q
ax = −kv 2 = −kvvx = −kvx vx2 + vy2 (9.3)
v
(9.4)

58
Along y-direction

Force due to Earth’s gravity = -mg

1
Air drag = −Fd sinϕ = − Cπρr2 v 2 sinϕ
2

Cπρr2 v 2 sinϕ
Net acceleration at any instant = −g − = −g − kv 2 sinϕ
2m
where,
1
k= Cπρr2
2m
vy
At any instant, sinϕ =
q v
where, v = vx + vy2
2

vy q
ay = −g − kv 2 = −g − kvvy = −g − kvy vx2 + vy2
v
From acceleration we can estimate the velocity and position along x axis and y axis at any instant
by Euler method.

q
ax(i+1) = −kvx vx2 + vy2
q
ay(i+1) = −g − kvy vx2 + vy2
vx(i+1) = vxi + haxi
vy(i+1) = vyi + hayi
x(i+1) = xi + hvxi
y(i+1) = yi + hvyi

TRAJECTORY WITHOUT CONSIDERING AIR DRAG

Algorithm
Step 1: Start.

Step 2: Import pylab and numpy modules.

Step 3: Initialise x,y

Step 4: Initialise t and set step size, dt as 0.001

Step 5: Give velocity of projection and angle of projection.

59
 Step 6: Find out the velocity in x direction and y direction.

Step 7: Give acceleration due to gravity in x and y directions.

Step 8: Initialize the array for time, velocity and position.

Step 9: Introduce while loop with the condition expression y≥ 0.0.

Step 10: Give the equations to find acceleration, velocity and position inside the loop.

Step 11: Append the values of velocity, position and time inside the loop .

Step 12: Outside the loop

find maximum time, range and height.

Step 13: Print all the above data and plot the graph between x and y position.

Step 14: Show.

Step 15: Stop.

Program
from numpy import *
from pylab import *
x , y =0.0 ,0.0
t , dt =0 ,0.001
u=e v a l ( i n p u t ( ’ Enter t he v e l o c i t y o f p r o j e c t i o n ’ ) )
ang=e v a l ( i n p u t ( ’ Enter t he a n g l e o f p r o j e c t i o n i n d e g r e e s : ’ ) )
vx=u * c o s ( r a d i a n s ( ang ) )
vy=u * s i n ( r a d i a n s ( ang ) )
ax =0.0
ay= =9.8
x1 , y1 , vx1 , vy1 , t 1 =[x ] , [ y ] , [ vx ] , [ vy ] , [ t ]
w h i l e y >=0.0:
vy=vy+ay * dt
x=x+vx * dt
y=y+vy * dt
t=t+dt
vy1 . append ( vy )
x1 . append ( x )
y1 . append ( y )
t 1 . append ( t )
p r i n t ( ’ Time o f f l i g h t = %8.3 f s ’%max( t 1 ) )
p r i n t ( ’ H o r i z o n t a l range = %8.3 f m’%max( x1 ) )
p r i n t ( ’ Maximum h e i g h t = %8.3 f m’%max( y1 ) )
p l o t ( x1 , y1 )
xlabel ( ’x ’)
ylabel ( ’y ’)

60
t i t l e ( ’ P r o j e c t i l e Motion ( without a i r drag ) = T r a j e c t o r y ’ )
show ( )

Output
Output 1:
Enter t h e v e l o c i t y o f p r o j e c t i o n 5
Enter t h e a n g l e o f p r o j e c t i o n i n d e g r e e s : 4 5
Time o f f l i g h t = 0.721 s
H o r i z o n t a l range = 2.549 m
Maximum h e i g h t = 0.636 m

Output 2:
Enter t he v e l o c i t y o f p r o j e c t i o n 5
Enter t he a n g l e o f p r o j e c t i o n i n d e g r e e s : 3 0
Time o f f l i g h t = 0.510 s
H o r i z o n t a l range = 2.208 m
Maximum h e i g h t = 0.318 m

61
TRAJECTORY WITH CONSIDERING AIR DRAG

Algorithm
Step 1: Start.

Step 2: Import pylab and numpy modules.

Step 3: Initialise x,y

Step 4: Initialise t and set step size, dt as 0.01

Step 5: Give velocity of projection and angle of projection.

Step 6: Find out the velocity in x direction and y direction.

Step 7: Give acceleration due to gravity.

Step 8: Set c=0.5

Step 9: Initialize the array for time, velocity and position.

Step 10: Introduce while loop with the condition expression y≥ 0

Step 11: Give the equations to

find acceleration, velocity and position inside the loop.

Step 12: Append the values of velocity, position and time inside the loop.

Step 13: Outside the loop

find maximum time, range and height.

Step 14: Print all the above data and plot the graph between x and y position.

Step 15: Show.

Step 16: Stop.

Program
from numpy import *
from pylab import *
x , y =0.0 ,0.0
t , dt =0 ,0.01
u=e v a l ( i n p u t ( ’ Enter t he v e l o c i t y o f p r o j e c t i o n ’ ) )
ang=e v a l ( i n p u t ( ’ Enter a n g l e o f p r o j e c t i l e i n d e g r e e s : ’ ) )
vx=u * c o s ( r a d i a n s ( ang ) )
vy=u * s i n ( r a d i a n s ( ang ) )
g =9.8
c =0.5
x1 , y1 , vx1 , vy1 , t 1 =[x ] , [ y ] , [ vx ] , [ vy ] , [ t ]

62
w h i l e y>= 0 :
x=x+vx * dt
y=y+vy * dt
v=s q r t ( vx **2+vy * * 2)
vx=vx=c * v * vx * dt
vy=vy =(g + c * v * vy ) * dt
t=t+dt
x1 . append ( x )
y1 . append ( y )
vx1 . append ( vx )
vy1 . append ( vy )
t 1 . append ( t )
p r i n t ( ’ Time o f f l i g h t = %8.3 f s ’ % max( t 1 ) )
p r i n t ( ’ H o r i z o n t a l Range = %8.3 f m’ % max( x1 ) )
p r i n t ( ’ Maximum Height = %8.3 f m’ % max( y1 ) )
p l o t ( x1 , y1 )
xlabel ( ’x ’)
ylabel ( ’y ’)
t i t l e ( ’ P r o j e c t i l e Motion ( with a i r drag ) = T r a j e c t o r y ’ )
show ( )

Output
Output 1:
Enter t h e v e l o c i t y o f p r o j e c t i o n 5
Enter a n g l e o f p r o j e c t i l e i n d e g r e e s : 4 5
Time o f f l i g h t = 0.610 s
H o r i z o n t a l Range = 1.387 m
Maximum Height = 0.449 m

Output 2:
Enter t he v e l o c i t y o f p r o j e c t i o n 5
Enter a n g l e o f p r o j e c t i l e i n d e g r e e s : 3 0

63
Time o f f l i g h t = 0.450 s
H o r i z o n t a l Range = 1.339 m
Maximum Height = 0.245 m

64
Experiment 10

Schr̈odinger Equation

Aim
To plot wave function and probability density of a particle in the box by solving Schr̈odinger
equation.

Principle
Particle in a box is a simple quantum mechanical problem, in which a particle is trapped in a
box with infinitely hard walls. Let the motion of a particle is restricted to x-axis between x=0
and x=L. Particle does not lose energy when it collides with infinity hard walls, so that its total
energy is [Link] potential energy, U of the particle is infinity on both sides of the box and
is zero inside the well. Therefore wave function of the particle outside the box ψ=0 for x ≤0 and
x ≥L. Within the box Schr̈odinger equation becomes,

d2 ψ 2m
+ 2 Eψ = 0 (10.1)
dx2 ℏ
solution of Eqn 10.1 is,
√ √
2mE 2mE
ψ = A sin x + B cos x (10.2)
ℏ ℏ
where, A and B are constants to be evaluated. Boundary condition is given by ψ=0 for x=0 and
x=L. At x=0, cos0=1. The second term cannot describe the particle because it does not vanish
at x=0. Hence we conclude that B=0. The sine term always yields ψ=0 at x=0.
At x=L, ψ=0 only when
√
2mE
L = nπ where n = 1, 2, 3, .. (10.3)
ℏ
Eigen values constituting the energy levels of the [Link] Eqn 10.3,

(n2 π 2 ℏ2 )
En = where n = 1, 2, 3, .. (10.4)
2mL2

65
WAVE EQUATION

The particle in a box eigenfunctions are given by,

nπx
ψn (x) = A sin( ) (10.5)
L
The constant A , thus far arbitrary, can be adjusted so that is ψn (x) normalized. The normal-
ization condition over the domain of the particle 0 ≤ x ≤ a is given by,
Z a
|ψn (x)|2 dx = 1 (10.6)
0

Substituting in Eqn 10.5 for a=L, we get,

Z L Z nπ
2 2 nπx 2 L
A sin ( )dx = A sin2 θdθ (10.7)
0 L nπ 0
Finally we can write the normalized eigenfunctions,
r
2 nπx
ψn (x) = sin( ), where n = 1, 2, 3, .. (10.8)
L L
The probability density of a wave function = ψn2 . Although ψn may be positive as well as negative,
ψn2 is always positive.

Algorithm
Step 1: Start.

Step 2: Import pylab

Step 3: Read the length of the box as L

Step 4: Initialise x as 0.0

Step 5: Read the value n

Step 6: Initialise x1 and y1

Step 7: Initialise h as 0.01

Step 8: Initialise y2
r
2 nπx
Step 9: Define f(x) as sin( )
L L
Step 10: Define function f1 (x) as f(x) * f(x)

Step 11: Begin while loop for x ≤ L

Step 12: Set y as f(x)

66
Step 13: Set x = x + h

Step 14: Append x and y to x1 and y1

Step 15: Set y1 as f1 (x)

Step 16: Append y to y2

Step 17: End while loop

Step 18: Plot(x1 , y1 )

Step 19: Plot(x1 , y2 )

Step 20: Stop.

Program
from pylab import *
def f (x ) :
r e t u r n ( s q r t (2 /L) * s i n ( ( n * p i * x ) /L ) )
def f1 (x ) :
return ( f ( x )* f ( x ) )
L=e v a l ( i n p u t ( ’ e n t e r t h e l e n g t h o f t he box ’ ) )
x =0.0
n=e v a l ( i n p u t ( ’ e n t e r t he n : ’ ) )
x1 = [ ]
y1 = [ ]
y2 = [ ]
h=0.01
w h i l e x<=L :
y=f ( x )
x=x+h
x1 . append ( x )
y1 . append ( y )
y=f 1 ( x )
y2 . append ( y )
subplot (1 ,2 ,1)
g r i d ( True )
t i t l e ( ’WAVE FUNCTION’ )
xlabel ( ’x ’)
y l a b e l ( ’ wave f u n c t i o n ’ )
p l o t ( x1 , y1 )
subplot (1 ,2 ,2)
g r i d ( True )
t i t l e ( ’PROBABILITY FUNCTION’ )
xlabel ( ’x ’)
y l a b e l ( ’ p r o b a b i l i t y functon ’ )

67
p l o t ( x1 , y2 )
show ( )

Output
Output 1:
e n t e r t h e l e n g t h o f t h e box1 . 9
enter the n :16

Output 2:
e n t e r t h e l e n g t h o f t h e box1 . 9
enter the n : 1

68
Experiment 11

Rutherford Scattering

Aim
To plot the trajectory of a particle moving in a Coulomb field and to determine the angle of
deflection as a function of impact parameter.

Principle
Rutherford’s experiment is to measure the deflection of a beam of particle of gold nuclei due to
Coulomb repulsion. The electrostatic force is given by
2eZe
F =
4πϵ0 r2
Resolving and putting r2 =x2 +y2 ,the components of acceleration are

2Ze2 x 2Ze2 y
Ax = , Ay =
4mπϵ0 r3 4mπϵ0 r3
putting

2Ze2
C= ,
4mπϵ0
d2 y
ay = 2 ,
dt
d2 x
ax = 2
dt
one get,

d2 x Cx d2 y Cy
2
= 2 , = 2
dt (x + y 2 )3/2 dt 2 (x + y 2 )3/2

69
Algorithm
Step 1: Start.

Step 2: Import pylab.

Step 3: Read the initial values of velocities and position of particles in the beam, impact
parameter and the time step dt.

Step 4: Initialize for loop using range function for the desired number. Repetitions from step
3 to step 7.

Step 5: Initialize for loop using range function for the desired number. Repetitions from step
4 to step 6.

Step 6: Calculate vx and vy using the equations: vx=vx+x[i]*c*dt/(x[i]*x[i]+y[i]*y[i])**1.5

vy=vy+y[i]*c*dt/(x[i]*x[i]+y[i]*y[i])**1.5

Step 7: Calculate x and y using equations: x[i+1]=x[i]+vx*dt

y[i+1]=y[i]+vy*dt

Step 8: Plot the path of α particle.

Step 9: Calculate cot(θ/2)A.

Step 10: Plot the relation between b and cot(θ/2).

Step 11: Show.

Step 12: Stop.

Program
from pylab import *
c =21.8274
dt =0.0001
x=z e r o s ( 1 0 0 0 0 1 , ’ f l o a t ’ )

70
y=z e r o s ( 1 0 0 0 0 1 , ’ f l o a t ’ )
b=l i n s p a c e ( 0 , 1 , 1 0 )
c o t t h e t a b y 2=z e r o s ( 1 0 , ’ f l o a t ’ )
f o r j i n range ( 1 0 ) :
x [ 0 ] , y [ 0 ] , t , vx , vy==5,b [ j ] , 0 , 1 0 , 0
f o r i i n range ( 1 0 0 0 0 ) :
vx+=x [ i ] * c * dt / ( x [ i ] * x [ i ]+y [ i ] * y [ i ] * * 1 . 5 )
vy +=y [ i ] * c * dt / ( x [ i ] * x [ i ]+y [ i ] * y [ i ] * * 1 . 5 )
x [ i +1]=x [ i ]+vx * dt
y [ i +1]=y [ i ]+vy * dt
subplot (1 ,2 ,1)
xlabel ( ’x ’)
ylabel ( ’y ’)
t i t l e ( ’ path o f alpha p a r t i c l e ’ )
plot (x , y)
c o t t h e t a b y 2 [ j ]=( x [ i ] = x [ 0 ] ) / ( y [ i ] = y [ 0 ] )
subplot (1 ,2 ,2)
x l a b e l ( ’ impact parameter b ’ )
y l a b e l ( ’ cot ( thetaby2 ) ’ )
t i t l e ( ’ r e l a t i o n between b and t h e t a ’ )
g r i d ( ’ on ’ )
plot (b , cotthetaby2 )
show ( )

Output

71
Experiment 12

Harmonic Oscillator

Aim
To find the solution of second-order differential equation using Runge-Kutta method.

Principle
The fourth order RK method can be used to solve numerically the higher order ordinary differ-
ential equations. Let us consider a second order differential equation of the form.

d2 y
= g(x, y, z) (12.1)
dx2
put
dy
=z (12.2)
dx
Then , the Eqn 12.1 has been reduced to first order simultaneous differential equation.
dy
= z = f (x, y, z) (12.3)
dx
dz
= g(x, y, z) (12.4)
dx
Now we can directly write down the RK fourth order formula for deriving the system (Eqn 12.3).
Let the initial connection of above system be given by

y(xn ) = yn (12.5)
y′(xn ) = z(xn ) (12.6)
k1 = hf (tn , yn , zn ) (12.7)
m1 = hg(tn , yn , zn ) (12.8)

72
 h k1 m1
k2 = hf (tn + , yn + , zn + ) (12.9)
2 2 2
h k1 m1
m2 = hg(tn + , yn + , zn + ) (12.10)
2 2 2
h k2 m2
k3 = hf (tn + , yn + , zn + ) (12.11)
2 2 2
h k2 m2
m3 = hg(tn + , yn + , zn + ) (12.12)
2 2 2
k4 = hf (tn + h, yn + k3 , zn + m3 ) (12.13)
m4 = hg(tn + h, yn + k3 , zn + m3 ) (12.14)

Now using initial conditions yn , zn and fourth order RK fourmula we compute,

1
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) (12.15)
6
1
zn+1 = zn + (m1 + 2m2 + 2m3 + m4 ) (12.16)
6
This method can be extended on similar tends to solve system of n first order differential equa-
tions. The differential of a harmonic oscillator with or without damping is

d2 y dy
2
= 2k0 − ω2y (12.17)
dx dx
Here we put

dy
= γ or z (12.18)
dx
f (t, y, v) = γ or z (12.19)
g(t, y, v) = - 2k0 γ - ω 2 y (12.20)

k0 is the damping coefficient and if k0 =0, then it represent the undamped case of harmonic
oscillator.

Algorithm
Step 1: Start.

Step 2: Import math module

Step 3: Import pylab

Step 4: Define the differential equation

Step 5: Define frequency (ω) , final value of x (xf) and stepsize h

Step 6: Set of initial conditions x,y,z

Step 7: Initialize lists x1, y1, z1 to store the values of x, y, z

73
 Step 8: Input the value of damping constant , k0

Step 9: Do steps 10 to steps 23 while x ≤ xf

Step 10: k1 ← hf(tn , yn , zn )

Step 11: m1 ← hg(tn , yn , zn )

Step 12: k2 ← hf(tn + h2 , yn + k1

2
, zn + m1
2
)

Step 13: m2 ← hg(tn + h2 , yn + k1

2
, zn + m1
2
)

Step 14: k3 ← hf(tn + h2 , yn + k2

2
, zn + m2
2
)

Step 15: m3 ← hg(tn + h2 , yn + k2

2
, zn + m2
2
)

Step 16: k4 ← hf(tn +h, yn +k3 , zn +m3 )

Step 17: m4 ← hg(tn +h, yn +k3 , zn +m3 )

Step 18: x ← x + h
1
Step 19: y ← y + 6
(k1 + 2k2 + 2k3 + k4 )
1
Step 20: z ← z + 6
(m1 + 2m2 + 2m3 + m4 )

Step 21: Append the value of x in the list x

Step 22: Append the value of y in the list y

Step 23: Append the value of z in the list z

Step 24: Print the value of frequency ω and stepsize h

Step 25: Plot the graph between time and position

Step 26: Plot the graph between time and velocity

Step 27: Plot the graph between velocity and position

Step 28: Show.

Step 29: Stop.

Program
from math import *
from pylab import *
def f (x , y , z ) : return z
d e f g ( x , y , z ) : r e t u r n ( = 2 * k0 * z=w** 2 * y )
w=20.2
x f=5
h=0.005

74
x =0.0
y =0.0
z =0.4
x1 , y1 , z1 = [ ] , [ ] , [ ]
k0=e v a l ( i n p u t ( ’ e n t e r th e v a l u e s o f k ’ ) )
w h i l e x<=x f :
k1=h * f ( x , y , z )
m1=h * g ( x , y , z )
k2=h * f ( x+(h / 2 ) , y+(k1 / 2 ) , z+(m1/ 2 ) )
m2=h * g ( x+(h / 2 ) , y+(k1 / 2 ) , z+(m1/ 2 ) )
k3=h * f ( x+(h / 2 ) , y+(k2 / 2 ) , z+(m2/ 2 ) )
m3=h * g ( x+(h / 2 ) , y+(k2 / 2 ) , z+(m2/ 2 ) )
k4=h * f ( x+h , y+k3 , z+m3)
m4=h * g ( x+h , y+k3 , z+m3)
x , y , z=x+h , y+(k1+2* k2+2* k3+k4 ) / 6 , z+(m1+2*m2+2*m3+m4)/ 6
x1 . append ( x )
y1 . append ( y )
z1 . append ( z )
p r i n t ( ’ v a l u e o f w=%6.4 f , h=%6.4 f ’%(w, h ) )
figure (1)
subplot (2 ,1 ,1)
p l o t ( x1 , y1 )
x l a b e l ( ’ time ==>’)
y l a b e l ( ’ p o s i t i o n ===>’)
title ( ’ oscillator ’)
subplot (2 ,1 ,2)
p l o t ( x1 , z1 )
x l a b e l ( ’ time ===>’)
y l a b e l ( ’ v e l o c i t y ===>’)
figure (2)
p l o t ( y1 , z1 )
x l a b e l ( ’ p o s i t i o n ===>’)
y l a b e l ( ’ v e l o c i t y ===>’)
t i t l e ( ’ o s c i l l a t o r phase s p a c e p l o t ’ )
show ( )

Output
Output 1:
e n t e r t h e v a l u e s o f k0
v a l u e o f w=20.2000 , h =0.0050
Output 2:
e n t e r t h e v a l u e s o f k0 . 8
v a l u e o f w=20.2000 , h =0.0050

75
Figure 12.1: Output 1

Figure 12.2: Output 2

76