Empowering Scientific Computing with

Python: Numerical Analysis and

Comparative Advantages

Sourabh Kumar Das1,∗

1 Dept of Physics, Raja N.L Khan Women’s College (Autonomous),

Gope Palace, Midnapore,721102,West Bengal, India.

skd@ rnlkwc. ac. in

1. Python for Numerical Analysis: A Comparative Perspective

Python has established itself as a leading programming language for numerical analy-
sis, offering distinct advantages over traditional languages such as C, C++, Fortran, or
proprietary tools like MATLAB. Its suitability can be highlighted through the following
points:

1. Simplicity and Readability: Python’s high-level, human-readable syntax re-

duces programming complexity, enabling researchers and students to implement
numerical methods with fewer lines of code compared to low-level languages.

2. Comprehensive Scientific Ecosystem: Python provides robust libraries such

as NumPy for array-based operations, SciPy for advanced numerical computations,
and SymPy for symbolic mathematics. These libraries are highly optimized and rival
the computational efficiency of traditional compiled languages.

3. Effective Data Visualization: Numerical analysis often requires graphical inter-

pretation. Python’s libraries, including Matplotlib, Seaborn, and Plotly, facili-
tate high-quality visualization, thereby supporting data exploration and presenta-
tion in research contexts.

4. Interdisciplinary Applicability: Beyond numerical computation, Python inte-

grates seamlessly with domains such as data science, artificial intelligence, and
engineering simulations. This interdisciplinary adaptability distinguishes it from
domain-specific languages.

1
 5. Open-Source and Community Support: Python benefits from continuous de-
velopment, extensive documentation, and a vibrant open-source community. This
ensures accessibility, reproducibility, and collaborative advancement, which may
not be as pronounced in proprietary or less widely adopted alternatives.

In summary, Python’s combination of readability, efficiency, specialized libraries, and

broad applicability establishes it as a superior choice for numerical analysis when com-
pared with other programming languages.

• open source programming language

• many built in libraries

• cross platform language

• large user network and free support

Q: What is the editor that we want to use to write Python programs and
execute them?

A: There are several editors and IDEs (Integrated Development Environments) available
for writing and executing Python programs. The choice depends on the purpose:

• If you are just learning Python basics: Use Thonny (very beginner-friendly) or
IDLE (comes with Python).

• If you are doing research, data science, or numerical analysis: Use Jupyter
Notebook / JupyterLab, as it allows step-by-step execution, equations, and visu-
alizations.

• If you want to develop larger applications or professional projects: Use VS

Code (lightweight, flexible) or PyCharm (a full professional IDE).

Recommended Choice: Since the focus is on numerical analysis and scientific comput-
ing, the best option is to use Jupyter Notebook for experimentation and documenta-
tion, together with VS Code for clean, reusable Python scripts. There are many editor
like Visual studio/ anaconda.
Q: why we use COLAB when we have many other editor of python ?

• integrated with Google drive

• connects to the cloud based run time.

• easy sharing

2
1.1. Uses of array:

One-Dimensional Array (List)

# One Dimensional Array

arr = [ 1 0 , 20 , 30 , 40 , 50]

print ( " One␣ Dimensional ␣ Array ␣ e l e m e n t s : " )

for i in a r r :
print ( i )

Two-Dimensional Array (Matrix)

# Two Dimensional Array

arr2D = [
[1 , 2 , 3] ,
[4 , 5 , 6] ,
[7 , 8 , 9]
]

print ( "Two␣ Dimensional ␣ Array ␣ e l e m e n t s : " )

for row in arr2D :
for element in row :
print ( element , end=" ␣ " )
print ( )

Output

One Dimensional Array elements:

10
20
30
40
50

Two Dimensional Array elements:

1 2 3
4 5 6
7 8 9

3
matrix multiplication

import numpy as np
from s c i p y import l i n a l g

A = np . a r r a y ( [ [ 1 , 2 , 3 ] , [ 4 , 5 , 6 ] , [ 7 , 8 , 9 ] ] )
B = np . a r r a y ( [ [ 9 , 8 , 7 ] , [ 6 , 5 , 4 ] , [ 3 , 2 , 1 ] ] )

C = l i n a l g . b l a s . dgemm ( 1 . 0 , A, B) # Matrix m u l t i p l i c a t i o n
print (C)

Explanation

1. import numpy as np Loads the NumPy library (renamed as np). It is used here
to create arrays (matrices A and B).

2. from scipy import linalg Imports SciPy’s linear algebra module, which provides
optimized routines for matrix operations.

3. Matrix definitions:    

1 2 3 
9 8 7
A = 4 5 6 , B = 6 5 4
   
   
7 8 9 3 2 1

4. [Link](1.0, A, B) Here, dgemm means “Double-precision General

Matrix-Matrix multiplication” (a BLAS routine). It computes:

C = 1.0 × (A · B)

5. print(C) Displays the resulting matrix.

Result

The multiplication yields:  


30 24 18
C =  84 69 54
 
 
138 114 90

4
2. Q: Why should we learn about Numerical methods and Sim-
ulation technique ?

• Most of the physical problems do not have analytical solution. In order to gethrer
understanding of the system , we need computational technique to arrive at the so-
lution. These solutions when analysed , rendered understanding of the key concept
that were embeded into the equations.

• They are powerful methods to handle large systems of equations with several un-
known , large matrices, complicated nonlinearites etc.

• Different computational techniques provide efficient use of the computer by learning

computer languages such as C, Python, Fortan etc.

• Many of the software packages which we use on the daily basis for our computing,
plotting, documentation needs employ numerical techniques.

• Numerical technique provide an idea of how errors propagates, especially in high

accuracy computation.

3. Roots of algebric and Transcendental equations:

Suppose for a function f (x) of variable x, we can find a root (x = ξ) of the equation when
f (ξ) = 0.
If the function f (x) is quadratic or cubic, biquadratic or some other form , we know how
to analytically solve the equation through the algebric formula.
If the function f (x) is a polynomial of higher degree or the equation involves some tran-
scendental or trigonometrix function like ex , lnx, sinx, tanx, sinhx etc transcendental
equation, then no formula exist and we find the roots through various approximate nu-
merical methods. In order to efficiently solve by the approximate methods, we need the
help of computer.

3.1. Bisection method:

If a function f (x) is a continuous between the two limits x = a and b, and f (a) and f (b)
are of opposite signs , then there exists a root between the a and b. Let us think that
the root is at x = ξ so that f (ξ) = 0.
The root is found between a and b if the following condition holds:

f (a) × f (b) < 0 (1)

5
 Figure 1: The root is bracketed in [a,b]

The root can be approximated by the midpoint , xm = (a+b) 2

. That means we bisect the
interval. Now in case , f (xm ) = 0 , then that is the root. If not, we search the root in
either of the intervals: [a, xm ] or [xm , b]. In this case ,

f (a) × f (xm ) < 0 (2)

The root lies in the first half , otherwise the root lies in the second half. Then as before
we bisect the new interval and repeat the process until the root is found according to the
desired accuracy.

