EP4210/24040 - Computational Physics
Kirit Makwana
Department of Physics, IIT Hyderabad
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 1 / 27
Numerical Integration
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 2 / 27
Motivation
Several examples of integration problems in Physics
Many quantities are defined in the form of density and require
integration to obtain global values
∫
M= ρdV (1)
∫V
W = F · dl (2)
∫ L
ϕ = B · dA (3)
S
Integration also comes up in averages
∫
⟨f(x)⟩ = f(x)|ψ(x)|2 dx (4)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 3 / 27
1d integration
First fundamental theorem of calculus
∫ x=b
f(x)dx = F(b) − F(a) (5)
x=a
where F(x) is the anti-derivative of function f(x)
Riemann sum
∫ x=b ∑
N−1
f(x)dx = lim f(xi )h (6)
x=a h→0
i=0
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 4 / 27
First approximation
The Riemann sum itself provides a way to approximate the integral
which is the area under the curve
∫ x=b ∑
N−1
f(x)dx ≈ ∆xf(xi ) = Ĩ (7)
x=a i=0
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 5 / 27
Error
Lets calculate integral of f(x) = x2 with a = 0
∑
N−1
Ĩ = ∆x (i∆x)2 (8)
i=0
∑
N−1
3 2
= (∆x) i (9)
i=0
Now consider
(i − 1)3 = i3 − 3i2 + 3i − 1 (10)
3 3 2
=⇒ i − (i − 1) = 3i − 3i + 1 (11)
Now sum over
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 6 / 27
Error
This gives
∑
N ∑
N ∑
N
3 3 2
i − (i − 1) = 3 i −3 i+N+1 (12)
i=0 i=0 i=0
∑
N
=⇒ 1 + N3 = 3 i2 − 3(N)(N + 1)/2 + N + 1 (13)
i=0
From this we get
∑
N
N3 N2 N
i2 = + + (14)
3 2 6
i=0
This gives
( )
(N − 1)3 (N − 1)2 N − 1
Ĩ = ∆x3 + + (15)
3 2 6
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 7 / 27
The integral is
( )
N3 N2 N
3
Ĩ = ∆x − + (16)
3 2 6
b 3 b2 b
= − ∆x + ∆x2 (17)
3 2 6
The exact solution of the integral is
b3
I= (18)
3
From this we can see that the error in the approximation is O(∆x)
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 8 / 27
Trapezoidal approximation
Instead of taking a rectangle in the region, take a trapezoid at the top
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 9 / 27
Trapezoidal integration
The approximation is
∑
N−1
fi+1 − fi ∑ fi + fi+1 N−1
Ĩ = ∆x (fi + ) = ∆x (19)
2 2
i=0 i=0
Take the example of fi = i2 ∆x2
∑(
N−1
1
)
3 2
Ĩ = ∆x i +i+ (20)
2
i=0
( )
(N − 1)3 (N − 1)2 N − 1 N(N − 1) N
= ∆x3 + + + + (21)
3 2 6 2 2
( 3 )
N N
= ∆x3 + (22)
3 6
b3 b
= + ∆x2 (23)
6 6
The error is O(∆x2 )
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 10 / 27
Algorithm
Exercise: Evaluate
∫ 1
sin(x)
I= dx (24)
0 x
Use the Trapezoidal integration method
#d e f i n e t h e f u n c t i o n and v e c t o r i z e i t
v e c f u n c = np . v e c t o r i z e ( f u n c )
Accurate answer: 0.946083070367183
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 11 / 27
Interpolation
Lets interpolate the function between the grid points to get a better
approximation
Lets build polynomials in some smaller set of grid points, and then
use these polynomials for integration
Sub-divide the N + 1 grid points into groups of n + 1 points
If there are N + 1 grid points, choose a number n + 1 < N such that
N is divisible by n
Then we can have a series of sets of n + 1 points each
How many sets of points do we have?
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 12 / 27
Interpolation
The number of sets m will be
N
m= (25)
n
Now in each of these m sets, there are n + 1 points
Lets label these points j = 0 to j = n, so that xj = xk + j∆x
Here xk will be k ∗ ∆x ∗ n for k = 0 to k = m − 1
Example, take N = 20 and n = 4
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 13 / 27
Polynomial interpolation
Consider a set of n + 1 points xj , j = 0 to j = n
We know the value of the function at each of these points, i.e.
fj = f(xj )
We can find a polynomial of power n such that it passes through each
of the points fj
This is because a n power polynomial has n + 1 free coefficients
In fact we can construct such a polynomial explicitly
∑
n
Pn (x) = fj Lj (x) (26)
j=0
= f0 L0 (x) + f1 L1 (x) + ... (27)
This is called Lagrange’s interpolation formula
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 14 / 27
Lagrange interpolation
Here Lj (x) is an nth power polynomial given by
∏
n
x − xk
Lj (x) = (28)
xj − xk
k=0,k̸=j
x − x0 x − xj−1 x − xj+1
= ... (29)
xj − x0 xj − xj−1 xj − xj+1
Check that this is power n polynomial
It is not singular
It satisfies Pn (xj ) = fj
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 15 / 27
Polynomial Integration
Now we can integrate over this polynomial
∫ b
I= f(x)dx (30)
a
∫ b
Ĩ = Pn (x)dx (31)
a
Define x = a + t∆x, =⇒ dx = dt∆x
∫ n
Ĩ = ∆x Pn (t)dt (32)
0
∑
n
Pn (t) = fj Lj (t) (33)
j=0
∏
n
(a + t∆x) − (a + k∆x)
Lj (t) = (34)
(a + j∆x − (a + k∆x))
k=0,k̸=j
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 16 / 27
Integration
This reduces to
∏ t−k
Lj (t) = (35)
j−k
k=0,k̸=j
Now Ĩ become
∫ n ∑
n
Ĩ = ∆x dt fj Lj (t) (36)
0 j=0
Consider the integral of just one of these terms
∫ n ∫ n ∏ t−k
dtLj (t) = dt (37)
0 0 j−k
k=0,k̸=j
Now these integrals can be evaluated, because they are simply
polynomials, with no information about the original function f
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 17 / 27
Weights of integral
These define what are known as the weights of the integral
∫ n ∫ n ∏ t−k
wj = dtLj (t) = dt (38)
0 0 j−k
k=0,k̸=j
These weights only depend on the choice of n
Once these are determined, then the approximate intergral can be
evaluated as
∑
n
Ĩ = ∆x fj wj (39)
j=0
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 18 / 27
Trapezoidal rule
If we take n = 1, that means just taking two neighbouring points,
then a = xj and b = xj+1
∫ 1
1
w0 = dt(1 − t) = (40)
0 2
∫ 1
1
w1 = dtt = (41)
0 2
Using this we get
∑
n
Ĩ = ∆x fj wj (42)
j=0
fj + fj+1
= ∆x (43)
2
This is simply the trapezoidal rule
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 19 / 27
Simpson’s rule
Consider n = 2. Then a = xj and b = xj+2
Lets calculate the weights
∫ 2
1 1
w0 = (1 − t)(2 − t)dt = (44)
0 2 3
∫ 2
4
w1 = t(2 − t)dt = (45)
0 3
∫ 2
t(t − 1) 1
w2 = dt = (46)
0 2 3
Using these weights we get
∆x
Ĩ = (fj + 4fj+1 + fj+2 ) (47)
3
This is called as Simpson’s rule
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 20 / 27
Larger interval
Suppose there is an interval [a, b] divided into N segments, so j = 0
to j = N
Then the integral with n = 1 is (Trapezoidal rule)
∑
j=N−1
fj + fj+1
Ĩ = ∆x (48)
2
j=0
Similarly, if we take n = 2 (Simpsons’s rule), in which case N has to
be even
∑
j=N−2
∆x
Ĩ = (fj + 4fj+1 + fj+2 ) (49)
3
j=0,2,4
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 21 / 27
Algorithm
#d e f i n e N wh i ch s h o u l d be e v e n
d e l t a x = ( b−a ) /N
x a r r = np . l i n s p a c e ( a , b , N)
#d e f i n e f u n c
v e c t o r _ f u n c = np . v e c t o r i z e ( f u n c )
func_arr = vector_func ( xarr )
i n t e g r a l =0.0
f o r j i n range ( 0 ,N− 1 , 2 ) :
^^ I i n t e g r a l = i n t e g r a l +( d e l t a x / 3 . 0 ) ∗ ( f u n c _ a r r [ j ] \
^^ I ^^ I +4∗ f u n c _ a r r [ j +1]+ f u n c _ a r r [ j +2])
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 22 / 27
Accuracy
Consider the Taylor series expansion of the function around point j + 1
f(x) = f(xj ) + (x − xj )f′ (xj ) + (x − xj )2 f′′ (xj ) + O(∆x)3 (50)
If we ignore the higher order term, this is actually a second order
polynomial in x
If f(x) is a second order polynomial, then the n = 2 approximation will
give us the exact solution of the problem
From this we can argue that the error of the Simpson method is
O(∆x)3
Actually it can be showed that the error is even better or higher order
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 23 / 27
2D integration
Consider an integration in 2 dimensions
∫ bx [ ∫ by ]
I= f(x, y)dy dx (51)
ax ay
Lets say the grid is ∆x = (bx − ax )/Nx and ∆y = (by − ay )/Ny
Define the first integral as
∫ by
Iy (xi ) = f(xi , yj )dyj (52)
ay
This can be evaluated with the Simpson or Trapezoidal rule, with xi
just being a parameter. Then again Iy (xi ) can be integrated with
Simpson rule
∫ bx
I= Iy (xi )dx (53)
ax
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 24 / 27
2D integration
If the limit depends on the other variable, for, ex, if the integration
domain is a unit circle x2 + y2 = 1
√ √
Then ay = − 1 − x2 and by = 1 − x2
∫ √ 2 1−xi
Iy (xi ) = √ f(xi , y)dy (54)
− 1−x2i
This integral can be done by Simpson method, depending on how the
data is provided, interpolation and grid spacing can be adjusted
And then this is followed by
∫ 1
I= Iy (xi )dx (55)
−1
which can again be done with Simpson method
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 25 / 27
3D integration
We have volume integral
∫ bz {∫ by [∫ bx ]}
I= dz dy dxf(x, y, z) (56)
az ay ax
This can be treated as 3 successive 1d integrals, or as a 2D integral
and a 1D integral
{ ∫ by [ ∫ bx ]}
Ĩx,y (z) = dy dxf(x, y, z) (57)
ay ax
∫ bz
Ĩ = dzĨx,y (z) (58)
az
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 26 / 27
Exercise
A thin circular disk of radius unity has a surface density of
σ = exp(−x2 − y2 ) (59)
Calculate its mass by implementing a 2D numerical integration with
Simpson’s rule
Compare with the analytical result
Kirit MakwanaDepartment of Physics, IIT Hyderabad
EP4210/24040 - Computational Physics 27 / 27