Computational Physics Lab Solutions

Satyam Kumar Singh

May 11, 2026

Problem 1: Bisection Method

Question: Use the bisection method to find the positive root of x3 − 4x − 9 = 0 correct
up to 10−5 .
Explanation: Let f (x) = x3 − 4x − 9. We first find an interval [a, b] such that
f (a)f (b) < 0. Testing values: f (2) = 8 − 8 − 9 = −9 and f (3) = 27 − 12 − 9 = 6.
Since the signs differ, a positive root exists in [2, 3]. The bisection method iteratively halves
this interval using xm = a+b 2
, updating the boundaries based on the sign of f (xm ) until
−5
|b − a|/2 ≤ 10 .
Python Code:

1 def f ( x ) :
2 return x **3 - 4* x - 9
3

4 a , b = 2.0 , 3.0
5 tol = 1e -5
6

7 if f ( a ) * f ( b ) > 0:
8 print ( " Root is not bracketed . " )
9 else :
10 while ( b - a ) / 2.0 > tol :
11 mid = ( a + b ) / 2.0
12 if f ( mid ) == 0:
13 break
14 elif f ( a ) * f ( mid ) < 0:
15 b = mid
16 else :
17 a = mid
18 root = ( a + b ) / 2.0
19 print ( f " Root : { root :.5 f } " )

1
 Problem 2: Newton-Raphson Method
Question: Apply the Newton-Raphson method to determine the root of e−x − x = 0
starting from x0 = 0.5.
Explanation: Let f (x) = e−x − x. The first derivative is f ′ (x) = −e−x − 1. The
Newton-Raphson update formula is:

f (xn )
xn+1 = xn −
f ′ (xn )

We iterate this formula starting from x0 = 0.5 until the difference between successive ap-
proximations is negligible.
Python Code:

1 import math
2

3 def f ( x ) :
4 return math . exp ( - x ) - x
5

6 def df ( x ) :
7 return - math . exp ( - x ) - 1
8

9 x0 = 0.5
10 tol = 1e -6
11

12 while True :
13 x1 = x0 - f ( x0 ) / df ( x0 )
14 if abs ( x1 - x0 ) < tol :
15 break
16 x0 = x1
17

18 print ( f " Root : { x1 :.5 f } " )

Problem 3: Secant Method

Question: Use the secant method to solve cos x − x = 0 with initial guesses x0 = 0 and
x1 = 1.
Explanation: Let f (x) = cos x − x. The Secant method approximates the derivative
using two previous points, resulting in the update formula:
xn − xn−1
xn+1 = xn − f (xn )
f (xn ) − f (xn−1 )

2
 We use the initial guesses x0 = 0 and x1 = 1 and iterate until convergence.
Python Code:

1 import math
2

3 def f ( x ) :
4 return math . cos ( x ) - x
5

6 x0 = 0.0
7 x1 = 1.0
8 tol = 1e -6
9

10 while abs ( x1 - x0 ) > tol :

11 x2 = x1 - f ( x1 ) * ( x1 - x0 ) / ( f ( x1 ) - f ( x0 ) )
12 x0 = x1
13 x1 = x2
14

15 print ( f " Root : { x1 :.5 f } " )

Problem 4: Quantum Well Energy Levels

p
Question: Solve k tan(ka) = k02 − k 2 for a = 1 and k0 = 5 to determine the lowest
allowed value of k numerically. √
Explanation: Rearranging the equation yields f (k) = k tan(k) − 25 − k 2 = 0. The
tangent function has an asymptote at k = π/2 ≈ 1.57. The lowest allowed k must lie in the
interval (0, 1.57). We can use the bisection method on an interval like [1.0, 1.5] to locate this
root.
Python Code:

1 import math
2

3 def f ( k ) :
4 return k * math . tan ( k ) - math . sqrt (25 - k **2)
5

6 a , b = 1.0 , 1.5
7 tol = 1e -5
8

9 while ( b - a ) / 2.0 > tol :

10 mid = ( a + b ) / 2.0
11 if f ( a ) * f ( mid ) < 0:
12 b = mid
13 else :
14 a = mid

3
15

16 k_lowest = ( a + b ) / 2.0
17 print ( f " Lowest allowed k : { k_lowest :.5 f } " )

Problem 5: Gaussian Elimination

Question: Solve the system 2x + 2y − z = 8, −3x − y + 2z = −11, −2x + y + 2z = −3
using Gaussian elimination.
Explanation: We represent the system as an augmented matrix [A|B]. Gaussian
elimination involves two main steps: 1. Forward Elimination: Applying row operations
to convert the matrix into an upper triangular form. 2. Back Substitution: Solving for
the variables starting from the bottom equation (which now has only one variable) up to the
top.
Python Code:

1 import numpy as np
2

3 A = np . array ([[ 2.0 , 2.0 , -1.0] ,

4 [ -3.0 , -1.0 , 2.0] ,
5 [ -2.0 , 1.0 , 2.0]])
6 B = np . array ([8.0 , -11.0 , -3.0])
7 n = len ( B )
8

9 # Forward Elimination
10 for i in range ( n ) :
11 for j in range ( i +1 , n ) :
12 factor = A [ j ][ i ] / A [ i ][ i ]
13 for k in range (i , n ) :
14 A [ j ][ k ] -= factor * A [ i ][ k ]
15 B [ j ] -= factor * B [ i ]
16

17 # Back Substitution
18 x = np . zeros ( n )
19 x [n -1] = B [n -1] / A [n -1][ n -1]
20

21 for i in range (n -2 , -1 , -1) :

22 sum_ax = 0
23 for j in range ( i +1 , n ) :
24 sum_ax += A [ i ][ j ] * x [ j ]
25 x [ i ] = ( B [ i ] - sum_ax ) / A [ i ][ i ]
26

27 print ( f " Solution : x ={ x [0]:.2 f } , y ={ x [1]:.2 f } , z ={ x [2]:.2 f } " )

4
 Problem 6: Gauss-Seidel Method
Question: Solve the system 10x − y + 2z = 6, −x + 11y − z = 25, 2x − y + 10z = −11
using the Gauss-Seidel method. Perform four iterations starting from (0, 0, 0).
Explanation: The Gauss-Seidel method is an iterative technique. First, we rewrite the
system by isolating x, y, and z on the left side:
1
x= (6 + y − 2z)
10
1
y = (25 + x + z)
11
1
z = (−11 − 2x + y)
10
Starting with initial guesses x0 = 0, y0 = 0, z0 = 0, we update each variable sequentially,
using the most recently calculated values of the other variables immediately in the next step.
Python Code:

1 x , y , z = 0.0 , 0.0 , 0.0

2 iterations = 4
3

4 for i in range ( iterations ) :

5 x = (6.0 + y - 2.0* z ) / 10.0
6 y = (25.0 + x + z ) / 11.0
7 z = ( -11.0 - 2.0* x + y ) / 10.0
8 print ( f " Iter { i +1}: x ={ x :.4 f } , y ={ y :.4 f } , z ={ z :.4 f } " )

Problem 7: Linear Interpolation