6
Algorithm 1 Bisection Method for Root Finding
1: Input: Continuous function f (x), interval [a, b], tolerance ϵ
2: Output: Approximate root c
3: if f (a) · f (b) ≥ 0 then
4: Print “Invalid interval: Root may not exist”
5: Stop
6: end if
7: while b−a 2
> ϵ do
a+b
8: c← 2
9: if f (c) = 0 then
10: return c {Exact root found}
11: else if f (a) · f (c) < 0 then
12: b←c
13: else
14: a←c
15: end if
16: end while
17: return a+b 2
{Approximate root}

import numpy as np
from [Link] import bisect

# Define the function

def f(x):
return x**3 - (2*x) - 5 # f(x) = x^3 - 2x - 5

# Interval [a, b] where f(a)*f(b) < 0

a, b = 2, 3 # root lies between 2 and 3

# Use SciPy’s bisect method

root = bisect(f, a, b, xtol=1e-6)

print("Root found:", root)

import numpy as np

def f(x):
return x**3 - 2*x - 5 # f(x) = x^3 - 2x - 5

a, b = 2, 3 # root lies between 2 and 3

tol = 1e-6

7
while (b - a)/2 > tol:
c = (a + b) / 2
if f(c) == 0:
break
elif f(a) * f(c) < 0:
b = c
else:
a = c

root = (a + b) / 2
print("Root found:", root)

from sympy import symbols, Eq, nsolve

x = symbols(’x’)
f = x**3 - (2*x) - 5 # f(x) = x^3 - 2x - 5

# Use nsolve with an initial guess

root = nsolve(Eq(f, 0), 2)

print("Root found:", root)

We want to find the root of the equation

f (x) = x3 − 2x − 5

using SymPy. The program can be explained step by step as follows:

1. Import required functions

from sympy import symbols, Eq, nsolve

• symbols is used to define symbolic variables.

• Eq represents an equation, e.g. f (x) = 0.

• nsolve is a numerical solver that finds approximate solutions.

2. Define the variable

8
 x = symbols(’x’)

Here, x is declared as a symbolic variable (not just a number), required for building
symbolic expressions.

3. Define the function

f = x**3 - 2*x - 5

This defines f (x) = x3 − 2x − 5 as a symbolic expression.

4. Call the solver

root = nsolve(Eq(f, 0), 2)

• Eq(f, 0) tells SymPy to solve the equation f (x) = 0.

• The initial guess is 2, since the real root lies near 2.094.

• nsolve computes a numerical solution.

5. Print the result

print("Root found:", root)

This displays the computed root, approximately

Root found: 2.09455148154233.

3.2. Secant Method

Theory:
To search the root, we start from two arbitrary points: a and b. We draw a secant
between the two coordinates:[a, f (a)] and [b, f (b)].
The equation for any line passing through the above two points in the plane:

y − f (a) f (b) − f (a)

= (3)
x−a b−a

Now this line must intersect the x-axis (where y = 0), and we can approximate that as
the root of the given equation:
f (x) = 0

9
So we can write for the approximate root (x = ξ) :

f (a) · (b − a)
ξ−a=−
f (b) − f (a)

(b − a) (b − a)
⇒ ξ = a − f (a) = b − f (b) (4)
f (b) − f (a) f (b) − f (a)
Now we can write a = x0 and b = x1 and think that ξ = x2 is the next approximation of
the root. So we can write an iterative expression:,

x1 − x 0
x2 = x1 − f (x1 ) · (5)
f (x1 ) − f (x0 )

Thus, to evaluate the approximate root , we use the following general iterative for-
mula:
xn − xn−1
xn+1 = xn − f (xn ) · (6)
f (xn ) − f (xn−1 )

Figure 2: The successive points and secants are shown

import numpy as np
from [Link] import root_scalar

# Define the function

def f(x):
return x**3 - x - 2 # Example: solve f(x)=0

10
# Apply secant method with two initial guesses
sol = root_scalar(f, method=’secant’, x0=1, x1=2)

# Display results
if [Link]:
print("Root found:", [Link])
else:
print("Method did not converge")

Explanation of the Python Program (Secant Method using SciPy)

1. Importing Libraries: The program imports two libraries:

• numpy as np for numerical operations.

• root_scalar from [Link] which provides the secant method.

2. Defining the Function: The function is defined as

f (x) = x3 − x − 2

We want to find the root of this nonlinear equation.

3. Applying Secant Method: The function root_scalar is used with

• method=’secant’ to specify the secant method.

• Initial guesses x0 = 1 and x1 = 2.

This function iteratively updates the values until it finds the root.

4. Checking Convergence: After the iterations, the program checks whether the
solution has converged.

• If converged, it prints the root value.

• Otherwise, it prints that the method did not converge.

def f(x):
return x**3 - 2*x - 5

x0, x1, tol = eval(input(’x0, x1, tol \n’))

while abs(f(x1)) >= tol:

11
 x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
x0, x1 = x1, x2
print(’x1 = {}, f(x1) = {}’.format(x1, f(x1)))

print(’Root =’, x1)

1. Function Definition: The function f (x) = x3 − 2x − 5 is defined in Python using

def f(x): return x**3 - 2*x - 5. This is the equation whose root we want to
find using the secant method.

2. Input: The user provides three values: the initial guesses x0 , x1 , and the tolerance
value tol. This is done using: x0, x1, tol = eval(input(’x0, x1, tol \n’)).

3. Iteration Process: A while loop runs until the condition |f (x1 )| ≥ tol is satisfied.
This ensures the process continues until the approximate root is accurate within
the given tolerance.

4. Secant Formula Update: In each iteration, the next approximation x2 is calcu-

lated using the secant method formula:

(x1 − x0 )
x2 = x1 − f (x1 ) · .
f (x1 ) − f (x0 )

Then the values are updated: x0, x1 = x1, x2.

5. Intermediate Output: After each iteration, the program prints the current value
of x1 and f (x1 ) using: print(’x1={}, f(x1)={}’.format(x1, f(x1))).

6. Final Output: Once the loop ends, the program prints the root approximation
using: print(’Root =’, x1). This is the solution of the equation within the given
tolerance.

Note:

It can be observed that the secant method is faster than Bisection method. It does not
require that the root remain bracketed, like bisection method [Link] it does not
always converge. The regula falsi method (or the false position method) uses the same
formula as the secant method with s difference.

12
3.3. Newton-Raphson Method:

This method is a improved version of the other two methods we discussed so far. If x0 is
the approximate root of the equation:

f (x) = 0

we can expand the function in Taylor’s series around that point. Suppose, x1 = x0 + h
be the correct root of the equation, so that f (x1 ) = 0. Taylor series expansion:

h2 ′′
f (x0 + h) = 0 = f (x0 ) + h · f ′ (x0 ) +
f (x0 ) + · · · (7)
2!
Neglecting the terms with second and higher order derivatives , we have

f (x0 )
f (x0 ) + h · f ′ (x0 ) = 0 ⇒ h = − (8)
f ′ (x0 )

Thus an approximation of the root can be written as,

f (x0 )
x1 = x0 + h = x0 − (9)
f ′ (x0 )

Given an initial approximation x0 , we can generate x1 and then successively, x2 , x3 · · ·

in order to reach closer and closer to the root.
In general, we get the iterative expression:

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

This is the Newton-Raphson (RP) formula.

