Indian Institute of Technology Gandhinagar

PH 509: Computational Physics

Assignment - 2

Due date: 17th Aug

Name: [SHUBHAM GARG] Roll number: [16510081]

Instructions :
Submit Assignment based on this format
Your File Name for submission should be [N ame] [RollN o] Assignment 2.pdf
Please copy the problem, write your approach, your python code , results (Figures, tables) etc.

1 Problem 2
You are given a function f (x) = x3 + 2x2 + 3x + 1. Plot this function using matplotlib when x [5, 5]
with 512 points. Get the input from the user to give a number x0 between [-2,2]. At the point (x0 , f (x0 ),
draw a tangent to this curve. Label the axes and curve. Save the output as a png file.

1.1 Approach
Plotting a function with x can be easily done by using plot function of pyplot library. For drawing
tangent at some point of this curve, we need slope of curve at that point. for this we import sympy
library and use its diff and subs function. Using this slope and the initial condition , that is our x0 ,
we can make an equation for straight line and plot the same.

1.2 Program
Python script for the above problem is written as:
import m a t p l o t l i b . p y p l o t a s p l t
import numpy a s np
from sympy import

x=symbols ( x )
y=x3+2x2+3x+1
m= d i f f ( y , x )
X=input ( X )
s l o p e=m. s u b s ( x ,X)
c=y . s u b s ( x ,X)
print m, s l o p e , c
x=np . l i n s p a c e ( 5 , 5 , 5 1 2 )
t=s l o p e ( xX)+c
y=x3+2x2+3x+1
p l t . p l o t ( x , t , r -- , l a b e l= tangent )
p l t . show ( )
p l t . p l o t ( x , y , b , l a b e l=f(x) )
p l t . x l a b e l ( x )
p l t . y l a b e l ( y )
plt . legend ()
p l t . g r i d ( True )

1.3 Results

1
Figure 1: F(x) and its tangent at x=1.75

Figure 2: F(x) and its tangent at x=-1.75

2
2 Problem 3
Pn sin(2(2k1)x)
Suppose f (x, n) = k=1 2k1 . The use will give the value of n. Plot f (x, n) for that n. See the
trend as n increases.

2.1 Approach
First we create 1D array for x. Now for every value of x in this array, we need to evaluate f(x,n) where
n is given by user. For summation we use for loop and run it from 0 to n. Then plot f(x,n) vs x.

2.2 Program
Python script for above problem.

import m a t p l o t l i b . p y p l o t a s p l t
import numpy a s nm
n=input ( n= )
x=nm . a r a n g e ( 0 , 1 0 0 0 , 2 )
def f ( x , n ) :
sum=0
f o r k in range ( 1 , n +1):
y=(nm . s i n ( 2 nm . p i ( 2 k 1.0) x ) ) / ( 2 k 1.0)
sum=sum+y
return sum
plt . close ()
p l t . p l o t ( x , f ( x , n ) , y -- )
p l t . x l a b e l ( X )
p l t . y l a b e l ( Y )
p l t . t i t l e ( f(x ,1000) vs x )
p l t . show ( )

3
2.3 Results
We plotted the above program for three different values of n.

Figure 3: for n=10

Figure 4: for n=100

4
Figure 5: for n=1000

3 Problem 4
Problem 4 (From: Computational Physics, Nicolas Giordano- Ch 1, Problem 6)

3.1 Approach
We consider N to be the population at any time t. So our function is the rate of change of population with
time and is first order differential equation. First we use euler method to solve the diffrential equation.
Given initial polpulation, we can solve for the population at time t0 + timestep. So our accuracy depends
on the timestep we are taking. Programming part is simple. We make an array for time,t. We already
have initial population at t = 0. Running a for loop, we calculate population by using previous value from
the list. Thus we get population array corresponding to time array. Using plot, we get time evolution
of population. We then compare this result from euler method by solving our function analytically and
plotting the same.

3.2 Program
Python script

import m a t p l o t l i b . p y p l o t a s p l t
import numpy a s np
from sympy import
plt . close ()
def f ( x ) :
y=a xbx 2
return y

#i n i t i a l c o n d i t i o n s
t 0=0
t f =10
h=0.1
a=1
b=0

5
c =1000 #i n i t i a l p o p u l a t i o n
N=[ c ]
#a r r a y f o r time
t=np . a r a n g e ( t0 , t f , h )

#E u l e r method

f o r i in range ( len ( t ) ) :
y=N[ 1]+h f (N[ 1 ] )
i f len (N)<len ( t ) :
N. append ( y )
N=np . a r r a y (N)
p l t . p l o t ( t , N, l a b e l= Euler method )
p l t . x l a b e l ( time )
p l t . y l a b e l ( Population )
p l t . t i t l e ( Population vs Time )
p l t . g r i d ( True )
p l t . show ( )

#a n a l y t i c a l s o l u t i o n
M= [ ]
f o r i in t :
y=c np . exp ( a i )
M. append ( y )
M=np . a r r a y (M)
p l t . p l o t ( t ,M, r , l a b e l= Analytic solution )
p l t . l e g e n d ( t i t l e = time step =0.1 )
p l t . show ( )

#A n a l y t i c s o l u t i o n f o r ( 1 0 , 3 )
import m a t p l o t l i b . p y p l o t a s p l t
import numpy a s np
t=np . a r a n g e ( 0 , 1 5 , 0 . 0 1 )
plt . close ()
G= [ ]
f o r i in t :
y=300np . exp ( i / 1 0 ) / ( 3 0 0 np . exp ( i /10) 299)
G. append ( y )
G=np . a r r a y (G)

plt . p l o t ( t , G, y )
plt . x l a b e l ( Time )
plt . y l a b e l ( Population )
plt . g r i d ( True )
plt . t i t l e ( Population vs Time )
plt . show ( )

3.3 Results
To show the accuracy in eulers method due to change in time step, we plotted the euler method for two
different time steps.

6
Figure 6: a=1,b=0,Euler method

Figure 7: a=1,b=0,Analytic method

7
Figure 8: a=10,b=3,Analytic method