Question: The temperature distribution in a rod is given at x = {0, 1, 2, 3, 4} as T =
{100, 80, 65, 50, 40}. Estimate the temperature at x = 2.5 m.
Explanation: The target point x = 2.5 lies between the data points (x1 , T1 ) = (2, 65)
and (x2 , T2 ) = (3, 50). The linear interpolation formula uses the equation of a straight line
connecting these two adjacent points:
T2 − T1
T (x) = T1 + (x − x1 )
x2 − x 1
Python Code:

1 x1 , T1 = 2.0 , 65.0
2 x2 , T2 = 3.0 , 50.0
3 x_target = 2.5

5
4

5 T_target = T1 + (( T2 - T1 ) / ( x2 - x1 ) ) * ( x_target - x1 )
6 print ( f " Estimated Temperature at x =2.5: { T_target :.1 f } C " )

Problem 8: Lagrange Interpolation

Question: Estimate the value of f (1.5) from the data: (1, 2), (2, 3), (3, 5).
Explanation: The Lagrange interpolating polynomial for n points is given by P (x) =
Pn−1
i=0 yi Li (x), where the basis polynomials Li (x) are:
n−1
Y x − xj
Li (x) =
j=0
xi − x j
j̸=i

For our three points, we calculate L0 (1.5), L1 (1.5), and L2 (1.5) and sum their products with
the corresponding y values.
Python Code:

1 x_data = [1.0 , 2.0 , 3.0]

2 y_data = [2.0 , 3.0 , 5.0]
3 x_target = 1.5
4

5 result = 0.0
6 n = len ( x_data )
7

8 for i in range ( n ) :
9 term = y_data [ i ]
10 for j in range ( n ) :
11 if j != i :
12 term *= ( x_target - x_data [ j ]) / ( x_data [ i ] - x_data [ j ])
13 result += term
14

15 print ( f " f (1.5) approx = { result :.3 f } " )

Problem 9: Finite Differences

Question: A particle trajectory is given by t = {0, 1, 2, 3, 4} and x = {0, 2, 8, 18, 32}.
Estimate velocity and acceleration at t = 2 s.
Explanation: We use the central finite difference approximations for the first and
second derivatives. With a uniform step size of h = 1 s:

6
 • Velocity (1st derivative): v(t) ≈ x(t+h)−x(t−h)
2h

• Acceleration (2nd derivative): a(t) ≈ x(t+h)−2x(t)+x(t−h)

h2

At t = 2, we use x(3) = 18, x(2) = 8, and x(1) = 2.

Python Code:

1 t = [0 , 1 , 2 , 3 , 4]
2 x = [0 , 2 , 8 , 18 , 32]
3 h = 1.0
4 index = 2 # t =2 is at index 2
5

6 velocity = ( x [ index +1] - x [ index -1]) / (2 * h )

7 acceleration = ( x [ index +1] - 2* x [ index ] + x [ index -1]) / ( h **2)
8

9 print ( f " Velocity at t =2: { velocity :.2 f } m / s " )

10 print ( f " Acceleration at t =2: { acceleration :.2 f } m / s ^2 " )

Problem 10: Trapezoidal Rule

R1 dx
Question: Evaluate 0 1+x2
using n = 4 intervals. Discuss the accuracy and how to decide
n for fourth decimal place accuracy.
Explanation: The trapezoidal rule approximates the integral as:

" n−1
#
h X
I≈ f (x0 ) + 2 f (xi ) + f (xn )
2 i=1

For n = 4, h = (1 − 0)/4 = 0.25. To determine the required n for 10−4 accuracy, we use the
3
error bound |E| ≤ (b−a)
12n2
max |f ′′ (x)|. Since max |f ′′ (x)| for 1+x
1
2 on [0, 1] is 2 (at x = 0), we
2 −4
solve 12n2 ≤ 10 , which yields n ≥ 41.
Python Code:

1 def f ( x ) :
2 return 1.0 / (1.0 + x **2)
3

4 a , b = 0.0 , 1.0
5 n = 4
6 h = (b - a) / n
7

8 integral = f ( a ) + f(b)
9 for i in range (1 , n):
10 x_i = a + i * h
11 integral += 2 * f ( x_i )

7
12

13 integral *= h / 2.0
14

15 print ( f " Integral with n =4: { integral :.4 f } " )

Problem 11: Simpson’s 1/3 Rule

Rπ
Question: Evaluate 0
sin x dx using Simpson’s 1/3 rule with n = 6 intervals. Compare
the result with exact value of integral.
π−0
Explanation: The step size is h = 6
= π6 . The points are xi = i π6 for i = 0 to 6.
Simpson’s 1/3 rule approximates the integral as:
" #
h X X
I≈ f (x0 ) + 4 f (xi ) + 2 f (xi ) + f (xn )
3 i=1,3,5 i=2,4

Rπ
The exact analytical value is 0
sin x dx = [− cos x]π0 = −(−1) − (−1) = 2.
Python Code:

1 import math
2

3 def f ( x ) :
4 return math . sin ( x )
5

6 a , b = 0.0 , math . pi
7 n = 6
8 h = (b - a) / n
9

10 integral = f ( a ) + f ( b )
11 for i in range (1 , n ) :
12 x_i = a + i * h
13 if i % 2 == 0:
14 integral += 2 * f ( x_i )
15 else :
16 integral += 4 * f ( x_i )
17

18 integral *= h / 3.0
19 exact = 2.0
20 error = abs ( exact - integral )
21

22 print ( f " Simpson ’s Integral : { integral :.6 f } " )

23 print ( f " Exact Value : { exact :.6 f } " )
24 print ( f " Error : { error :.6 e } " )

8
 Problem 12: Normalization Integral of a Gaussian
R∞ 2
Question: Compute the normalization integral I = −∞
e−x dx using the trapezoidal
method over the interval [−4, 4]. Choose suitable h to get the result accurate up to third
decimal place.
Explanation: We approximate the infinite limits with [−4, 4] since e−16 is negligible.
To guarantee an accuracy of 5 × 10−4 (third decimal place), we use the error bound for
2 2
the trapezoidal rule: |E| ≤ (b−a)h
12
max |f ′′ (x)|. For f (x) = e−x , the second derivative is
2
f ′′ (x) = e−x (4x2 − 2), which has a maximum absolute value of 2 at x = 0.

(4 − (−4))h2 16h2 √
· 2 ≤ 0.0005 =⇒ ≤ 0.0005 =⇒ h ≤ 0.000375 ≈ 0.019
12 12
Choosing h = 0.01
√ (or n = 800) will safely achieve the desired accuracy. The exact value of
the integral is π ≈ 1.772.
Python Code:

1 import math
2

3 def f ( x ) :
4 return math . exp ( - x **2)
5

6 a , b = -4.0 , 4.0
7 h = 0.01 # Chosen to ensure 3 rd decimal accuracy
8 n = int (( b - a ) / h )
9

10 integral = f ( a ) + f(b)
11 for i in range (1 , n):
12 x_i = a + i * h
13 integral += 2 * f ( x_i )
14

15 integral *= h / 2.0
16 exact = math . sqrt ( math . pi )
17

18 print ( f " Calculated Integral : { integral :.4 f } " )

19 print ( f " Exact Value ( sqrt ( pi ) ) : { exact :.4 f } " )

Problem 13: Euler’s Method

dy
Question: Solve the ordinary differential equation dx
= x + y, with y(0) = 1 using Euler’s