Algorithm 2 Newton-Raphson
1: Define: f (x), f ′ (x)
2: Give the starting point a and tolerance value
3: Calculate: ff′(x)
(x)
and update x until |f (x)| > tol

def f(x): return x**3 - 2*x - 5 # FUNCTION

def h(x): return 3*x**2 - 2 # Derivative

x, tol = eval(input(’initial point, tolerance \n’))

while abs(f(x)) >= tol:
x = x - f(x)/h(x)
print(x)

13
print(’Root =’, x)

Program Output:

initial point, tolerance

3.0,0.001
2.36
2.1271967801588163
2.095136036933634
2.0945516738242675
Root= 2.0945516738242675

alternate program

from [Link] import newton

# Define function
f = lambda x: x**3 - 2*x - 5

# Newton’s method (derivative calculated internally if not provided)

root = newton(f, x0=3.0, tol=1e-3) # initial guess=3.0, tolerance=0.001

print("Root =", root)

Program Output:

Root = 2.094555400381212

Explanation of the SciPy Program

1. Importing library: The command from [Link] import newton im-

ports the Newton method function from the SciPy optimization module.

2. Defining the function: The function

f (x) = x3 − 2x − 5

is defined using a lambda expression in Python as f = lambda x: x**3 - 2*x - 5.

3. Calling Newton’s method: The built-in method newton(f, x0=3.0, tol=1e-3)

applies the Newton–Raphson algorithm to find the root of f (x) = 0 starting from
the initial guess x0 = 3.0 with tolerance 10−3 .

14
 4. Output: Finally, the root of the equation is displayed using print("Root =", root).

4. Integration

Definite integral of a function with single variable:

Z b
I= f (x)dx (11)
a

we often fail to integrate a function analytically even if the function is simply looking.
Also, sometimes the function f (x) is not explicitly given, instead some tabular values
are provided. In such cases we need to take help of numerical methods and computer
programs to evaluate the required integral. A definite integral is the ’area under the

Figure 3: The Area under the graph

curve’, S.
For the numeric integration, we proceed by approximating the function, f (x) between
the two limits by a suitable interpolating function, ϕ(x).
Z b Z b
∴ f (x)dx ≈ ϕ(x)dx (12)
a a

We can approximate the function f (x) by the Newton-Gregory forward difference for-
mula.

