UNIVERSITY OF ENGINEERING AND TECHNOLOGY, TAXILA
FACULTY OF TELECOMMUNICATION AND INFORMATION ENGINEERING
University of COMPUTER ENGINEERING DEPARTMENT
ENGINEERING & TECHNOLOGY
TAXILA
A PROGRAMMER'S PERSPECTIVE
NUMERICAL METHODS
& PROBABILITY
ASSIGNMENT
03
MINAL FATIMA
21-CP-55
UNIVERSITY OF ENGINEERING AND TECHNOLOGY, TAXILA
FACULTY OF TELECOMMUNICATION AND INFORMATION ENGINEERING
COMPUTER ENGINEERING DEPARTMENT
CODE:
import numpy as np
import matplotlib
matplotlib.use('TkAgg') # or any other backend that works for you
import matplotlib.pyplot as plt
from scipy.integrate import quad
def sor_solver(A, b, x0=None, omega=1.25, tol=1e-6, max_iter=1000):
n = len(b)
x = x0 if x0 is not None else np.zeros(n)
for k in range(max_iter):
x_old = np.copy(x)
for i in range(n):
sigma = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i + 1:], x_old[i + 1:])
x[i] = (1 - omega) * x_old[i] + (omega / A[i, i]) * (b[i] - sigma)
if np.linalg.norm(x - x_old, ord=np.inf) < tol:
return x, k + 1
raise ValueError("SOR method did not converge")
def gauss_seidel_solver(A, b, x0=None, tol=1e-6, max_iter=1000):
n = len(b)
x = x0 if x0 is not None else np.zeros(n)
for k in range(max_iter):
x_old = np.copy(x)
for i in range(n):
sigma = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i + 1:], x_old[i + 1:])
x[i] = (b[i] - sigma) / A[i, i]
if np.linalg.norm(x - x_old, ord=np.inf) < tol:
return x, k + 1
raise ValueError("Gauss-Seidel method did not converge")
def trapezoidal_rule(func, a, b, n):
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = func(x)
result = h * (y[0] + 2 * np.sum(y[1:-1]) + y[-1]) / 2
return result
def simpsons_rule(func, a, b, n):
h = (b - a) / n
x = np.linspace(a, b, n + 1)
y = func(x)
result = h * (y[0] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]) + y[-1]) / 3
return result
def exact_integral(func, a, b):
return quad(func, a, b)[0]
UNIVERSITY OF ENGINEERING AND TECHNOLOGY, TAXILA
FACULTY OF TELECOMMUNICATION AND INFORMATION ENGINEERING
COMPUTER ENGINEERING DEPARTMENT
def calculate_error(func, a, b, n, exact_integral_value):
trapezoidal_result = trapezoidal_rule(func, a, b, n)
simpsons_result = simpsons_rule(func, a, b, n)
errors = {
'Trapezoidal': np.abs(exact_integral_value - trapezoidal_result),
'Simpson\'s 1/3': np.abs(exact_integral_value - simpsons_result)
}
return errors
def plot_error_graph(errors):
methods = list(errors.keys())
values = list(errors.values())
plt.plot(methods, values, marker='o', linestyle='-', color='blue')
plt.xlabel('Integration Methods')
plt.ylabel('Error')
plt.title('Error Comparison between Trapezoidal and Simpson\'s 1/3 Methods')
plt.show()
if __name__ == "__main__":
option = int(input("Choose an option (1 for Linear Equations, 2 for Numerical
Integration): "))
if option == 1:
# Linear Equations
A = np.array([[4, -1, 0, 0], [-1, 4, -1, 0], [0, -1, 4, -1], [0, 0, -1, 3]])
b = np.array([15, 10, 10, 10])
sor_solution, sor_iterations = sor_solver(A, b)
gs_solution, gs_iterations = gauss_seidel_solver(A, b)
print("\nSOR Method:")
print("Solution:", sor_solution)
print("Number of Iterations:", sor_iterations)
print("\nGauss-Seidel Method:")
print("Solution:", gs_solution)
print("Number of Iterations:", gs_iterations)
elif option == 2:
# Numerical Integration
a = 0
b = 2
n = 100 # Number of subintervals
func = lambda x: x ** 2
trapezoidal_result = trapezoidal_rule(func, a, b, n)
simpsons_result = simpsons_rule(func, a, b, n)
exact_integral_value = exact_integral(func, a, b)
print("\nTrapezoidal Method:")
UNIVERSITY OF ENGINEERING AND TECHNOLOGY, TAXILA
FACULTY OF TELECOMMUNICATION AND INFORMATION ENGINEERING
COMPUTER ENGINEERING DEPARTMENT
print("Result:", trapezoidal_result)
print("\nSimpson's 1/3 Method:")
print("Result:", simpsons_result)
errors = calculate_error(func, a, b, n, exact_integral_value)
plot_error_graph(errors)
print("\nExact Integral Value:", exact_integral_value)
else:
print("Invalid option. Please choose 1 or 2.")
OUTPUT:
UNIVERSITY OF ENGINEERING AND TECHNOLOGY, TAXILA
FACULTY OF TELECOMMUNICATION AND INFORMATION ENGINEERING
COMPUTER ENGINEERING DEPARTMENT