9
 method up to x = 0.5 with step size h = 0.1.
Explanation: Euler’s method estimates the next value using the tangent at the current
point:
yn+1 = yn + h · f (xn , yn )
Given f (x, y) = x + y, x0 = 0, y0 = 1, h = 0.1, we perform 5 iterations to reach x = 0.5.
Python Code:

1 def f (x , y ) :
2 return x + y
3

4 x0 , y0 = 0.0 , 1.0
5 x_end = 0.5
6 h = 0.1
7

8 x , y = x0 , y0
9 n_steps = int (( x_end - x0 ) / h )
10

11 print ( f " x_0 ={ x :.1 f } , y_0 ={ y :.4 f } " )

12 for i in range ( n_steps ) :
13 y = y + h * f (x , y )
14 x = x + h
15 print ( f " x_ { i +1}={ x :.1 f } , y_ { i +1}={ y :.4 f } " )

Problem 14: Fourth-Order Runge-Kutta (RK4) Method

dy
Question: Solve dt
= −2y, with y(0) = 1 using the fourth-order Runge-Kutta method for
one step of size h = 0.2.
dy
Explanation: For a given ODE dt
= f (t, y), the RK4 method uses a weighted average
of four slopes to advance the solution:

k1 = hf (tn , yn )
 
h k1
k2 = hf tn + , yn +
2 2
 
h k2
k3 = hf tn + , yn +
2 2
k4 = hf (tn + h, yn + k3 )
1
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )
6
Python Code:

10
1 def f (t , y ) :
2 return -2.0 * y
3

4 t0 , y0 = 0.0 , 1.0
5 h = 0.2
6

7 k1 = h * f ( t0 , y0 )
8 k2 = h * f ( t0 + h /2.0 , y0 + k1 /2.0)
9 k3 = h * f ( t0 + h /2.0 , y0 + k2 /2.0)
10 k4 = h * f ( t0 + h , y0 + k3 )
11

12 y1 = y0 + (1.0/6.0) * ( k1 + 2* k2 + 2* k3 + k4 )
13 print ( f " y (0.2) approx = { y1 :.5 f } " )

Problem 15: Systems of ODEs for a Damped Oscillator

d2 x
Question: The equation of motion of a damped oscillator is dt2
+ 0.5 dx
dt
+ 4x = 0. Convert
this equation into two first-order differential equations and outline a numerical scheme to
solve it using Runge-Kutta methods.
Explanation: To convert a second-order ODE into a system of first-order ODEs, we
introduce a new variable for the first derivative. Let the velocity be v = dxdt
. Then the
dv d2 x
acceleration is dt = dt2 . Substituting this into the original equation yields the coupled
system:
dx
=v
dt
dv
= −0.5v − 4x
dt
   
x v
We define a state vector Y⃗ = and its derivative function F⃗ (t, Y⃗ ) = . The
v −0.5v − 4x
standard RK4 scheme can then be applied to vectors.
Python Code:

1 import numpy as np
2

3 # System of ODEs returns an array [ dx / dt , dv / dt ]

4 def F (t , Y ) :
5 x , v = Y [0] , Y [1]
6 dxdt = v
7 dvdt = -0.5* v - 4.0* x
8 return np . array ([ dxdt , dvdt ])
9

11
10 def rk4_step (F , t , Y , h ) :
11 k1 = h * F (t , Y )
12 k2 = h * F ( t + h /2 , Y + k1 /2)
13 k3 = h * F ( t + h /2 , Y + k2 /2)
14 k4 = h * F ( t + h , Y + k3 )
15

16 Y_next = Y + (1.0/6.0) * ( k1 + 2* k2 + 2* k3 + k4 )
17 return Y_next
18

19 # Example initialization
20 t0 = 0.0
21 Y0 = np . array ([1.0 , 0.0]) # Initial position x =1 , velocity v =0
22 h = 0.1
23

24 Y1 = rk4_step (F , t0 , Y0 , h )
25 print ( f " After one step : x ={ Y1 [0]:.4 f } , v ={ Y1 [1]:.4 f } " )

Problem 16: Poisson Equation (Initial Value Problem)

d2 V
Question: Solve the one-dimensional Poisson equation dx2
= −ρ(x) on a discrete lattice
using finite difference methods for ρ(x) = 1 in the range [0, 2] with boundary conditions
V (0) = 0 and V ′ (0) = 1.0.
Explanation: Using the central finite difference approximation for the second deriva-
tive, we get:
Vi+1 − 2Vi + Vi−1
2
= −1 =⇒ Vi+1 = 2Vi − Vi−1 − h2
h
Since this is set up as an initial value problem, we are given V0 = 0 and V ′ (0) = 1.0. We
can use the forward difference approximation for the first derivative to find V1 :
V1 − V0
= 1.0 =⇒ V1 = V0 + h
h
From V0 and V1 , we can step forward to calculate the rest of the lattice points.
Python Code:

1 import numpy as np
2

3 h = 0.2
4 x = np . arange (0 , 2.0 + h , h )
5 V = np . zeros ( len ( x ) )
6

7 # Initial Conditions
8 V [0] = 0.0

12
9 V [1] = V [0] + h * 1.0 # From V ’(0) = 1.0
10

11 # Forward Stepping
12 for i in range (1 , len ( x ) - 1) :
13 V [ i +1] = 2* V [ i ] - V [i -1] - h **2
14

15 for i in range ( len ( x ) ) :

16 print ( f " x = { x [ i ]:.1 f } , V = { V [ i ]:.4 f } " )

Problem 17: Radioactive Decay Simulation

Question: A radioactive sample obeys N (t) = N0 e−λt . Using λ = 0.3 day−1 and N0 =
1000, generate numerical values of N (t) from t = 0 to 10 days with interval ∆t = 1 day.
Explanation: This is a straightforward evaluation of an exact analytical function over
a set of discrete time steps. We simply loop through t = 0, 1, 2, . . . , 10 and calculate N (t) at
each step.
Python Code:

1 import math
2

3 N0 = 1000.0
4 lam = 0.3
5

6 print ( " Time ( days ) \ tN ( t ) " )

7 for t in range (11) :
8 N = N0 * math . exp ( - lam * t )
9 print ( f " { t }\ t \ t { N :.1 f } " )

Problem 18: Discretizing the Schrödinger Equation

2
Question: A particle in a one-dimensional box satisfies − 12 ddxψ2 = Eψ. Discretize the
equation using finite differences and construct the corresponding Hamiltonian matrix.
Explanation: We approximate the second derivative using central finite differences:

d2 ψ ψi+1 − 2ψi + ψi−1

2
≈
dx h2
Substituting this into the Schrödinger equation gives:
1 1 1
− 2
ψi−1 + 2 ψi − 2 ψi+1 = Eψi
2h h 2h
13
 This forms an eigenvalue problem H ψ⃗ = E ψ,
⃗ where H is a tridiagonal matrix. The diagonal
elements are 1/h , and the adjacent off-diagonal elements are −1/(2h2 ).
2

Python Code:

1 import numpy as np
2

3 N = 5 # Number of interior lattice points

4 h = 1.0 / ( N + 1)
5 H = np . zeros (( N , N ) )
6