t(t − 1) 2 t(t − 1)(t − 2) 3 t(t − 1)(t − 2) · · · (t − n + 1 n

ϕ(x) = y0 +t·∆y0 + ∆ y0 + ∆ y0 +· · ·+ ∆ y0
2! 3! n!
(13)

15
where, t = (x−x
h
0)

To compute the integral numerically, with the interpolating function with forward dif-
ferences, we may truncate the series at some point according to the level of accuracy we
may set and to the computer time we may want to devote.

4.1. Rectangular Rule: (n=0)

This is the simplest case. Here we have , h = b − a and f (x) ≈ ϕ(x) = y0

Z b
∴I≈ ϕ(x)dx = hy0 (14)
a

The value of the integral is the area of the rectangle with width h = b − a and height
y0 .

Composite Formula

The error of the integral can be reduced if the interval (the limit of the integration) ,
[a, b] is made smaller. This means if we make, h = (b − a) sufficiently small , we can
get better result even by this crude method of integration. It will be clear that for a
given value of h, the refined formula by Trapezoidal rule and Simpson’s rules give more
accurate results.

With any numerical scheme , we can achieve sufficient accuracy if we can make the value
of h smaller and smaller. To achieve better accuracy , we can devide the interval [a, b]
into n number of equal sub intervals such that,

(a = x0 ) < x1 < x2 < x3 · · · < (xn = b)

Thus we can write the definite integral in the composite form,

Z b Z x1 Z x2 Z x3 Z xn
I= f (x)dx = ydx + ydx + ydx + · · · + ydx (15)
a x0 x1 x2 xn−1

(b−a)
In the above equation, y = f (x) and width of each sub intervals, h = n
.

Using, the rectangular rule for integration over each sub interval, we get a composite
formula,:
Z xn
I= ydx = h[y0 + y1 + y2 + · · · + yn−1 ] (16)
x0

Error: The error for the approximation for the function is ϵ(x) = h · tf ′ (ξ) i,e for the
interpolation. The error for this integration wil be proportional to h2 .

16
4.2. Trapezoidal rule :(n=1)

Here, x0 = a, x1 = b and f (x) ≈ ϕ(x) = y0 + t · ∆y0 .

Z b Z b
I≈ ϕ(x) = [y0 + t(y1 − y0 )]hdt
a a
1
= hy0 + ht(y1 − y0 ) ·
2
y0 + y 1
 
=h (17)
2

The value of the integral is the area of the trapezium with base h = (b − a) and bounded
by the ordinates y0 and y1 .

Composite formula

As before, we can apply the trapezoidal rule over each of the sub-intervals and achieve
the composite formula sa the following:
Z xn
h
I− ydx = [(y0 + y1 ) + (y1 + y2 ) + · · · + (yn−1 + yn )]
x0 2

h
I= [y0 + 2(y1 + y2 + · · · + yn−1 ) + yn ] (18)
2

To evaluate the integral : Z 1

x2 dx
0

import numpy as np
x,h=[Link](0,1,101,retstep=True)
y=x**2
I=0.5*h*(y[0]+y[-1]+2*sum(y[1:-1]))
print (’Integral=’,I)

Program Output:

Integral= 0.33335000000000004

Explanation of the Program

1. Importing library: The command

import numpy as np

17
 imports the NumPy library, which provides numerical and mathematical tools for
array operations.

2. Creating grid points: The command

x,h = [Link](0,1,101,retstep=True)

creates 101 equally spaced points between 0 and 1.

• x stores the array of values.

• h stores the spacing between two consecutive points, i.e. step size.

3. Defining the function: The line

y = x**2

computes y = x2 for each element of x.

4. Applying Trapezoidal Rule: The line

I = 0.5*h*(y[0] + y[-1] + 2*sum(y[1:-1]))

applies the trapezoidal rule for numerical integration:

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

where h is the step size, f (x) is the function, x0 is the first point, and xn is the last
point.

5. Displaying the result: Finally,

print(’Integral=’, I)

prints the approximate value of the integral.

Alternate method of the above program:

import numpy as np
from [Link] import quad

f = lambda x: x**2
I, _ = quad(f, 0, 1)

18
print("Integral=", I)

Program Output:

Integral= 0.33333333333333337

Explanation of the Program

1. Importing libraries:
import numpy as np

imports the NumPy library (though in this code it is not directly used).

from [Link] import quad

imports the function quad from the [Link] module, which is used for
numerical integration.

2. Defining the function: The line

f = lambda x: x**2

defines the function f (x) = x2 using a lambda expression.

3. Performing integration: The command

I, _ = quad(f, 0, 1)

calls the quad function to evaluate

Z 1 Z 1
I= f (x) dx = x2 dx
0 0

The result is stored in I, while the second returned value (an estimate of error) is
ignored by assigning it to _.

4. Displaying the result: Finally,

print("Integral=", I)

prints the computed value of the integral. The expected result is

1
I= ≈ 0.3333
3

19
4.3. Simpson’s 1/3 rule (n=2):

For n = 2 , we consider up to ∆2 y0 term. So, we consider the three points : (x0 , y0 ),

(x1 , y1 ) and (x2 , y2 ).
a+b
If we consider x0 = a and x2 = b, we have x1 = 2
(the mid point).
Z b" #
Z b
t(t − 1) 2
∴I≈ ϕ(x)dx = y0 + t · ∆y0 + ∆ y0 dx
a a 2!
Z 2" #
(t2 − t)
= y0 + t · (y1 − y0 ) + (y2 − 2y1 + y0 ) hdt
0 2
h
= [y0 + 4y1 + y2 ] (19)
3

Note: The above integral is basically the area under the quadratic curve that is passing
through the three points (x0 , y0 ), (x1 , y1 ) and (x2 , y2 ).The quadratic function is due to
the consideration of the interpolating polynomial up to the 2nd degree.

Independent derivation of Simpson’s formula:

As a mathematical exercise , we may also derive the Simpson’s 1/3 rule by considering a
quadratic function,
y = f (x) = a0 + a1 x + a2 x2
 
and choosing the three points (a, f (a)), ( a+b
2
,f a+b
2
) and (b, f (b)) to evaluate the pa-
 
a+b
rameters a0 , a1 and a2 . The parameters are expressed in the terms of f (a), f 2
and
f (b).
Z b Z b 
∴I= f (x)dx ≈ a0 + a1 x + a2 x2 dx
a a
b2 − a 2 b3 − a 3
= a0 (b − a) + a1 + a2 (20)
2 3

Putting the values of a0 , a1 and a2 we find out that,

" ! #
b−a a+b
I= f (a) + 4f + f (b) (21)
6 2

b−a
Considering h = 2
, the formula takes the form as derived before.

Simpson’s rule : Composite formula

For Simpson’s 1/3 rule, to be applied between two points , we require two equally
spaced sub intervals each of length h. So the rule requires the division of the whole range
[a, b] into an even number of sub intervals with width h.

20
To be explicit, we have two sub intervals given by the three points (x0 , x1 , x2 ) over which
we apply Simpson’s rule.
Z x2
h
ydx == (y0 + 4y1 + y2 )
x0 3

Similarly, for the next two sub intervals given by,(x2 , x3 , x4 ) , we get
Z x4
h
ydx = (y2 + 4y3 + y4 )
x2 3

· · · and so on,
Z xn
h
∴I= ydx = [(y0 + 4y1 + y2 ) + (y2 + 4y3 + y4 ) + · · · + (yn−2 + 4yn−1 + yn )]
x0 3
h
= [y0 + 4 (y1 + y3 + · · · + yn−1 ) + 2 (y2 + y4 + · · · + yn−2 )] (22)
3

Python script : Composite Simpson’s formula (Using Numpy array Slic-

ing)

import numpy as np
x,h=[Link](0,1,101,retstep=True)
y=x**2
I=h/3*(y[0]+y[-1]+4*sum(y[1:-1:2])+2*sum(y[2:-2:2]))
print (’Integral=’,I)

Program Output:

Integral= 0.33333333333333337

Alternate program for Python script : Composite Simpson’s formula:

import numpy as np
from [Link] import simpson

# define grid points

x = [Link](0, 1, 101)

# function values
y = x**2

# compute integral using Simpson’s rule

I = simpson(y, x)

21
print("Integral=", I)

Program Output:

Integral= 0.3333333333333333

Step-by-step explanation of the program

1. Import libraries

import numpy as np
from [Link] import simpson

• numpy (as np) provides array creation and vectorized operations.

• simpson is SciPy’s implementation of the composite Simpson’s rule for numer-

ical integration.

2. Define grid points

x = [Link](0, 1, 101)

This creates N = 101 equally spaced sample points on the interval [0, 1]. The
uniform step size (spacing) is

xN −1 − x0 1−0
h= = = 0.01.
N −1 100

3. Evaluate the function at the grid points

y = x**2

This computes yi = f (xi ) = x2i for every sample point xi .

4. Integrate using Simpson’s rule

I = simpson(y, x)

The composite Simpson’s 1/3 rule (for equally spaced points) approximates the
integral by
h
 X X 
I≈ y0 + yN −1 + 4 yi + 2 yi ,
3 i odd i even, 1≤i≤N −2

22
 where the sums run over interior indices: the 4-coefficient multiplies values at odd
indices and the 2-coefficient multiplies values at even interior indices. Note: Simp-
son’s rule requires an even number of intervals (i.e. an odd number of sample points
N ). In this code N = 101 (odd), so there are 100 intervals (even) — the rule applies
directly.

5. Print the result

print("Integral=", I)

This prints the computed numerical value. For f (x) = x2 on [0, 1] the exact integral
is Z 1
1
x2 dx = ≈ 0.333333 . . . ,
0 3
and the Simpson result with 101 points will match this to numerical precision.

4.4. Simpson’s 3/8 rule: (n=3)

We can go further Retaining one more term in the interpolating formula we try even
better approximation,
Z bZ b" #
Z b
t(t − 1) 2 t(t − 1)(t − 2) 3
I≈ ϕ(x)dx = y0 + t · ∆y0 + ∆ y0 + ∆ y0 dx (23)
a a a 2! 3!

It can be shown that,

Z x3
3h
I= ydx = [y0 + 3y1 + 3y2 + y3 ] (24)
x0 8

This is based on the cubic interpolation (rather than on the quadratic interpolation as
on Simpson’s 1/3 rule). So the 3rd order polynomial goes through the four points :
x0 , x1 , x2 , x3 .

The composite rule of Simpson’s 3/8 formula becomes,

Z xn
I= ydx
x0
3h
= [y0 + 3 (y1 + y2 + y4 + y5 + · · · + yn−1 ) + 2 (y3 + y6 + y9 + · · · + yn−1 ) + yn ]
8
(25)

Write the python script of the Simpson’s 3/8 formula:

import numpy as np

23
def simpsons_38(f, a, b, n):
if n % 3 != 0: raise ValueError("n must be a multiple of 3")
h = (b - a) / n
x = [Link](a, b, n+1)
y = f(x)
return (3*h/8) * (y[0] + y[-1] +
3*[Link](y[1:-1][[Link](1,n) % 3 != 0]) +
2*[Link](y[3:-1:3]))

# Example
f = lambda x: x**2
print(simpsons_38(f, 0, 1, 6))

Program Output:

0.3333333333333333

Alternative python script of the above program:

import sympy as sp

x = [Link](’x’)
f = x**2 # Example function

I = [Link](f, (x, 0, 1))

print("Integral =", [Link]())

Program output:

Integral = 0.333333333333333

5. Ordinary Differential Equation(ODE):

Let us consider a general the first order differential equation,

dy
= f (x, y) (26)
dx

With the initial condition: y (x0 ) = y0 .

In the following , we will describe the methods in which the step-by-step solutions are

24
obtained at different points staring from a given initial value. To begin with, the first order
differential equation is either integrated by an approximation method or the function is
approximated by the Taylor’s series expansion.

5.1. Euler’s Method:

We want to solve the differential equation (26) for y at different values of x. Consider
these points are: xn = x0 + n · h, n = 1, 2, 3 · · · .

Between the first two points: [x0 , x1 ] we have,

Z x1
y1 = y0 + f (x, y)dx (27)
x0

If we assume , f (x, y) = f (x0 , y0 ), we can get an approximate formula:

y1 = y0 + h · f (x0 , y0 ) (28)

, where h = x1 − x0 . Note that the integration is done here by the rectangle rule.

Similarly, for the next range, [x1 , x2 ]:

y2 = y1 + h · f (x1 , y1 ) (29)

So, we can write the general formula which is called Euler’s formula:

yn+1 = yn + h · f (xn , yn ) (30)

Note: The formula could also be obtained by retaining up to first order term in the
Taylor series expansion:
h2
y(x) = y0 + hy0′ + y0′′ + · · ·
2!

About Error:

Note that the Euler’s formula is a crude method. To achieve a good accuracy, the step size
h must be very small. However, for increasing number of grid points n, the truncation
error may also increase. In that case, there is a possibility that the method becomes
unstable. So, we must be careful that this method may not work well for any arbitrarily
defined function f , especially for the non-linear differential equations.

In the Euler’s method, the error is clearly φ(h2 ). It is customery in numerical solution
of differential equations ; we call a method of order n if the error is φ(hn+1 ). Thus Euler’s
method is a first order method.

25
Algorithm 3 Euler method
1: Define: The function f (x, y)
2: Set the initial value y0 and step size h.
3: Update y = y + h ∗ f (x, y) and x = x + h. Iterate this step in a loop.

dy
Problem: dx
= 2xy = f (x, y), with initial condition y(−2) = 1

import numpy as np
import [Link] as plt

# dy/dx = f(x,y)
f = lambda x, y: 2*x*y

# Euler’s method in compact form

def euler(f, x0, y0, h, n):
x = [Link](x0, x0 + n*h, n+1)
y = [Link](n+1)
y[0] = y0
for i in range(n):
y[i+1] = y[i] + h*f(x[i], y[i])
return x, y

# Initial condition
x, y = euler(f, -2, 1, 0.1, 40)

# Exact solution
y_exact = [Link](x**2 - 4)

# Print few values

print(np.column_stack((x[:6], y[:6], y_exact[:6])))

# Plot
[Link](x, y, ’ro-’, label="Euler’s method")
[Link](x, y_exact, ’b-’, label="Exact solution")
[Link](), [Link](), [Link]()

Program output

[[-2. 1. 1. ]
[-1.9 0.6 0.67705687]
[-1.8 0.372 0.46766643]

26
 [-1.7 0.23808 0.32955896]
[-1.6 0.1571328 0.23692776]
[-1.5 0.1068503 0.17377394]]

Explanation of the Program

Step 1: Import Libraries

We import the necessary Python libraries:

• numpy for numerical computations and array handling.

• [Link] for plotting graphs.

Step 2: Define the Differential Equation

The given first-order ordinary differential equation is

dy
= f (x, y) = 2xy.
dx

Step 3: Euler’s Method

Euler’s method is used to approximate the solution of the differential equation. The
iterative formula is:
yi+1 = yi + h f (xi , yi ),

where

27
 • h = step size,

• xi = the current value of x,

• yi = the current value of y.

In Python, this is implemented by initializing arrays for x and y, and updating them in
a loop.

Step 4: Initial Conditions

We are given the initial condition

y(−2) = 1.

Thus, the starting point is (x0 , y0 ) = (−2, 1). The program computes n = 40 steps with
step size h = 0.1.

Step 5: Exact Solution

The exact solution of the differential equation is obtained by solving the ODE analyti-
cally:
dy dy
= 2xy ⇒ = 2x dx.
dx y
Integrating both sides:
ln y = x2 + C.

Applying the initial condition y(−2) = 1, we find

ln(1) = (−2)2 + C ⇒ C = −4.

Hence, the exact solution is

2 −4
y(x) = ex .

Step 6: Output of the Program

The program prints the first six values of:

• x (grid points),

• y (Euler’s approximate solution),

• yexact (exact solution).

Step 7: Graphical Comparison

Finally, the program plots:

28
 • The numerical solution using Euler’s method (red circles with connecting lines).

• The exact solution (blue continuous curve).

This allows us to visually compare the approximation with the analytical solution.

5.2. Modified Euler Method:

We consider the integral: y1 = y0 + xx01 f (x, y)dx. In Euler’s method , we approximated

R

this by Rectangular rule and obtained the expression : y1 ≈ y0 + hf (x0 , y0 ).

Now we can do a better approximation of the integral by Trapezoidal rule.

h
y1 = y0 + [f (x0 , y0 ) + f (x1 , y1 )] (31)
2

The equation (31) can be solved iteratively. Given (x0 , y0 ) we may begin with an approx-
imate value of y1 and put it on the right to obtain a new value of y1 on the left. This new
y1 is put back on the right again and obtain a better approximation and so on. Thus ,
we solve for y1 in a self-consistent way.
(0)
Step-1: We predict y1 by Euler’s formula:

(0)
y1 = y0 + hf (x0 , y0 ) (32)

Step-2: We correct it by:

(1) hh (0)
i
y1 = y0 + f (x0 , y0 ) + f (x1 , y1 ) (33)
2
(1)
Putting (32) into (33) , we get y1 . This is also called predictor-corrector method.
(1) (2)
Next we use y1 to get the next approximate value y1 and so on. We can carry on this
according to the level of accuracy , we set.

The general iteration formula:

(n) hh (n−1)
i
y1 = y0 + f (x0 , y0 ) + f (x1 , y1 ) (34)
2

5.3. Runge-Kutta Methods: Runge-Kutta 2ND Order (RK2):

We started by plugging Euler’s original formula, y1 = y0 + hf (x0 , y0 ) into the modified

formula and obtained:

h
y1 = y0 + [f0 + f (x0 + h, y0 + hf0 )] (35)
2

29
where, f0 = f (x0 , y0 ).

Now we designate two steps, k1 = hf0 and k2 = hf (x0 + h, y0 + k1 ).Thus we can write
the following steps:
k1 = hf0
k2 = hf (x0 + h, y0 + k1 )
1
y1 = y0 + [k1 + k2 ]
2

python script for RK2

import numpy as np
import [Link] as plt

# Differential equation dy/dx = f(x,y)

def f(x, y):
return 2*x*y # Example: dy/dx = 2xy

# Runge-Kutta 2nd order method

def rk2(f, x0, y0, h, n):
x = [Link](n+1)
y = [Link](n+1)
x[0], y[0] = x0, y0

for i in range(n):
k1 = h*f(x[i], y[i])
k2 = h*f(x[i] + h, y[i] + k1)
y[i+1] = y[i] + 0.5*(k1 + k2)
x[i+1] = x[i] + h
return x, y

# Initial conditions and parameters

x0, y0 = -2, 1 # y(-2) = 1
h, n = 0.1, 40 # step size and number of steps

# Apply RK2
x, y = rk2(f, x0, y0, h, n)

# Exact solution: y(x) = exp(x^2 - 4)

y_exact = [Link](x**2 - 4)

30
# Print first few values
print(" x RK2 approx Exact")
print(np.column_stack((x[:6], y[:6], y_exact[:6])))

# Plot results
[Link](x, y, ’ro-’, label="RK2 Method")
[Link](x, y_exact, ’b-’, label="Exact Solution")
[Link]("x"); [Link]("y")
[Link](); [Link](); [Link]()

Program output:

x RK2 approx Exact

[[-2. 1. 1. ]
[-1.9 0.686 0.67705687]
[-1.8 0.4791024 0.46766643]
[-1.7 0.34073763 0.32955896]
[-1.6 0.24683034 0.23692776]
[-1.5 0.18216079 0.17377394]]

Explanation of the Python Program

1. Importing Libraries

import numpy as np

31
import [Link] as plt

We import the numpy library for numerical computations and [Link] for
plotting graphs.

2. Defining the Differential Equation

def f(x, y):

return 2*x*y

This defines the differential equation

dy
= f (x, y) = 2xy.
dx

3. Runge-Kutta 2nd Order Method

def rk2(f, x0, y0, h, n):

x = [Link](n+1)
y = [Link](n+1)
x[0], y[0] = x0, y0

We create arrays for x and y values, initialize with initial condition (x0 , y0 ).

for i in range(n):
k1 = h*f(x[i], y[i])
k2 = h*f(x[i] + h, y[i] + k1)
y[i+1] = y[i] + 0.5*(k1 + k2)
x[i+1] = x[i] + h
return x, y

The RK2 iteration uses:

k1 = hf (xi , yi ), k2 = hf (xi + h, yi + k1 ),

yi+1 = yi + 21 (k1 + k2 ), xi+1 = xi + h.

4. Initial Conditions and Parameters

x0, y0 = -2, 1 # y(-2) = 1

h, n = 0.1, 40 # step size and number of steps

We start from x = −2, y(−2) = 1, step size h = 0.1, and number of steps n = 40.

32
5. Applying the RK2 Function

x, y = rk2(f, x0, y0, h, n)

This computes approximate values of y(x) using RK2.

6. Exact Solution

y_exact = [Link](x**2 - 4)

The analytical solution is

2 −4
y(x) = ex .

7. Displaying Results

print(np.column_stack((x[:6], y[:6], y_exact[:6])))

This prints the first six values of x, approximate y, and exact y(x) for comparison.

8. Plotting the Solution

[Link](x, y, ’ro-’, label="RK2 Method")

[Link](x, y_exact, ’b-’, label="Exact Solution")
[Link]("x"); [Link]("y")
[Link](); [Link]();
[Link](’[Link]’)
[Link]()

The RK2 approximate solution is plotted with red circles and compared to the exact
solution shown as a blue curve. The figure is also saved as [Link].

5.4. Runge-Kutta 4th order (RK4):

The formula for the higher order Runge-Kutta method can be derived. In the following,
we give an outlined how the 4th order Runge-Kutta formula can be obtained.

The steps of RK4 method can be written in a general ways,

k1 = hf0
k2 = hf (x0 + αh, y0 + βk1 )
y1 = y0 + Ak1 + Bk2

To begin with the parameters α, β and the weight factors A and B are taken as arbitrary.
Basically, the second slope is constructed at an arbitrary point away from (x0 , y0 ) .

33
The general steps of 4th order RK4 METHOD :

k1 = hf (x0 , y0 )
k2 = hf (x0 + α1 h, y0 + β1 k1 )
k3 = hf (x0 + α2 h, y0 + β2 k1 + β3 k2 )
k4 = hf (x0 + α3 h, y0 + β4 k1 + β5 k2 + β6 k3 )
y1 = y0 + Ak1 + Bk2 + Ck3 + Dk4

The parameters are determined through the expanding both sides of the last equation by
Taylor series and then comparing the coefficients of the similar terms on the both sides
as done for the RK2.

The choice of the parameter values can be arbitrary and so we can have different forms
of the RK4 formula. The commonly used RK4 formula are given:

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

python script:

import numpy as np
from [Link] import solve_ivp
import [Link] as plt

# Example ODE: dy/dx = 2*x*y

def f(x, y):
return 2 * x * y

# Initial condition and domain

x0, y0 = -2, 1
xf = 2

# Solve using RK45 (which is adaptive but based on Runge-Kutta 4/5)

sol = solve_ivp(f, [x0, xf], [y0], method=’RK45’, dense_output=True)

34
# Generate points for plotting
x_vals = [Link](x0, xf, 100)
y_vals = [Link](x_vals)[0]

# Plot solution
[Link](x_vals, y_vals, label="RK4 Solution")
[Link]("x")
[Link]("y")
[Link]()
[Link]()
[Link](’[Link]’)
[Link]()

Program output:

6. Second order Differential Equation Solution:

To solve a a second order differential equation in computer,We split it into two first
order differential equations.
d2 y dy
2
+λ + ky = 0 (36)
dx dx

35
split into two coupled first order equations:

dy
= z = f1 (x, y, z)
dx
dz
= −λz − ky = f2 (x, y, z)
dx

Python script:

import numpy as np
import [Link] as plt

# Parameters
lam, k = 1.0, 2.0
h = 0.1
x0, xf = 0, 10

# System of ODEs
def f(x, Y):
y1, y2 = Y
return [Link]([y2, -lam*y2 - k*y1])

# RK4 loop
x_vals = [Link](x0, xf+h, h)
Y = [Link]((len(x_vals), 2))
Y[0] = [1, 0] # initial conditions y(0)=1, y’(0)=0

for i in range(len(x_vals)-1):
k1 = h*f(x_vals[i], Y[i])
k2 = h*f(x_vals[i] + h/2, Y[i] + k1/2)
k3 = h*f(x_vals[i] + h/2, Y[i] + k2/2)
k4 = h*f(x_vals[i] + h, Y[i] + k3)
Y[i+1] = Y[i] + (k1 + 2*k2 + 2*k3 + k4)/6

# Plot solution
[Link](x_vals, Y[:,0], label="RK4 Solution y(x)")
[Link]("x")
[Link]("y")
[Link]()
[Link]()
[Link](’[Link]’)

36
[Link]()

Program Output:

Explanation of the RK4 Program

We are solving the second order differential equation

d2 y dy
2
+ λ + ky = 0
dx dx

using the Runge–Kutta 4 (RK4) method.

1. Reduction to first-order system: Define two variables:

dy
y1 = y, y2 =
dx

Then the system becomes:

dy1 dy2
= y2 , = −λy2 − ky1
dx dx

2. Define parameters and initial conditions: We choose λ = 1.0, k = 2.0, step

size h = 0.1, and domain x ∈ [0, 10]. The initial conditions are:

y(0) = 1, y ′ (0) = 0

37
 3. Define the system in Python: A function f(x,Y) is written that returns
 
y2
f (x, Y ) =  
−λy2 − ky1

4. Apply the RK4 method: For each step, compute

k1 = hf (xn , Yn )
 
k2 = hf xn + h2 , Yn + k1
2
 
k3 = hf xn + h2 , Yn + k2
2

k4 = hf (xn + h, Yn + k3 )

Then update:
1
Yn+1 = Yn + (k1 + 2k2 + 2k3 + k4 )
6

5. Iterate over the grid: Store all computed values of y(x) and proceed step by
step across the interval.

6. Plot the solution: After computing, we plot y(x) using matplotlib.

alternate python script:

import numpy as np
from [Link] import odeint
import [Link] as plt

# Parameters
lam, k = 1.0, 2.0

# System of ODEs
def system(Y, x):
y1, y2 = Y
return [y2, -lam*y2 - k*y1]

# Initial conditions and domain

Y0 = [1, 0] # y(0)=1, y’(0)=0
x = [Link](0, 10, 200)

# Solve using odeint

sol = odeint(system, Y0, x)

38
# Plot solution
[Link](x, sol[:,0], label="y(x)")
[Link]("x"); [Link]("y")
[Link](); [Link]();
[Link](’[Link]’)
[Link]()

Explanation of the odeint() Program

We are solving the second order differential equation

d2 y dy
2
+ λ + ky = 0
dx dx

using the odeint() solver from the [Link] module.

1. Reduction to a system of first-order equations: Introduce new variables:

dy
y1 = y, y2 =
dx

Then the system becomes:

dy1 dy2
= y2 , = −λy2 − ky1
dx dx

2. Define parameters: We choose λ = 1.0 and k = 2.0 as constants.

3. Write the system in Python: In the program, the function

system(Y, x) → [y2 , −λy2 − ky1 ]

encodes the two coupled equations.

4. Set initial conditions: At x = 0, we specify

y(0) = 1, y ′ (0) = 0

which correspond to Y (0) = [1, 0].

5. Choose the domain of solution: The solution is computed for x ∈ [0, 10] using
200 equally spaced points.

39
 6. Call odeint(): The solver integrates the system over the chosen domain:

sol = odeint(system, Y0, x)

Here, sol[:,0] contains the values of y(x) and sol[:,1] contains y ′ (x).

7. Plot the result: Finally, the solution y(x) is plotted against x using the matplotlib
library.

Q1: Solve the following coupled differential equations: write the python script of the
following problem,
dx x2 dy
=y+x− , = −x
dt 3 dx
for the following different initial conditions: x(0) = 0, y(0) = −1, −2, −3, −4. Obtain
(x, y) plot to see. use scipy/ sympy module with optimized steps of the programs,

# ---------------- Part 1: Solve the differential equations ----------------

import numpy as np
from [Link] import solve_ivp
import [Link] as plt

# Define the system of ODEs

def system(t, z):
x, y = z
dxdt = y + x - x**2 / 3
dydt = -x * dxdt
return [dxdt, dydt]

# Time span for the solution

t_span = (0, 10)
t_eval = [Link](t_span[0], t_span[1], 500)

# Initial conditions
initial_ys = [-1, -2, -3, -4]

# Dictionary to store solutions

solutions = {}

# Solve the system for each initial condition

for y0 in initial_ys:

40
 sol = solve_ivp(system, t_span, [0, y0], t_eval=t_eval, method=’RK45’)
solutions[y0] = sol

# Print the solutions

print(f"\nSolution for initial condition y0 = {y0}:")
print("t x(t) y(t)")
for t, x, y in zip(sol.t, sol.y[0], sol.y[1]):
print(f"{t:.2f} {x:.5f} {y:.5f}")

# ---------------- Part 2: Plot the solutions ----------------

[Link](figsize=(8,6))
colors = [’b’, ’g’, ’r’, ’m’]

for y0, color in zip(initial_ys, colors):

sol = solutions[y0]
[Link](sol.y[0], sol.y[1], color=color, label=f’y0={y0}’)

[Link](’x(t)’)
[Link](’y(t)’)
[Link](’Phase plot of (x, y) for different initial conditions’)
[Link]()
[Link](True)
[Link](’[Link]’)
[Link]()

Program Output:

Solution for initial condition y0 = -1:

t x(t) y(t)
0.00 0.00000 -1.00000
0.02 -0.02024 -1.00020
0.04 -0.04091 -1.00084
0.06 -0.06203 -1.00192
0.08 -0.08362 -1.00350
0.10 -0.10570 -1.00559
0.12 -0.12832 -1.00825
0.14 -0.15153 -1.01161
0.16 -0.17534 -1.01569
0.18 -0.19979 -1.02047
0.20 -0.22488 -1.02596

41
0.22 -0.25066 -1.03221
0.24 -0.27716 -1.03925
0.26 -0.30442 -1.04715
0.28 -0.33249 -1.05600
0.30 -0.36144 -1.06588
0.32 -0.39131 -1.07691
0.34 -0.42218 -1.08924
0.36 -0.45413 -1.10299
0.38 -0.48724 -1.11835
0.40 -0.52158 -1.13548
0.42 -0.55726 -1.15460
0.44 -0.59437 -1.17591
0.46 -0.63303 -1.19965
0.48 -0.67333 -1.22607
0.50 -0.71541 -1.25544
0.52 -0.75938 -1.28804
0.54 -0.80537 -1.32416
0.56 -0.85354 -1.36419
0.58 -0.90418 -1.40924
0.60 -0.95747 -1.45957
0.62 -1.01356 -1.51538
0.64 -1.07266 -1.57719
0.66 -1.13507 -1.64579
0.68 -1.20113 -1.72228
0.70 -1.27125 -1.80808
0.72 -1.34592 -1.90489
0.74 -1.42567 -2.01472
0.76 -1.51111 -2.13988
0.78 -1.60291 -2.28298
0.80 -1.70181 -2.44693
0.82 -1.80860 -2.63493
0.84 -1.92432 -2.85215
0.86 -2.05067 -3.10836
0.88 -2.18847 -3.40418
0.90 -2.33912 -3.74464
0.92 -2.50480 -4.14109
0.94 -2.68846 -4.61118
0.96 -2.89382 -5.17886
0.98 -3.12535 -5.87441

42
1.00 -3.38833 -6.73440
1.02 -3.68878 -7.80169
1.04 -4.03555 -9.15173
1.06 -4.43977 -10.86849
1.08 -4.91673 -13.08348
1.10 -5.49135 -16.06175
1.12 -6.19822 -20.20246
1.14 -7.08764 -26.17566
1.16 -8.23277 -34.90586
1.18 -9.77387 -48.69250
1.20 -11.97162 -72.69638
1.22 -15.32519 -118.33955
1.24 -21.12847 -224.71975
1.26 -33.50851 -562.51563
1.28 -78.65122 -3089.96824

Solution for initial condition y0 = -2:

t x(t) y(t)
0.00 0.00000 -2.00000
0.02 -0.04049 -2.00082
0.04 -0.08186 -2.00335
0.06 -0.12419 -2.00772
0.08 -0.16759 -2.01425
0.10 -0.21216 -2.02302
0.12 -0.25794 -2.03407
0.14 -0.30502 -2.04751
0.16 -0.35351 -2.06351
0.18 -0.40353 -2.08231
0.20 -0.45522 -2.10424
0.22 -0.50874 -2.12967
0.24 -0.56428 -2.15907
0.26 -0.62203 -2.19297
0.28 -0.68222 -2.23196
0.30 -0.74509 -2.27670
0.32 -0.81090 -2.32795
0.34 -0.87992 -2.38650
0.36 -0.95246 -2.45324
0.38 -1.02884 -2.52911
0.40 -1.10950 -2.61584

43
0.42 -1.19506 -2.71592
0.44 -1.28584 -2.82977
0.46 -1.38235 -2.95876
0.48 -1.48528 -3.10545
0.50 -1.59551 -3.27354
0.52 -1.71414 -3.46789
0.54 -1.84242 -3.69455
0.56 -1.98184 -3.96069
0.58 -2.13406 -4.27469
0.60 -2.30092 -4.64604
0.62 -2.48468 -5.08767
0.64 -2.68907 -5.62611
0.66 -2.91655 -6.26816
0.68 -3.17126 -7.03705
0.70 -3.45971 -7.98062
0.72 -3.79079 -9.17129
0.74 -4.17574 -10.70609
0.76 -4.62820 -12.70663
0.78 -5.16695 -15.36440
0.80 -5.81895 -18.94637
0.82 -6.62547 -23.92845
0.84 -7.65712 -31.30008
0.86 -9.02181 -42.78546
0.88 -10.89826 -61.35626
0.90 -13.67571 -95.46542
0.92 -18.19151 -167.46418
0.94 -26.84913 -363.24331
0.96 -50.11909 -1258.39422
0.98 -331.41945 -54956.41522

Solution for initial condition y0 = -3:

t x(t) y(t)
0.00 0.00000 -3.00000
0.02 -0.06075 -3.00185
0.04 -0.12285 -3.00755
0.06 -0.18653 -3.01768
0.08 -0.25194 -3.03246
0.10 -0.31922 -3.05201
0.12 -0.38854 -3.07666

44
0.14 -0.46013 -3.10689
0.16 -0.53426 -3.14336
0.18 -0.61122 -3.18692
0.20 -0.69136 -3.23859
0.22 -0.77507 -3.29955
0.24 -0.86277 -3.37118
0.26 -0.95493 -3.45503
0.28 -1.05205 -3.55281
0.30 -1.15470 -3.66643
0.32 -1.26351 -3.79841
0.34 -1.37956 -3.95406
0.36 -1.50348 -4.13471
0.38 -1.63614 -4.34313
0.40 -1.77880 -4.58485
0.42 -1.93313 -4.86821
0.44 -2.10119 -5.20433
0.46 -2.28543 -5.60711
0.48 -2.48870 -6.09325
0.50 -2.71427 -6.68223
0.52 -2.96624 -7.40223
0.54 -3.25081 -8.30108
0.56 -3.57304 -9.40041
0.58 -3.94216 -10.77072
0.60 -4.37227 -12.54183
0.62 -4.88240 -14.90281
0.64 -5.49648 -18.10205
0.66 -6.24971 -22.56072
0.68 -7.19347 -28.87324
0.70 -8.42075 -38.42041
0.72 -10.08933 -53.94886
0.74 -12.47499 -80.80560
0.76 -16.19959 -134.16795
0.78 -22.82440 -263.24292
0.80 -37.97313 -723.13532
0.82 -108.07264 -5835.21173

Solution for initial condition y0 = -4:

t x(t) y(t)
0.00 0.00000 -4.00000

45
0.02 -0.08101 -4.00328
0.04 -0.16390 -4.01357
0.06 -0.24901 -4.03165
0.08 -0.33656 -4.05776
0.10 -0.42684 -4.09239
0.12 -0.52019 -4.13638
0.14 -0.61704 -4.19092
0.16 -0.71788 -4.25755
0.18 -0.82329 -4.33817
0.20 -0.93389 -4.43499
0.22 -1.05038 -4.55062
0.24 -1.17356 -4.68798
0.26 -1.30427 -4.85035
0.28 -1.44360 -5.04295
0.30 -1.59306 -5.27342
0.32 -1.75364 -5.54399
0.34 -1.92693 -5.86145
0.36 -2.11518 -6.23792
0.38 -2.32134 -6.69085
0.40 -2.54904 -7.24302
0.42 -2.80256 -7.92255
0.44 -3.08688 -8.76289
0.46 -3.40863 -9.81676
0.48 -3.77690 -11.15597
0.50 -4.20099 -12.83734
0.52 -4.69807 -15.02314
0.54 -5.29372 -17.98936
0.56 -6.02195 -22.12576
0.58 -6.93115 -28.05846
0.60 -8.09715 -36.78011
0.62 -9.66136 -50.63184
0.64 -11.88072 -74.72343
0.66 -15.25026 -120.12653
0.68 -21.07599 -226.41876
0.70 -33.49437 -565.84796
0.72 -78.72360 -3100.41543

46
Q2: Solve the following 2nd order non-homogeneous equation with RK4:

y ′′ − 2y ′ − 3y = 3x2 , y(0) = 9, y ′ (0) = −4

import numpy as np
import [Link] as plt

# Define the system

def f(x, Y):
y1, y2 = Y
dy1dx = y2
dy2dx = 2*y2 + 3*y1 + 3*x**2
return [Link]([dy1dx, dy2dx])

# RK4 step
def rk4_step(f, x, Y, h):
k1 = h * f(x, Y)

47
 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)
return Y + (k1 + 2*k2 + 2*k3 + k4)/6

# Parameters
x0, x_end = 0, 5
h = 0.1 # step size (minimal steps)
n = int((x_end - x0)/h)

# Initial conditions
Y = [Link]([9, -4]) # [y(0), y’(0)]
x_values = [x0]
y_values = [Y[0]]

# RK4 loop
x = x0
for i in range(n):
Y = rk4_step(f, x, Y, h)
x += h
x_values.append(x)
y_values.append(Y[0])

# Plot
[Link](x_values, y_values, ’b-’, label=’y(x) using RK4’)
[Link](’x’)
[Link](’y’)
[Link](’Solution of y\’\’ - 2y\’ - 3y = 3x^2 using RK4’)
[Link](True)
[Link]()
[Link](’[Link]’)
[Link]()

Program Output:

48
7. Curve fitting and Regression Analysis:

Q1: you have got the following data in your experiment: If you now want to do a least

x 9.56 10.62 11.45 12.01 13.89 14.23 17.5

y 7.05 6.88 6.75 6.66 6.38 6.33 5.82

square fit of the data with y = mx + c , then find the values of m and c.

import numpy as np
import [Link] as plt
from scipy import stats

# Data from the table

x = [Link]([9.56, 10.62, 11.45, 12.01, 13.89, 14.23, 17.5])
y = [Link]([7.05, 6.88, 6.75, 6.66, 6.38, 6.33, 5.82])

# Linear regression using [Link]

slope, intercept, r_value, p_value, std_err = [Link](x, y)

print(f"Slope (m) = {slope:.4f}")

print(f"Intercept (c) = {intercept:.4f}")
print(f"R-squared = {r_value**2:.4f}")

49
# Optional: Plotting
[Link](x, y, color=’b’, label=’Data points’)
[Link](x, slope*x + intercept, ’r-’, label=f’Fit: y = {slope:.2f}x + {intercept:.2f}
[Link](’x’)
[Link](’y’)
[Link](’Linear Regression Fit’)
[Link]()
[Link](True)
[Link](’[Link]’)
[Link]()

Program Output:

Slope (m) = -0.1541

Intercept (c) = 8.5177
R-squared = 0.9999

50