7 diagonal_val = 1.0 / ( h **2)

8 off_diagonal_val = -1.0 / (2 * h **2)
9

10 for i in range ( N ) :
11 H [i , i ] = diagonal_val
12 if i > 0:
13 H [i , i -1] = off_diagonal_val
14 if i < N - 1:
15 H [i , i +1] = off_diagonal_val
16

17 print ( " Hamiltonian Matrix ( H ) : " )

18 print ( np . round (H , 2) )

Problem 19: Numerical Differentiation for Electric Field

Question: Use numerical differentiation to estimate the electric field E(x) = − dV
dx
from
the potential data: x = {0, 1, 2, 3, 4}, V = {10, 8, 5, 2, 1}.
Explanation: The electric field is the negative gradient of the potential. We use
numerical derivatives with h = 1:

• Forward difference for the first point (x = 0): E0 ≈ − V1 −V

h
0

−Vi−1
• Central difference for interior points (x = 1, 2, 3): Ei ≈ − Vi+12h

• Backward difference for the last point (x = 4): E4 ≈ − V4 −V

h
3

Python Code:

1 x = [0 , 1 , 2 , 3 , 4]
2 V = [10.0 , 8.0 , 5.0 , 2.0 , 1.0]
3 h = 1.0
4 E = [0.0] * len ( x )
5

6 # Forward difference for endpoints

14
7 E [0] = -( V [1] - V [0]) / h
8 E [ -1] = -( V [ -1] - V [ -2]) / h
9

10 # Central difference for interior points

11 for i in range (1 , len ( x ) - 1) :
12 E [ i ] = -( V [ i +1] - V [i -1]) / (2 * h )
13

14 for i in range ( len ( x ) ) :

15 print ( f " x = { x [ i ]} , E = { E [ i ]} V / m " )

Problem 20: Projectile Trajectory (Euler Method)

Question: Simulate the trajectory of a projectile using Euler method under gravity. As-
sume: g = 9.8 m/s2 , v0 = 20 m/s, θ = 45◦ .
Explanation: The equations of motion are broken down into first-order ODEs:

dx dvx
= vx , =0
dt dt
dy dvy
= vy , = −g
dt dt
Using Euler’s method with a time step ∆t, the updates are: xnew = x+vx ∆t, ynew = y+vy ∆t,
and vy,new = vy − g∆t. We stop the simulation when the projectile hits the ground (y < 0).
Python Code:

1 import math
2

3 g = 9.8
4 v0 = 20.0
5 theta = math . radians (45)
6 dt = 0.05
7

8 x , y = 0.0 , 0.0
9 vx = v0 * math . cos ( theta )
10 vy = v0 * math . sin ( theta )
11 t = 0.0
12

13 print ( " Time ( s ) \ t x ( m ) \ t y ( m ) " )

14 while y >= 0:
15 print ( f " { t :.2 f }\ t { x :.2 f }\ t { y :.2 f } " )
16 x = x + vx * dt
17 y = y + vy * dt
18 vy = vy - g * dt
19 t = t + dt

15
 Problem 21: Monte Carlo Estimation of π
Question: Generate 104 random numbers uniformly distributed between 0 and 1. Use
them to estimate the value of π using the Monte Carlo method.
Explanation: We can estimate π by generating random points (x, y) in a 1 × 1 square
where 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1. We then check if the point falls inside the quarter-circle
of radius r = 1 by evaluating the condition x2 + y 2 ≤ 1. The ratio of points inside the
quarter-circle (M ) to the total number of points (N ) will approximate the ratio of their
areas:
Area of quarter circle π(1)2 /4 π
= =
Area of square 12 4
Thus, π ≈ 4 M
N
.
Python Code:

1 import random
2

3 N = 10000
4 M = 0
5

6 for _ in range ( N ) :
7 x = random . uniform (0 , 1)
8 y = random . uniform (0 , 1)
9

10 if x **2 + y **2 <= 1.0:

11 M += 1
12

13 pi_estimate = 4.0 * M / N
14 print ( f " Estimated value of pi : { pi_estimate :.4 f } " )

Problem 22: Canonical Ensemble Calculations

e−Ei /kB T
Question: A canonical ensemble has probabilities Pi = Z
. For energy levels E =
{0, 1, 2, 3} eV and kB T = 0.5 eV, calculate numerically: (a) The partition function Z. (b)
The average energy ⟨E⟩.
Explanation: The partition function Z is the sum of the Boltzmann factors over all
accessible states: X
Z= e−Ei /kB T
i

16
 The average energy (expectation value of energy) is given by:
X 1X
⟨E⟩ = Ei Pi = Ei e−Ei /kB T
i
Z i

We evaluate these sums discretely using the given arrays.

Python Code:

1 import math
2

3 E = [0.0 , 1.0 , 2.0 , 3.0]

4 kT = 0.5
5

6 Z = 0.0
7 for e in E :
8 Z += math . exp ( - e / kT )
9

10 avg_E = 0.0
11 for e in E :
12 P_i = math . exp ( - e / kT ) / Z
13 avg_E += e * P_i
14

15 print ( f " Partition Function Z = { Z :.5 f } " )

16 print ( f " Average Energy <E > = { avg_E :.5 f } eV " )

Problem 23: Newton-Raphson Method

Question: Use the Newton-Raphson method to determine the root of xe−x − 0.2 = 0
correct up to 10−6 .
Explanation: Let f (x) = xe−x − 0.2. We compute the derivative using the product
rule:
f ′ (x) = e−x − xe−x = e−x (1 − x)
The Newton-Raphson update step is:

f (xn ) xn e−xn − 0.2

xn+1 = xn − = xn −
f ′ (xn ) e−xn (1 − xn )

We must choose an initial guess x0 . Since f (0) = −0.2 and f (1) = e−1 − 0.2 ≈ 0.16, a root
exists between 0 and 1. We will use x0 = 0.5.
Python Code:

1 import math
2

17
3 def f ( x ) :
4 return x * math . exp ( - x ) - 0.2
5

6 def df ( x ) :
7 return math . exp ( - x ) * (1.0 - x )
8

9 x0 = 0.5
10 tol = 1e -6
11

12 while True :
13 x1 = x0 - f ( x0 ) / df ( x0 )
14 if abs ( x1 - x0 ) < tol :
15 break
16 x0 = x1
17

18 print ( f " Root : { x1 :.6 f } " )

Problem 24: Bisection Method for x = cos x

Question: Determine the solution of the nonlinear equation x = cos x using the bisection
iteration method. Discuss the convergence of the iteration.
Explanation: Let f (x) = x − cos x. Since f (0) = −1 and f (1) = 1 − cos(1) ≈ 0.459,
the root is bracketed in [0, 1]. The bisection method iteratively halves this interval. The
convergence of the bisection method is guaranteed and linear. The error bound decreases by
exactly a factor of 2 at each step, so the maximum absolute error after n iterations is b−a
2n
.
Python Code:

1 import math
2

3 def f ( x ) :
4 return x - math . cos ( x )
5

6 a , b = 0.0 , 1.0
7 tol = 1e -5
8 iterations = 0
9

10 while ( b - a ) / 2.0 > tol :

11 mid = ( a + b ) / 2.0
12 if f ( mid ) == 0:
13 break
14 elif f ( a ) * f ( mid ) < 0:
15 b = mid
16 else :

18
17 a = mid
18 iterations += 1
19

20 root = ( a + b ) / 2.0
21 print ( f " Root : { root :.5 f } found in { iterations } iterations . " )

Problem 25: Newton’s Forward Interpolation

Question: Use Newton’s forward interpolation formula to estimate y(2.5) from the data.
Estimate at which value of x, you will get y = 0.88.
Explanation: Newton’s Forward Interpolation formula is:

u(u − 1) 2 u(u − 1)(u − 2) 3

y = y0 + u∆y0 + ∆ y0 + ∆ y0 + . . .
2! 3!
where u = x−xh
0
. Here, h = 1 and x0 = 0. To find y(2.5), we construct the forward difference
table and apply the formula. To find x when y = 0.88, we can use a root-finding method
(like the secant method) on our constructed polynomial function P (x) − 0.88 = 0. Since
y(1) = 0.84 and y(2) = 0.91, the value y = 0.88 should occur for an x between 1 and 2.
Python Code:

1 import numpy as np
2 from scipy . optimize import fsolve
3

4 x = [0.0 , 1.0 , 2.0 , 3.0 , 4.0]

5 y = [0.0 , 0.84 , 0.91 , 0.14 , -0.76]
6 n = len ( x )
7

8 # Constructing the forward difference table

9 diff_table = np . zeros (( n , n ) )
10 diff_table [: , 0] = y
11

12 for j in range (1 , n ) :
13 for i in range ( n - j ) :
14 diff_table [ i ][ j ] = diff_table [ i +1][ j -1] - diff_table [ i ][ j -1]
15

16 def newto n_ f o rw a r d_ e v al ( target_x ) :

17 h = x [1] - x [0]
18 u = ( target_x - x [0]) / h
19 result = diff_table [0][0]
20 u_term = 1.0
21 factorial = 1.0
22

23 for i in range (1 , n ) :

19
24 u_term *= ( u - ( i - 1) )
25 factorial *= i
26 result += ( u_term * diff_table [0][ i ]) / factorial
27 return result
28

29 y_25 = new t o n_ f o rw a r d_ e v al (2.5)

30 print ( f " Estimated y (2.5) : { y_25 :.4 f } " )
31

32 # Finding x for y = 0.88

33 # We know from data that y =0.88 is roughly between x =1 and x =2
34 def root_func ( guess_x ) :
35 return n ew t o n_ f o rw a r d_ e v al ( guess_x ) - 0.88
36

37 x_target = fsolve ( root_func , 1.5) [0]

38 print ( f " Estimated x for y =0.88: { x_target :.4 f } " )

Problem 26: Least Squares Linear Fitting

Question: Fit a straight line y = ax + b to the following data using the least squares
method: x = {1, 2, 3, 4, 5}, y = {2.1, 4.2, 5.9, 8.1, 9.8}.
Explanation: To fit a line y = ax + b, we minimize the sum of squared errors. This
leads to the normal equations:
X X
b·n+a xi = yi
X X X
b xi + a x2i = xi y i

By calculating the sums from the given data (n = 5), we can solve this 2 × 2 system of linear
equations to find the slope a and intercept b. Analytically, the solutions are:
P P P P P
n xy − x y y−a x
a= P P , b=
n x2 − ( x)2 n
Python Code:

1 x = [1.0 , 2.0 , 3.0 , 4.0 , 5.0]

2 y = [2.1 , 4.2 , 5.9 , 8.1 , 9.8]
3 n = len ( x )
4

5 sum_x = sum ( x )
6 sum_y = sum ( y )
7 sum_xx = sum ( xi **2 for xi in x )
8 sum_xy = sum ( xi * yi for xi , yi in zip (x , y ) )
9

20
10 # Using the analytical formulas
11 denominator = n * sum_xx - sum_x **2
12 a = ( n * sum_xy - sum_x * sum_y ) / denominator
13 b = ( sum_y - a * sum_x ) / n
14

15 print ( f " Fitted Line : y = { a :.2 f } x + { b :.2 f } " )

Problem 27: Exponential Curve Fitting

Question: The decay of current in an RL circuit is measured. Fit the data numerically to
the model I(t) = I0 e−Rt/L .
Explanation: We can linearize the exponential model by taking the natural logarithm
of both sides:  
R
ln(I) = ln(I0 ) − t
L
Let Y = ln(I), A = −R/L, and B = ln(I0 ). The equation becomes a linear relation
Y = At + B. We calculate Yi = ln(Ii ) for the given data and perform a standard linear
least-squares fit to find A and B. Finally, I0 = eB and the decay constant is −A.
Python Code:

1 import math
2

3 t = [0.0 , 1.0 , 2.0 , 3.0 , 4.0]

4 I = [10.0 , 6.1 , 3.7 , 2.2 , 1.4]
5 n = len ( t )
6

7 # Linearize the data

8 Y = [ math . log ( i ) for i in I ]
9

10 sum_t = sum ( t )
11 sum_Y = sum ( Y )
12 sum_tt = sum ( ti **2 for ti in t )
13 sum_tY = sum ( ti * Yi for ti , Yi in zip (t , Y ) )
14

15 # Least squares for Y = At + B

16 denominator = n * sum_tt - sum_t **2
17 A = ( n * sum_tY - sum_t * sum_Y ) / denominator
18 B = ( sum_Y - A * sum_t ) / n
19

20 I0 = math . exp ( B )
21 decay_constant = -A # This represents R / L
22

21
23 print ( f " Fitted Model : I ( t ) = { I0 :.2 f } * exp ( -{ decay_constant :.4 f } *
t)")

Problem 28: Simpson’s 1/3 Rule Variation

R3 dx
Question: Use Simpson 1/3 rule to evaluate 0 1+x
. Check the variation of results with
no of intervals.
Explanation: Simpson’s 1/3 rule requires an even number of intervals n. We will write
a loop to evaluate the integral for multiple values of n (e.g., n = 2, 4, 6, 10) to observe how
the numerical result converges toward the exact analytical value, which is ln(1 + 3) − ln(1) =
ln(4) ≈ 1.38629.
Python Code:

1 import math
2

3 def f ( x ) :
4 return 1.0 / (1.0 + x )
5

6 a , b = 0.0 , 3.0
7 exact = math . log (4.0)
8

9 print ( " n \ t Integral \ t Error " )

10 for n in [2 , 4 , 6 , 8 , 10]:
11 h = (b - a) / n
12 integral = f ( a ) + f ( b )
13

14 for i in range (1 , n ) :
15 x_i = a + i * h
16 if i % 2 == 0:
17 integral += 2 * f ( x_i )
18 else :
19 integral += 4 * f ( x_i )
20

21 integral *= h / 3.0
22 error = abs ( exact - integral )
23 print ( f " { n }\ t { integral :.6 f }\ t { error :.2 e } " )

Problem 29: Trapezoidal Error Bound

R2
Question: Evaluate 0
x2 e−x dx using the trapezoidal rule with h = 0.25. Is this h sufficient

22
 for 3rd decimal place accuracy? If not, find the required h.
Explanation: First, we evaluate the integral with h = 0.25. To check accuracy, we
2
analyze the error bound |E| ≤ (b−a)h 12
max |f ′′ (x)|. For f (x) = x2 e−x , the second derivative
is f ′′ (x) = (x2 −4x+2)e−x . On the interval [0, 2], max |f ′′ (x)| = 2 (at x = 0). The maximum
2
error for h = 0.25 is E ≤ (2)(0.25)
12
(2) ≈ 0.0208, which is much larger than 0.0005 (required for
2
3rd decimal place accuracy). To find the required h: 2·h 12
· 2 ≤ 0.0005 =⇒ h2 ≤ 0.0015 =⇒
h ≤ 0.0387. A step size of h = 0.02 or smaller is required.
Python Code:

1 import math
2

3 def f ( x ) :
4 return ( x **2) * math . exp ( - x )
5

6 a , b = 0.0 , 2.0
7 h = 0.25
8 n = int (( b - a ) / h )
9

10 integral = f ( a ) + f ( b )
11 for i in range (1 , n ) :
12 integral += 2 * f ( a + i * h )
13 integral *= h / 2.0
14

15 exact = 2.0 - 10.0 * math . exp ( -2.0)

16 print ( f " Integral ( h =0.25) : { integral :.4 f } " )
17 print ( f " Exact Value : { exact :.4 f } " )
18 print ( " h =0.25 is NOT sufficient for 3 rd decimal place accuracy . " )
19

20 # Required h calculation
21 M = 2.0 # Max of | f ’ ’( x ) | on [0 , 2]
22 desired_error = 0.0005
23 required_h = math . sqrt ((12 * desired_error ) / (( b - a ) * M ) )
24 print ( f " Required step size h <= { required_h :.4 f } " )

Problem 30: Modified Euler Method

dy
Question: Solve dx
= y − x2 with y(0) = 1 using the modified Euler method up to x = 0.4
with h = 0.1.
Explanation: The modified Euler method (also known as Heun’s method) improves
upon standard Euler by taking an average of the slopes at the beginning and the estimated
end of the interval.

23
 (0)
1. Predictor: Estimate the next y using standard Euler: yn+1 = yn + hf (xn , yn )
(0)
2. Corrector: Refine the estimate: yn+1 = yn + h2 [f (xn , yn ) + f (xn+1 , yn+1 )]

Python Code:

1 def f (x , y ) :
2 return y - x **2
3

4 x = 0.0
5 y = 1.0
6 h = 0.1
7 x_target = 0.4
8

9 print ( f " x = { x :.1 f } , y = { y :.5 f } " )

10

11 while round (x , 1) < x_target :

12 # Predictor step
13 slope1 = f (x , y )
14 y_predict = y + h * slope1
15

16 # Corrector step
17 x_next = x + h
18 slope2 = f ( x_next , y_predict )
19 y_next = y + ( h / 2.0) * ( slope1 + slope2 )
20

21 # Update variables
22 x = x_next
23 y = y_next
24 print ( f " x = { x :.1 f } , y = { y :.5 f } " )

Problem 31: Coupled ODEs (RK4 Method)

dx dy
Question: Solve the coupled equations dt
= y, dt
= −x using the Runge-Kutta fourth-
order method for one time step in the range [0, π]. Check the change in solution with the
step size and learn to choose suitable h.
Explanation: This system describes a simple harmonic oscillator. We can represent
⃗
it as a vector ODE: ddtY = F⃗ (t, Y⃗ ), where Y⃗ = [x, y]T and F⃗ (t, Y⃗ ) = [y, −x]T . Since initial
conditions are not explicitly provided in the problem text, we will assume standard initial
conditions x(0) = 1 and y(0) = 0 (which corresponds to x(t) = cos(t) and y(t) = − sin(t)).
We will evaluate one step for different values of h to observe the truncation error.
Python Code:

24
1 import numpy as np
2 import math
3

4 def F (t , Y ) :
5 return np . array ([ Y [1] , -Y [0]])
6

7 def rk4_step (F , t , Y , h ) :
8 k1 = h * F (t , Y )
9 k2 = h * F ( t + h /2.0 , Y + k1 /2.0)
10 k3 = h * F ( t + h /2.0 , Y + k2 /2.0)
11 k4 = h * F ( t + h , Y + k3 )
12 return Y + (1.0/6.0) * ( k1 + 2* k2 + 2* k3 + k4 )
13

14 Y0 = np . array ([1.0 , 0.0])

15 t0 = 0.0
16

17 print ( " h \ t \ t x_approx \ t exact_x \ t Error " )

18 for h in [0.1 , 0.5 , 1.0 , math . pi ]:
19 Y1 = rk4_step (F , t0 , Y0 , h )
20 exact_x = math . cos ( h )
21 error = abs ( exact_x - Y1 [0])
22 print ( f " { h :.4 f }\ t \ t { Y1 [0]:.5 f }\ t { exact_x :.5 f }\ t { error :.2 e } " )

Problem 32: Logistic Equation

dN N
Question: The logistic equation is given by dt
= rN (1 − K
). Using r = 0.5, K = 100,
and N (0) = 10, solve numerically up to t = 10.
Explanation: The logistic equation models population growth with a carrying capacity
K. We can solve this first-order ODE numerically using Euler’s method or the Runge-Kutta
method. For simplicity and sufficient accuracy over this range, we use
 the standard

 Euler
Ni
method with a small step size (h = 0.1). Update rule: Ni+1 = Ni + h rNi 1 − K .
Python Code:

1 r = 0.5
2 K = 100.0
3 N = 10.0
4 t = 0.0
5 t_end = 10.0
6 h = 0.1
7

8 print ( " Time \ t Population ( N ) " )

9 while t <= t_end + 1e -5:

25
10 if round (t , 1) % 1.0 == 0: # Print every 1 time unit
11 print ( f " { t :.1 f }\ t { N :.2 f } " )
12

13 # Euler method step

14 dNdt = r * N * (1.0 - N / K )
15 N = N + h * dNdt
16 t = t + h

Problem 33: Time-Independent Heat Equation

d2 T
Question: Solve the time-independent heat equation dx2
= 0 for a rod of length L = 1 m
with boundary conditions T (0) = 100◦ C, T ′ (0) = 2.0. Plot the temperature distribution.
2
Explanation: The equation ddxT2 = 0 implies that the temperature gradient is constant.
Integrating twice gives the analytical solution T (x) = C1 x + C2 . Applying the boundary
conditions: 1. T (0) = 100 =⇒ C2 = 100 2. T ′ (0) = 2.0 =⇒ C1 = 2.0 Thus, the exact
temperature distribution is purely linear: T (x) = 2.0x + 100.
Python Code:

1 import numpy as np
2 # Note : In a real environment , uncomment the matplotlib lines to see
the plot
3 # import matplotlib . pyplot as plt
4

5 L = 1.0
6 x = np . linspace (0 , L , 100)
7

8 # Exact solution based on initial conditions

9 T = 2.0 * x + 100.0
10

11 print ( " x ( m ) \ t T ( C ) " )

12 for i in [0 , 25 , 50 , 75 , 99]: # Print a few sample points
13 print ( f " { x [ i ]:.2 f }\ t { T [ i ]:.2 f } " )
14

15 # plt . plot (x , T , label = ’ T ( x ) = 2 x + 100 ’)

16 # plt . xlabel ( ’ Position x ( m ) ’)
17 # plt . ylabel ( ’ Temperature T ( C ) ’)
18 # plt . title ( ’ Temperature Distribution in Rod ’)
19 # plt . grid ( True )
20 # plt . show ()

26
 Problem 34: 1D Random Walk Simulation
Question: A particle performs a one-dimensional random walk. Write an algorithm to
simulate 1000 steps and calculate: (a) Mean displacement. (b) Mean square displacement.
Explanation: To find the statistical mean displacement ⟨x⟩ and mean square dis-
placement ⟨x2 ⟩, we must simulate an ensemble of many independent random walkers (e.g.,
M = 1000 particles), each taking N = 1000 steps. At each step, a particle moves either
+1 or −1 with equal probability. For an unbiased random walk, theory predicts: 1. Mean
displacement ⟨x⟩ ≈ 0 2. Mean square displacement ⟨x2 ⟩ ≈ N = 1000
Python Code:

1 import random
2

3 n_steps = 1000
4 n_particles = 1000 # Number of trials to average over
5

6 sum_x = 0.0
7 sum_x_squared = 0.0
8

9 for _ in range ( n_particles ) :

10 x = 0
11 for _ in range ( n_steps ) :
12 # Move +1 or -1 with 50% probability
13 step = 1 if random . random () > 0.5 else -1
14 x += step
15

16 sum_x += x
17 sum_x_squared += x **2
18

19 mean_x = sum_x / n_particles

20 mean_x_sq = sum_x_squared / n_particles
21

22 print ( f " Number of steps : { n_steps } " )

23 print ( f " Mean Displacement <x >: { mean_x :.2 f } " )
24 print ( f " Mean Square Displacement <x ^2 >: { mean_x_sq :.2 f } " )

Problem 35: Monte Carlo Integration

R1√
Question: Use the Monte Carlo integration method to evaluate 0
1 − x2 dx and hence
estimate the value of π.
Explanation: The integral represents the area of a quarter circle of radius 1, so
π
analytically I = 4
, meaning π = 4I. Using the Monte Carlo Mean Value Theorem for

27
 integration, the integral of f (x) over [a, b] can be approximated by generating N uniformly
distributed random numbers xi in [a, b] and computing the average function value:
N
1 X
I ≈ (b − a) f (xi )
N i=1
√
Here, f (x) = 1 − x2 , a = 0, and b = 1.
Python Code:

1 import random
2 import math
3

4 N = 100000
5 sum_f = 0.0
6

7 for _ in range ( N ) :
8 x = random . uniform (0 , 1)
9 f_x = math . sqrt (1.0 - x **2)
10 sum_f += f_x
11

12 integral = sum_f / N
13 pi_estimate = 4.0 * integral
14

15 print ( f " Evaluated Integral : { integral :.5 f } " )

16 print ( f " Estimated value of pi : { pi_estimate :.5 f } " )

Problem 36: Two-Level System Thermodynamics

Question: The partition function of a two-level system is Z = 1 + e−ϵ/kB T . For ϵ = 1
eV, calculate numerically: (a) Partition function, (b) Probability of occupation, (c) Average
energy at T = 300 K and T = 1000 K.
Explanation: For a two-level system with energies E0 = 0 and E1 = ϵ:

1. Partition function: Z = 1 + e−ϵ/kB T

1 e−ϵ/kB T
2. Probabilities: P0 = Z
and P1 = Z

3. Average Energy: ⟨E⟩ = E0 P0 + E1 P1 = 0 · P0 + ϵ · P1 = ϵP1

We use the Boltzmann constant kB ≈ 8.617 × 10−5 eV/K to ensure consistent units (eV).
Python Code:

28
1 import math
2

3 kB = 8.617 e -5 # eV / K
4 epsilon = 1.0 # eV
5 temperatures = [300 , 1000]
6

7 for T in temperatures :
8 print ( f " --- Temperature : { T } K ---" )
9 kT = kB * T
10

11 # ( a ) Partition function
12 Z = 1.0 + math . exp ( - epsilon / kT )
13 print ( f " Partition Function Z : { Z :.5 e } " )
14

15 # ( b ) Probabilities
16 P0 = 1.0 / Z
17 P1 = math . exp ( - epsilon / kT ) / Z
18 print ( f " P ( State 0) : { P0 :.5 f } , P ( State 1) : { P1 :.5 e } " )
19

20 # ( c ) Average Energy
21 avg_E = epsilon * P1
22 print ( f " Average Energy : { avg_E :.5 e } eV \ n " )

Problem 37: Radial Schrödinger Equation for Hydrogen

2
Question: The radial equation is − 12 ddru2 + [− 1r ]u = 0.5u. Solve this using Euler method in
the range [0, 5]. Given u(0) = 0 and u′ (0) = 1.0.
Explanation: We rearrange the equation to isolate the second derivative:

d2 u
   
1 2
= −2 0.5 + u=− 1+ u
dr2 r r

We convert this to a system of first-order ODEs by setting v = du . Then dv = − 1 + 2r u.

dr dr
At r = 0, the term ur is indeterminate (0/0). Using L’Hopital’s rule, limr→0 u(r)
r
= u′ (0) = 1.0.
To avoid a division-by-zero error in the code, we can start our numerical integration at a very
small non-zero radius r = h, using the initial conditions to approximate u(h) ≈ h · u′ (0) = h
and v(h) ≈ 1.0.
Python Code:

1 h = 0.01
2 r_end = 5.0
3

29
4 # Start at r = h to avoid division by zero
5 r = h
6 u = h * 1.0 # u ( h ) approx h * u ’(0)
7 v = 1.0 # v = du / dr
8

9 print ( " r \ t \ t u ( r ) " )

10 while r <= r_end :
11 if round (r , 2) % 1.0 == 0: # Print every 1.0 interval
12 print ( f " { r :.2 f }\ t { u :.4 f } " )
13

14 # Euler method step

15 du_dr = v
16 dv_dr = -(1.0 + 2.0 / r ) * u
17

18 u = u + h * du_dr
19 v = v + h * dv_dr
20 r = r + h

Problem 38: Minimum via Bisection Method

Question: Use the bisection method to determine the minimum of f (x) = x2 − 4x + 5 in
the interval 0 ≤ x ≤ 5 by solving df /dx = 0.
Explanation: To find the minimum, we find the root of the first derivative.

df
g(x) = = 2x − 4 = 0
dx
We check the bounds: g(0) = −4 and g(5) = 6. Since the signs are opposite, a root is
bracketed in [0, 5]. We apply the standard bisection algorithm to g(x).
Python Code:

1 def g ( x ) :
2 return 2 * x - 4
3

4 a , b = 0.0 , 5.0
5 tol = 1e -5
6

7 while ( b - a ) / 2.0 > tol :

8 mid = ( a + b ) / 2.0
9 if g ( mid ) == 0:
10 break
11 elif g ( a ) * g ( mid ) < 0:
12 b = mid
13 else :

30
14 a = mid
15

16 min_x = ( a + b ) / 2.0
17 print ( f " Location of minimum : x = { min_x :.5 f } " )

Problem 39: Minimum via Newton-Raphson Method

Question: Use the Newton-Raphson method to locate the minimum of f (x) = x4 −3x2 +2.
Start with the initial guess x0 = 1.
Explanation: To find the extremum, we must find the root of the first derivative f ′ (x).
The Newton-Raphson formula for finding the root of f ′ (x) = 0 requires the second derivative
f ′′ (x):
f ′ (xn )
xn+1 = xn − ′′
f (xn )
For our function:

f ′ (x) = 4x3 − 6x
f ′′ (x) = 12x2 − 6

We iterate this starting from x0 = 1. Note: f ′′ (x) > 0 at the root confirms it is a local
minimum.
Python Code:

1 def df ( x ) :
2 return 4 * x **3 - 6 * x
3

4 def ddf ( x ) :
5 return 12 * x **2 - 6
6

7 x0 = 1.0
8 tol = 1e -6
9 iterations = 0
10

11 while True :
12 x1 = x0 - df ( x0 ) / ddf ( x0 )
13 iterations += 1
14 if abs ( x1 - x0 ) < tol :
15 break
16 x0 = x1
17

18 print ( f " Minimum located at x = { x1 :.5 f } after { iterations }

iterations . " )

31
 Problem 40: Gradient Descent Optimization
Question: Use the gradient descent method to determine the minimum of f (x) = x2 +
3x + 2. Use learning rate η = 0.1 and initial value x0 = 5.
Explanation: Gradient descent minimizes a function by iteratively moving in the
direction of steepest descent (the negative gradient). The update rule is:

df
xn+1 = xn − η (xn )
dx
The derivative is f ′ (x) = 2x + 3. With η = 0.1 and x0 = 5, we iterate until the step size
becomes smaller than a defined tolerance.
Python Code:

1 def df ( x ) :
2 return 2 * x + 3
3

4 x = 5.0
5 eta = 0.1
6 tol = 1e -6
7 iterations = 0
8

9 while True :
10 gradient = df ( x )
11 step = eta * gradient
12 x_new = x - step
13 iterations += 1
14

15 if abs ( step ) < tol :

16 break
17 x = x_new
18

19 print ( f " Minimum found at x = { x :.5 f } after { iterations } iterations . "

)

Problem 41: Equilibrium Points and Stability

Question: The potential energy of a particle is given by U (x) = x4 − 2x2 + 1. Find all
equilibrium points numerically and determine which correspond to stable equilibrium.
Explanation: Equilibrium points occur where the force is zero, meaning the first

32
 derivative of the potential energy is zero: U ′ (x) = 4x3 − 4x = 0. To determine stability, we
evaluate the second derivative U ′′ (x) = 12x2 − 4 at these points. If U ′′ (x) > 0, the point
is a local minimum (stable equilibrium). If U ′′ (x) < 0, it is a local maximum (unstable
equilibrium). We can find the roots numerically by scanning the domain to bracket roots,
then using the Bisection or Newton-Raphson method.
Python Code:

1 def dU ( x ) :
2 return 4 * x **3 - 4 * x
3

4 def ddU ( x ) :
5 return 12 * x **2 - 4
6

7 # Known analytically to be around -1 , 0 , 1. We provide guesses .

8 guesses = [ -1.5 , -0.1 , 1.5]
9 tol = 1e -6
10

11 print ( " Equilibrium Points : " )

12 for guess in guesses :
13 x = guess
14 while abs ( dU ( x ) ) > tol :
15 # Newton - Raphson for finding root of dU ( x ) = 0
16 x = x - dU ( x ) / ddU ( x )
17

18 stability = " Stable " if ddU ( x ) > 0 else " Unstable "
19 print ( f " x = { x :.4 f } \ t U ’ ’( x ) = { ddU ( x ) :.4 f } \ t ({ stability }) " )

Problem 42: Minimum via Newton-Raphson

Question: Use Newton-Raphson to determine the minimum of f (x) = ex + x2 .
Explanation: To find the minimum of f (x), we must find the root of its first derivative,
f ′ (x) = ex + 2x = 0. The Newton-Raphson method applied to f ′ (x) requires the second
derivative, f ′′ (x) = ex + 2. The update formula becomes:
f ′ (xn ) exn + 2xn
xn+1 = xn − = xn −
f ′′ (xn ) e xn + 2
Python Code:

1 import math
2

3 def df ( x ) :
4 return math . exp ( x ) + 2 * x
5

33
6 def ddf ( x ) :
7 return math . exp ( x ) + 2
8

9 x0 = -0.5 # Initial guess

10 tol = 1e -6
11 iterations = 0
12

13 while True :
14 x1 = x0 - df ( x0 ) / ddf ( x0 )
15 iterations += 1
16 if abs ( x1 - x0 ) < tol :
17 break
18 x0 = x1
19

20 print ( f " Minimum located at x = { x1 :.6 f } after { iterations }

iterations . " )

Problem 43: Minimum via Gradient Descent

Question: Use the gradient descent method to locate the minimum of f (x) = x3 − 6x2 +
9x + 1 by solving f ′ (x) = 0.
Explanation: Gradient descent approaches the local minimum by taking steps propor-
tional to the negative of the gradient. The gradient is f ′ (x) = 3x2 − 12x + 9. The update
rule is xn+1 = xn − ηf ′ (xn ). We choose a small learning rate η and an initial guess. Based
on the coefficients, roots of the derivative are at x = 1 and x = 3. To find the minimum, we
start with a guess like x0 = 4.
Python Code:

1 def df ( x ) :
2 return 3 * x **2 - 12 * x + 9
3

4 x = 4.0 # Initial guess

5 eta = 0.05 # Learning rate
6 tol = 1e -6
7 iterations = 0
8

9 while True :
10 gradient = df ( x )
11 step = eta * gradient
12 x_new = x - step
13 iterations += 1
14

15 if abs ( step ) < tol :

34
16 break
17 x = x_new
18

19 print ( f " Minimum found at x = { x :.5 f } after { iterations } iterations . "

)

Problem 44: Minimum of a Periodic Potential

2
Question: A particle moves in the potential V (x) = sin x + x10 . Use numerical methods to
determine the position of the minimum potential between x = −4 and x = 4.
Explanation: To find the minimum, we must solve for V ′ (x) = cos x + x5 = 0. We can
plot or scan the function to find that the global minimum in [−4, 4] occurs at a negative x
value (around x ≈ −1.3). We will use the Newton-Raphson method starting from x0 = −1.0
to precisely locate this minimum. The second derivative is V ′′ (x) = − sin x + 0.2.
Python Code:

1 import math
2

3 def dV ( x ) :
4 return math . cos ( x ) + x / 5.0
5

6 def ddV ( x ) :
7 return - math . sin ( x ) + 0.2
8

9 x0 = -1.0 # Initial guess near the minimum

10 tol = 1e -6
11 iterations = 0
12

13 while True :
14 x1 = x0 - dV ( x0 ) / ddV ( x0 )
15 iterations += 1
16 if abs ( x1 - x0 ) < tol :
17 break
18 x0 = x1
19

20 min_V = math . sin ( x1 ) + ( x1 **2) / 10.0

21 print ( f " Minimum located at x = { x1 :.5 f } " )
22 print ( f " Potential energy at minimum : V ( x ) = { min_V :.5 f } " )

35