VISVESVARAYA TECHNOLOGICAL UNIVERSITY

“JNANA SANGAMA”, BELAGAVI - 590 018

ASSIGNMENT REPORT on

“Numerical Methods and Applications”

Submitted by

B K KIRAN KUMAR 4SF21CY015

In partial fulfillment of the requirements for the VII semester

BACHELOR OF ENGINEERING

in

COMPUTER SCIENCE & ENGINEERING

Under the Guidance of

Mr. Sunil Kumar K

Assistant Professor, Department of Civil Engineering

at

Sahyadri College of Engineering and Management,

An Autonomous Institution, Mangaluru 2024-25

 MODULE 1
1. Newton-Raphson Method
 Algorithm:
Step 1: Start with an initial guess 𝑥0 , a tolerance value ε and number of
interations.
Step 2: Compute the function value 𝑓(𝑥0) and its derivative 𝑓′(𝑥0)
𝑓(𝑥𝑛)
Step 3: Calculate the next approximation 𝑥𝑛+1 = 𝑥𝑛 −
𝑓′(𝑥𝑛)

Step 4: Check if the difference between 𝑥𝑛 and 𝑥𝑛+1 is less than ∈. If true,
declare 𝑥𝑛+1 as the root and stop .Or if iteration count is greater that
number of iteration, if true stop. Else continue update 𝑥𝑛=𝑥𝑛+1 and
continue the iteration process from step 2.
 Flowchart:
  Reference Materials:

https://www.codesansar.com/numerical-methods/newton-raphson-method-matlab-program-
output.htm

 Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
Sem & Dept: 7CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(f, df, x0, tol=1e-4, max_iter=100):

"""
Newton-Raphson method for solving f(x) = 0.

Parameters:
f : Function for which the root is to be found.
df : Derivative of the function f(x).
x0 : Initial guess.
tol : Tolerance for stopping criterion (default: 1e-4).
max_iter : Maximum number of iterations (default: 100).

Returns:
root : The approximate root.
"""
x = x0
x_values = [x0]
iterations = [0]

for i in range(max_iter):
fx = f(x)
dfx = df(x)

if dfx == 0:
raise ValueError("Derivative is zero. Newton-Raphson method fails.")

x_new = round(x - (fx / dfx), 4)

print(f"Iteration {i+1}: x{i} = {x} and x{i+1} = {x_new}")

# Store values for visualization

x_values.append(x_new)
iterations.append(i+1)

# Check convergence
if abs(x_new - x) < tol:
print(f"Converged in {i+1} iterations.")
return x_new, x_values, iterations

x = x_new

print("Maximum iterations reached without convergence.")

return x, x_values, iterations
def visualize_newton_raphson(f, df, x0, max_iter=100, tolerance=1e-4):
"""
Visualize Newton-Raphson method iterations
"""
# Perform Newton-Raphson method

root, x_values, iterations = newton_raphson(f, df, x0, tolerance, max_iter)

# Create figure with two subplots

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.plot(iterations, x_values, marker='o')
plt.title('Newton-Raphson: Iterations vs X')
plt.xlabel('Iterations')
plt.ylabel('X Value')
plt.grid(True)

# Plot 2: Function values

x_range = np.linspace(min(x_values)-1, max(x_values)+1, 200)
y_values = [f(x) for x in x_range]

plt.subplot(1, 2, 2)
plt.plot(x_range, y_values, label='f(x)')
plt.plot(x_values, [f(x) for x in x_values], 'ro', label='Iterations')
plt.title('Function and Iterations')
plt.xlabel('X')
plt.ylabel('f(x)')
plt.axhline(y=0, color='r', linestyle='--')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

return root

def main():
# Example with specific function
def f(x):
return x + 0.24674*x - 2467.4011

def df(x):
return 1.24674

# Get user input

while True:
try:
x0 = float(input("Enter initial guess: "))
max_iter = int(input("Enter maximum iterations (default 100): ") or 100)
tolerance=float(input("Enter convergence tol(default 1e-4):")or 1e-4)

# Visualize the method

root = visualize_newton_raphson(f, df, x0, max_iter, tolerance)
print(f"\nApproximate Root: {root}")
break
 except ValueError:
print("Invalid input. Please enter numeric values.")

if __name__ == "__main__":
main()

 Output:
2. Matrix Inversion by Gauss-Siedel Method:

 Algorithm:

Step 1: Input Requirements

 Input the coefficient matrix A (square matrix).
 Input the constant vector b.
 Input the initial guess vector x0.
 Input the maximum number of iterations.
 Input the convergence tolerance.
Step 2: Initialization
 Set the iteration counter k=0.
 Create an empty list to store error values.
Step 3: Iterative Solution
 WHILE k<maximum iterations:
a. For each equation (row i):
o Update x[i] using the Gauss-Seidel formula:

o Use the updated x[i] in subsequent calculations.

b. Calculate the relative error.
c. Store the error value.
d. Check for convergence:
o IF error<tolerance, break the loop.
Step 4: Convergence Criteria
 Relative error is calculated as the maximum absolute difference between old and new
values, normalized by the current solution values.
Step 5: Visualization
 Plot error values against iterations to demonstrate convergence behavior.
  Flowchart:

 Reference Material:

https://www.geeksforgeeks.org/gauss-seidel-method/
https://in.mathworks.com/matlabcentral/answers/86085-gauss-seidel-method-in-matlab

 Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
import numpy as np
import matplotlib.pyplot as plt

def get_matrix_input(matrix_type):
"""
Get matrix input from user interactively
"""
while True:
try:
 # Get matrix dimensions
rows = int(input(f"Enter number of rows for {matrix_type} matrix: "))
cols = int(input(f"Enter number of columns for {matrix_type} matrix: "))

# Create matrix to store inputs

matrix = np.zeros((rows, cols), dtype=float)

# Get matrix elements

print(f"\nEnter elements for {matrix_type} matrix:")
for i in range(rows):
for j in range(cols):

matrix[i, j] = float(input(f"Enter element at position [{i},{j}]: "))

return matrix
except ValueError:
print("Invalid input. Please enter numeric values.")

def gauss_seidel(A, b, x0=None, max_iter=100, tolerance=1e-6):

"""
Solve linear system Ax = b using Gauss-Seidel method with detailed calculations
"""
n = len(A)

# Use zero vector as default initial guess

if x0 is None:
x0 = np.zeros_like(b, dtype=float)

x = x0.copy()
errors = []
iterations_details = []

print("\n--- Gauss-Seidel Iteration Details ---")

for k in range(max_iter):
x_old = x.copy()

# Detailed iteration information

iteration_info = {
'iteration': k + 1,
'previous_solution': x_old.copy(),
'updates': []
}

# Gauss-Seidel update
for i in range(n):
# Calculate sigma (sum of known terms)
sigma_left = np.dot(A[i, :i], x[:i]) # Using already updated values
sigma_right = np.dot(A[i, i+1:], x_old[i+1:]) # Using old values for remaining terms

# Calculate new x[i]

x[i] = (b[i] - sigma_left - sigma_right) / A[i, i]

# Store detailed update information

iteration_info['updates'].append({
'index': i,
'sigma_left': sigma_left,
'sigma_right': sigma_right,
'new_value': x[i]

})
 # Calculate error
error = np.max(np.abs((x - x_old) / x)) if np.any(x != 0) else np.max(np.abs(x - x_old))
errors.append(error)

# Store iteration details

iteration_info['current_solution'] = x.copy()
iteration_info['error'] = error
iterations_details.append(iteration_info)

# Print iteration details

print(f"\nIteration {k+1}:")
print("Current Solution:", x)
print("Error:", error)

# Convergence check
if error < tolerance:
break

return x, errors, iterations_details

def plot_convergence(errors):
"""
Plot convergence of Gauss-Seidel method
"""
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(errors) + 1), errors, marker='o')
plt.title('Gauss-Seidel Method Convergence')
plt.xlabel('Iterations')
plt.ylabel('Relative Error')
plt.yscale('log')
plt.grid(True)
plt.show()

def main():
# Get user inputs
print("Gauss-Seidel Method Solver")
print("-------------------------")

# Get coefficient matrix

print("\n1. Coefficient Matrix (A)")
A = get_matrix_input("Coefficient")

# Ensure square matrix

if A.shape[0] != A.shape[1]:
print("Error: Coefficient matrix must be square!")
return

# Get constant vector

print("\n2. Constant Vector (b)")
b = get_matrix_input("Constant")

# Ensure vector matches matrix dimensions

if b.shape[1] != 1 or b.shape[0] != A.shape[0]:
print("Error: Constant vector dimensions must match coefficient matrix!")
return

# Flatten b to 1D array
b = b.flatten()
 # Get initial guess
print("\n3. Initial Guess Vector (x0)")
x0 = get_matrix_input("Initial Guess")

# Flatten x0 to 1D array
x0 = x0.flatten()

# Get additional parameters

max_iter = int(input("\nEnter maximum iterations (default 100): ") or 100)
tolerance = float(input("Enter convergence tolerance (default 1e-6): ") or 1e-6)

# Solve system
solution, errors, iterations = gauss_seidel(A, b, x0, max_iter, tolerance)

print("\n--- Final Results ---")

print("Solution:", solution)
print("Number of iterations:", len(errors))

# Plot convergence
plot_convergence(errors)

if __name__ == "__main__":
main()
 Output:
 MODULE 2
1 Newton’s Forward Difference:
 Algorithm:
Step 1: Start the process of constructing the forward difference table.
Step 2: Input the given data: 𝓍 values (independent variable), Corresponding
𝓎 values (dependent variable).
Step 3: Form the table structure: Arrange 𝓍 values in the first column and
𝓎 values in the first row of the table.
Step 4: Compute forward differences:
 For each column (starting from the second column):
o Compute the differences between consecutive entries in
the previous column using the formula:

Here, represents the 𝒿 -th forward difference of

𝓎 at index 𝓍.
o Continue until all forward differences for all orders are
calculated.
Step 5: Store results in the forward difference table: Populate the 𝑖 -th row and
𝑗 -th column with the computed forward differences.
Step 6: Calculate the parameter 𝑢,using the formula:
𝑎−𝑥[0]
𝑢= where ℎ = 𝑥[1] − 𝑥[0]
ℎ
Step 7: Calculate the value at the interpolation using the formula,
  Flowchart:

 Code:
"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""

import numpy as np
import matplotlib.pyplot as plt
from typing import List, Tuple
class NewtonForwardInterpolation:
def __init__(self, x_data: List[float], y_data: List[float]):
"""
Initialize Newton's Forward Interpolation

Parameters:
x_data: List of x-coordinates
y_data: List of corresponding y-coordinates
"""
self.x_data = np.array(x_data)
self.y_data = np.array(y_data)
self.n = len(x_data)
self.step_size = self.x_data[1] - self.x_data[0]

# Validate inputs
if len(set(np.diff(self.x_data))) > 1:
raise ValueError("Data points must have uniform step size")

def forward_difference_table(self) -> np.ndarray:

"""
Compute forward difference table

Returns:
2D numpy array representing forward differences
"""
# Initialize difference table
diff_table = np.zeros((self.n, self.n))
diff_table[:, 0] = self.y_data

# Compute forward differences

for j in range(1, self.n):
for i in range(self.n - j):
diff_table[i, j] = diff_table[i+1, j-1] - diff_table[i, j-1]

return diff_table

def interpolate(self, x: float) -> float:

"""
Perform Newton's Forward Interpolation

Parameters:
x: Point to interpolate

Returns:
Interpolated y value
"""
# Compute u
u = (x - self.x_data[0]) / self.step_size

# Compute forward difference table

diff_table = self.forward_difference_table()

# Interpolation polynomial
result = diff_table[0, 0]
factorial = 1
 for i in range(1, self.n):
factorial *= i
result += (u * diff_table[0, i]) / factorial
u *= (u - i)

return result

def plot_interpolation(self, interpolation_point: float):

"""
Visualize interpolation results

Parameters:
interpolation_point: x-value to interpolate
"""
plt.figure(figsize=(12, 6))

# Original data points

plt.scatter(self.x_data, self.y_data, color='red', label='Original Data')

# Interpolation point
interp_y = self.interpolate(interpolation_point)
plt.scatter(interpolation_point, interp_y, color='green',
s=200, marker='*', label='Interpolation Point')

# Generate smooth curve for visualization

x_smooth = np.linspace(self.x_data[0], self.x_data[-1], 200)
y_smooth = [self.interpolate(x) for x in x_smooth]
plt.plot(x_smooth, y_smooth, label='Interpolation Polynomial')

plt.title("Newton's Forward Interpolation")

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

# Print interpolation details

print(f"Interpolation Point: {interpolation_point}")
print(f"Interpolated Value: {interp_y}")

def get_user_inputs():
"""
Collect user inputs interactively
"""
# Data points input
print("Enter data points (x, y)")
x_data, y_data = [], []
while True:
try:
x = float(input("Enter x (or press Enter to finish): "))
y = float(input("Enter corresponding y: "))
x_data.append(x)
y_data.append(y)
except ValueError:
break
 # Interpolation point
interpolation_point = float(input("Enter point to interpolate: "))

return x_data, y_data, interpolation_point

def main():
# Get user inputs
x_data, y_data, interpolation_point = get_user_inputs()

# Create interpolation object

interpolator = NewtonForwardInterpolation(x_data, y_data)

# Compute and display forward difference table

diff_table = interpolator.forward_difference_table()
print("\nForward Difference Table:")
print(diff_table)

# Interpolate and plot

interpolator.plot_interpolation(interpolation_point)

if __name__ == "__main__":
main()

 Output:
2 Lagrange’s Interpolation Formula:
 Algorithm:
Step 1:Given a set of 𝑛 data points (𝑥0, 𝑦0),( 𝑥1, 𝑦1),...,( 𝑥𝑛−1, 𝑦𝑛−1) where 𝑥𝑖
are the known 𝑥-coordinates and 𝑦𝑖 are the corresponding 𝑦-coordinates.
Given the point 𝑎 (the value of 𝑥where you want to estimate 𝑦).
Step 2: Set Up Lagrange Polynomial: The Lagrange interpolation polynomial is
given by:
−1
𝑓(𝑎) = ∑𝑛𝑖=0 𝑦𝑖𝐿𝑖(𝑎)
Where 𝐿𝑖(𝑎) is given as

Or Lagrange Polynomial can also be given as:

Step 3: The result will give the interpolated value 𝑓(𝑎) at 𝑎.

 Flowchart:
  Code:
"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""

import numpy as np
import matplotlib.pyplot as plt

def lagrange_interpolation(x_points, y_points, x):

"""
Perform Lagrange Interpolation

Parameters:
x_points: list of x coordinates of known points
y_points: list of y coordinates of known points
x: x coordinate for interpolation

Returns:
Interpolated y value
"""
n = len(x_points)
y=0

for i in range(n):
# Calculate Lagrange basis polynomial
p=1
for j in range(n):
if i != j:
p *= (x - x_points[j]) / (x_points[i] - x_points[j])
y += p * y_points[i]

return y

def plot_lagrange_interpolation(x_points, y_points):

"""
Plot Lagrange Interpolation

Parameters:
x_points: list of x coordinates of known points
y_points: list of y coordinates of known points
"""
# Create dense interpolation points
x_interp = np.linspace(min(x_points), max(x_points), 200)
y_interp = [lagrange_interpolation(x_points, y_points, x) for x in x_interp]

# Plot original points and interpolation curve

plt.figure(figsize=(10, 6))
plt.plot(x_interp, y_interp, 'b-', label='Interpolated Curve')
plt.scatter(x_points, y_points, color='red', label='Original Points')

plt.title('Lagrange Interpolation')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.grid(True)
plt.show()
def main():
# Input data points
print("Lagrange Interpolation")
print("Enter the number of points:")
n = int(input())

x_points = []
y_points = []

print("Enter x and y coordinates:")

for i in range(n):
x = float(input(f"Enter x{i+1}: "))
y = float(input(f"Enter y{i+1}: "))
x_points.append(x)
y_points.append(y)

# Plot the interpolation

plot_lagrange_interpolation(x_points, y_points)

if __name__ == "__main__":
main()
 Output:
 MODULE 3
1. Trapezoid Method for Numerical Integration:
 Algorithm:
Step 1: Start the process.
Step 2: Input the required data.
 Enter the lower limit (𝑎) of integration
 Enter the upper limit (𝑏) of integration.
 Enter the number of subintervals (𝑛)
Step 3: Define the function 𝑓(𝑥)
Step 4: Calculate the step size (ℎ)
 𝑏−𝑎
ℎ=
𝑛

Step 5: Compute the initial sum of the boundary values.

 Initialize the result as: 𝑟𝑒𝑠𝑢𝑙𝑡 = 𝑓(𝑎) + 𝑓(𝑏)
Step 6: Iterate through intermediate points.
 Loop from 𝑖 = 1 to 𝑖 = 𝑛 − 1
o Calculate the intermediate 𝑥𝑖 as: 𝑥𝑖 = 𝑎 + 𝑖. ℎ
o Update the result by adding 2. 𝑓(𝑥𝑖): 𝑟𝑒𝑠𝑢𝑙𝑡+= 2𝑓(𝑥𝑖)
ℎ
Step 7: Multiply the result by the factor :
2

 Integration value =ℎ . 𝑟𝑒𝑠𝑢𝑙𝑡

2
  Flowchart:

 Code:

import numpy as np
import matplotlib.pyplot as plt
from typing import Callable, Union

class TrapezoidalIntegration:
def __init__(self, func: Callable[[float], float]):
 """
Initialize Trapezoidal Integration

Parameters:
func: Integrand function f(x)
"""
self.func = func

def integrate(self,
a: float,
b: float,
n: int = 100) -> dict:
"""
Compute definite integral using Trapezoidal Rule

Parameters:
a: Lower integration limit
b: Upper integration limit
n: Number of subintervals

Returns:
Dictionary with integration results
"""
# Validate inputs
if a > b:
a, b = b, a

# Compute step size

h = (b - a) / n

# Initialize sum
integral_sum = 0.5 * (self.func(a) + self.func(b))

# Compute intermediate points

for i in range(1, n):
x=a+i*h
integral_sum += self.func(x)

# Final integration
integral_value = h * integral_sum

# Error estimation (Simpson's Rule as comparison)

simpson_integral = self._simpson_integration(a, b, n)
absolute_error = abs(simpson_integral - integral_value)
relative_error = absolute_error / abs(simpson_integral) * 100

return {
'integral': integral_value,
'method': 'Trapezoidal Rule',
'subintervals': n,
'absolute_error': absolute_error,
'relative_error': relative_error
}
def _simpson_integration(self,
a: float,
b: float,
n: int) -> float:
"""
Simpson's Rule for error comparison
 Parameters:
a: Lower integration limit
b: Upper integration limit
n: Number of subintervals

Returns:
Integral approximation
"""
h = (b - a) / n

def simpson_sum():
integral_sum = self.func(a) + self.func(b)

# Even terms
for i in range(1, n, 2):
x=a+i*h
integral_sum += 4 * self.func(x)

# Odd terms
for i in range(2, n-1, 2):
x=a+i*h
integral_sum += 2 * self.func(x)

return integral_sum * h / 3

return simpson_sum()

def visualize_integration(self,
a: float,
b: float,
n: int = 100):
"""
Visualize Trapezoidal Integration

Parameters:
a: Lower integration limit
b: Upper integration limit
n: Number of subintervals
"""
# Generate smooth function points
x = np.linspace(a, b, 200)
y = [self.func(xi) for xi in x]

# Compute trapezoidal points

x_trapezoid = np.linspace(a, b, n+1)

y_trapezoid = [self.func(xi) for xi in x_trapezoid]

# Plotting
plt.figure(figsize=(12, 6))
plt.plot(x, y, label='Function f(x)', color='blue')
plt.fill_between(x_trapezoid, y_trapezoid, alpha=0.3,
color='green', label='Trapezoidal Regions')

plt.title('Trapezoidal Rule Integration')

plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
 plt.show()

def get_user_inputs():
"""
Collect user inputs interactively
"""
# Function selection
print("Select Integration Function:")
print("1. f(x) = x^2")
print("2. f(x) = sin(x)")
print("3. f(x) = exp(x)")
print("4. Custom function")

choice = int(input("Enter your choice (1-4): "))

# Function mapping
functions = {
1: lambda x: x**2,
2: np.sin,
3: np.exp,
}

if choice == 4:
# Custom function input
def custom_func(x):
return eval(input("Enter function in terms of x: "))
func = custom_func
else:
func = functions[choice]

# Integration limits
a = float(input("Enter lower limit (a): "))
b = float(input("Enter upper limit (b): "))

# Number of subintervals
n = int(input("Enter number of subintervals: "))

return func, a, b, n

def main():
# Get user inputs
func, a, b, n = get_user_inputs()

# Create integration object

integrator = TrapezoidalIntegration(func)

# Perform integration
result = integrator.integrate(a, b, n)

# Print results
print("\n--- Integration Results ---")
for key, value in result.items():
print(f"{key.replace('_', ' ').title()}: {value}")

# Visualize integration
integrator.visualize_integration(a, b, n)
if __name__ == "__main__":
"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
main()

 Output:
2. Simpsons 𝟏Ú𝟑 𝒓𝒅 rule for Numerical Integration:
 Algorithm:
Step 1: Start the process.
Step 2: Input the required data.
 Enter the lower limit (𝑎) of integration
 Enter the upper limit (𝑏) of integration.
 Enter the number of subintervals (𝑛)
Step 3: Define the function 𝑓(𝑥)
Step 4: Calculate the step size (ℎ)
 𝑏−𝑎
ℎ= 𝑛
Step 5: Compute the initial sum of the boundary values.
 Initialize the result as: 𝑟𝑒𝑠𝑢𝑙𝑡 = 𝑓(𝑎) + 𝑓(𝑏)
Step 6: Sum Odd Terms:
 Loop through all odd indices from 1 to 𝑛 − 1:
o Compute 𝑥𝑖 = 𝑎 + 𝑖. ℎ
o Add 4.𝑓(𝑥𝑖) to the 𝑟𝑒𝑠𝑢𝑙𝑡
Step 7: Sum Even Terms:
 Loop through all odd indices from 2 to 𝑛 − 2:
o Compute 𝑥𝑖 = 𝑎 + 𝑖. ℎ
o Add 2.𝑓(𝑥𝑖) to the 𝑟𝑒𝑠𝑢𝑙𝑡
ℎ
Step 8: Multiply the result by the factor :
3

 Integration value =ℎ . 𝑟𝑒𝑠𝑢𝑙𝑡

3
 Flowchart:

 Code:
"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
import numpy as np
import matplotlib.pyplot as plt
import math

def simpsons_rule(func, a, b, n):

 """
Perform Simpson's 1/3 Rule Integration

Parameters:
func: Integration function
a: Lower bound

b: Upper bound
n: Number of subintervals (must be even)

Returns:
Integral approximation
"""
# Ensure even number of subintervals
if n % 2 != 0:
n += 1

h = (b - a) / n
x = np.linspace(a, b, n+1)
y = [func(xi) for xi in x]

integral = h/3 * (y[0] + y[-1] +

4*sum(y[1:-1:2]) +
2*sum(y[2:-1:2]))

return integral, x, y

def plot_integration(func, a, b, n):

"""
Plot function and integration area
"""
integral, x, y = simpsons_rule(func, a, b, n)

# Generate dense points for smooth curve

x_smooth = np.linspace(a, b, 200)
y_smooth = [func(xi) for xi in x_smooth]

plt.figure(figsize=(10, 6))
plt.plot(x_smooth, y_smooth, 'b-', label='Function')
plt.fill_between(x, 0, y, alpha=0.3, color='red', label='Integral Area')

plt.title(f"Simpson's 1/3 Rule Integration\nIntegral ≈ {integral:.4f}")

plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.legend()
plt.show()

return integral

def main():
print("Simpson's 1/3 Rule Integration")

# Get function from user

print("Enter function (use math functions, x as variable):")
func_str = input()
func = lambda x: eval(func_str)
 # Get integration bounds
print("Enter lower bound:")

a = float(input())
print("Enter upper bound:")
b = float(input())

# Get number of subintervals

print("Enter number of subintervals (even number):")
n = int(input())

# Perform integration and plot

result = plot_integration(func, a, b, n)
print(f"Integral approximation: {result:.4f}")

if __name__ == "__main__":
main()

 Output:
 MODULE 4
1 Euler Series Method:
 Algorithm:

Step 1: Input Requirements

Define differential equation function f(x,y).

Set initial values: x0, y0, step size h, and endpoint xend.

Step 2: Initialization

Set x0, y0, and create lists to store x-values and y-values.

Step 3: Iterative Process

While x0<xend:
Update y0=y0+h⋅f(x0,y0)
Update x0=x0+h
Append x0 and y0 to the lists.

Step 4: Plotting

Plot x-values vs. y-values.

Step 5: Output

Print the final approximation of y at x=xend.

Step 6: End
 Flowchart:
  Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""

import matplotlib.pyplot as plt

def get_differential_equation():
"""
Get the differential equation from the user.
Returns:
callable: A lambda function representing dy/dx.
"""
print("Enter the differential equation dy/dx as a function of x and y (e.g., x + y + x * y):")
eq_str = input()
return lambda x, y: eval(eq_str)

def euler(func, x0, y0, h, x_end):

"""
Solve the differential equation using Euler's method.

Parameters:
func (callable): Differential equation dy/dx as a function of x and y.
x0 (float): Initial x value.
y0 (float): Initial y value at x0.
h (float): Step size.
x_end (float): Target x value.

Returns:
None
"""
x_values = [x0]
y_values = [y0]

# Iterating until reaching the specified x

while x0 < x_end:
y0 = y0 + h * func(x0, y0)
x0 = x0 + h
x_values.append(x0)
y_values.append(y0)

# Plotting the solution

plt.plot(x_values, y_values, label="Euler's Method", marker='o')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Numerical Solution using Euler's Method")
plt.legend()
plt.grid()
plt.show()
 # Printing the final approximation
print(f"Approximate solution at x = {x_end} is {y0:.6f}")

def main():
"""
Main function to get user inputs and solve the differential equation.
"""
print("Euler's Method Numerical Solver")

# Get the differential equation from the user

func = get_differential_equation()

# Get initial conditions and target x from the user

x0 = float(input("Enter the initial value of x (x0): "))
y0 = float(input("Enter the initial value of y (y0): "))
h = float(input("Enter the step size (h): "))
x_end = float(input("Enter the target value of x (x_end): "))

# Solve the differential equation

euler(func, x0, y0, h, x_end)

if __name__ == "__main__":
main()

 Output:
2 Runge Kutta 4th Order Method:
 Algorithm:
Step 1: Start
Step 2: Define function f(x,y)
Step 3: Read values of initial condition (x0,y0) number of steps nnn, and calculation
point xn
Step 4: Calculate step size h=(xn−x0)/n
Step 5: Set i=0
Step 6: Loop
 Compute k1=h⋅f(x0,y0)
 Compute k2=h⋅f(x0+h2,y0+k12)
 Compute k3=h⋅f(x0+h/2,y0+k2/2)
 Compute k4=h⋅f(x0+h,y0+k3)
 Compute k=(k1+2k2+2k3+k4)/6
 Compute yn=y0+k
 Increment i=i+1
 Update x0=x0+h
 Update y0=yn
 Repeat while i<ni < ni<n
Step 7: Display yny_nyn as result
Step 8: Stop
 Flowchart:

 Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
def get_differential_equation():
"""
Get the differential equation from the user.
Returns:
callable: A lambda function representing dy/dx.
"""
print("Enter the differential equation dy/dx as a function of x and y (e.g., (x - y)/2):")
eq_str = input()
return lambda x, y: eval(eq_str)

def runge_kutta(dydx, x0, y0, x, h):

"""
Solve the differential equation using the 4th order Runge-Kutta method.

Parameters:
dydx (callable): Differential equation dy/dx as a function of x and y.
x0 (float): Initial x value.
y0 (float): Initial y value at x0.
x (float): Target x value.
h (float): Step size.

Returns:
float: Approximate value of y at x.
"""
# Calculate the number of steps
n = int((x - x0) / h)
y = y0

# Perform Runge-Kutta iteration

for _ in range(n):
k1 = h * dydx(x0, y)
k2 = h * dydx(x0 + 0.5 * h, y + 0.5 * k1)
k3 = h * dydx(x0 + 0.5 * h, y + 0.5 * k2)
k4 = h * dydx(x0 + h, y + k3)

# Update y and x for the next iteration

y += (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
x0 += h

return y

def main():
"""
Main function to get inputs from the user and solve the differential equation.
"""
print("Runge-Kutta 4th Order Method Solver")

# Get the differential equation from the user

dydx = get_differential_equation()

# Get initial conditions and target x from the user

x0 = float(input("Enter the initial value of x (x0): "))
y0 = float(input("Enter the initial value of y (y0): "))
x = float(input("Enter the target value of x: "))
h = float(input("Enter the step size (h): "))
 # Solve the differential equation
result = runge_kutta(dydx, x0, y0, x, h)

# Display the result0

print(f"\nThe value of y at x = {x} is approximately: {result:.6f}")

if __name__ == "__main__":
main()

 Output:
 MODULE 5
1 Laplace Equation Method:
 Algorithm:
Step 1: Input Requirements
 Set boundary conditions (top, bottom, left, right).
 Initialize interior points (U1 to U9) to 0.
 Set tolerance (1e-4) and max iterations (1000).
Step 2: Initialize Grid
 Create a 3x3 grid for U.
Step 3: Iterative Process
 For each iteration:
o Copy current grid to U_old.
o Update interior points (U1 to U9) using Laplace's equation.
Step 4: Convergence Check
 If the maximum change between the old and new grid is less than tolerance, stop.
Step 5: Termination
 If max iterations are reached without convergence, stop.
Step 6: Output
 Print the final grid values.

 Flowchart:
 Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
import numpy as np

# Initialize the grid

# Boundary conditions (edges of the grid)
boundary_top = [1000, 1000, 1000]
boundary_bottom = [2000, 2000, 2000]
boundary_left = [500, 1000, 2000]
boundary_right = [1000, 500, 1500]

# Initialize interior points with zeros or guesses

U = np.zeros((3, 3)) # For U1 to U9 (interior points)

# Tolerance for convergence

tolerance = 1e-4
max_iterations = 1000

# Iterative process
for iteration in range(max_iterations):
U_old = U.copy()

# Update each interior point using Laplace's equation

U[0, 0] = 0.25 * (boundary_top[0] + boundary_left[0] + U[0, 1] + U[1, 0]) # U1
U[0, 1] = 0.25 * (boundary_top[1] + U[0, 0] + U[0, 2] + U[1, 1]) # U2
U[0, 2] = 0.25 * (boundary_top[2] + U[0, 1] + U[1, 2] + boundary_right[0])# U3
U[1, 0] = 0.25 * (U[0, 0] + boundary_left[1] + U[2, 0] + U[1, 1]) # U4
U[1, 1] = 0.25 * (U[0, 1] + U[1, 0] + U[1, 2] + U[2, 1]) # U5
U[1, 2] = 0.25 * (U[0, 2] + U[1, 1] + boundary_right[1] + U[2, 2]) # U6
U[2, 0] = 0.25 * (U[1, 0] + boundary_left[2] + boundary_bottom[0] + U[2, 1]) # U7
U[2, 1] = 0.25 * (U[2, 0] + U[1, 1] + boundary_bottom[1] + U[2, 2]) # U8
U[2, 2] = 0.25 * (U[2, 1] + U[1, 2] + boundary_bottom[2] + boundary_right[2]) # U9

# Check for convergence

if np.max(np.abs(U - U_old)) < tolerance:
print(f"Converged in {iteration + 1} iterations.")
break
else:
print("Maximum iterations reached without convergence.")

# Display the final grid

print("Final values of interior nodes:")
print(U)
  Output:

2 Poisson’s Equation Solver :

 Algorithm:
Step 1: Input Requirements

Get grid size (nx, ny) and source function f(x,y) from the user.

Step 2: Initialize Grid

Calculate grid spacings dx and dy, and create meshgrid X,Y.

Step 3: Setup Source Term

Initialize vector bb and populate it with the evaluated source term for interior points.

Step 4: Construct Sparse Matrix (A)

Build matrix A using main and off-diagonals based on dx and dy.

Step 5: Apply Boundary Conditions

Set boundary values in vector bb to zero.

Step 6: Solve Linear System

Solve A⋅solution=bA⋅solution=b using spsolve.

Step 7: Reshape Solution

Reshape the solution vector into a 2D grid.

Step 8: Plot Solution

Plot the solution using contour and surface plots.

Step 9: End

Display the solution and plots.

 Flowchart:

 Code:

"""
Name: B K Kiran Kumar
USN: 4SF21CY015
DEPT: 7 CY-CSE
Faculty Incharge: Mr. Sunil Kumar K (Assistant Professor)
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

"""
Poisson Equation Solver in 2D

This script solves the 2D Poisson equation on a rectangular grid using

the finite difference method. The solution is visualized with contour
and surface plots.

Functions:
1. poisson_solver_2d(nx, ny, source_func):
Solves the 2D Poisson equation using a sparse matrix representation
of the discretized grid.

Parameters:
nx : Number of grid points in the x-direction.
ny : Number of grid points in the y-direction.
 source_func : Function defining the source term f(x, y).

Returns:
X, Y : Mesh grid for the solution domain.
solution_grid : Solution of the Poisson equation as a 2D array.

2. plot_poisson_solution(X, Y, solution):
Visualizes the solution using contour and surface plots.

3. main():
Accepts user input for grid size and source function, solves the
Poisson equation, and visualizes the solution.

"""

def poisson_solver_2d(nx, ny, source_func):

dx = 1 / (nx - 1)
dy = 1 / (ny - 1)

x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)

b = np.zeros((nx*ny))
for i in range(1, nx-1):
for j in range(1, ny-1):
idx = i + j*nx
b[idx] = source_func(X[j,i], Y[j,i]) * (dx*dy)

main_diag = -2*(1/dx**2 + 1/dy**2) * np.ones(nx*ny)

off_diag_x = (1/dx**2) * np.ones(nx*ny - 1)
off_diag_y = (1/dy**2) * np.ones(nx*ny - nx)
diagonals = [main_diag, off_diag_x, off_diag_x, off_diag_y, off_diag_y]
offsets = [0, -1, 1, -nx, nx]

A = diags(diagonals, offsets, shape=(nx*ny, nx*ny))

b[0:nx] = 0
b[-nx:] = 0
b[::nx] = 0
b[nx-1::nx] = 0

solution = spsolve(A, b)
solution_grid = solution.reshape((ny, nx))

return X, Y, solution_grid

def plot_poisson_solution(X, Y, solution):

plt.figure(figsize=(10, 8))
plt.subplot(2, 1, 1)
plt.contourf(X, Y, solution, cmap='viridis')
plt.colorbar()
plt.title("Poisson Equation Solution - Contour")
plt.subplot(2, 1, 2)
plt.pcolormesh(X, Y, solution, cmap='viridis')
 plt.colorbar()
plt.title("Poisson Equation Solution - Surface")
plt.tight_layout()
plt.show()

def main():
print("Enter grid size (nx ny):")
nx, ny = map(int, input().split())

print("Enter source term function f(x,y):")

source_str = input()
source_func = lambda x, y: eval(source_str)

X, Y, solution = poisson_solver_2d(nx, ny, source_func)

plot_poisson_solution(X, Y, solution)

if __name__ == "__main__":
main()

 Code 2:

import numpy as np

# Function to calculate the source term based on x, y

def source_function(x, y):
return -10 * (x*2 + y*2 + 10)

# Input the grid size (N x N grid)

grid_size = int(input("Enter the grid size (e.g., 4 for a 4x4 grid): "))
h = int(input("Enter the step size (h): "))

# Take user inputs for boundary values

top_boundary = list(map(float, input(f"Enter {grid_size} top boundary values (space-separated): ").split()))
bottom_boundary = list(map(float, input(f"Enter {grid_size} bottom boundary values (space-separated):
").split()))
left_boundary = list(map(float, input(f"Enter {grid_size} left boundary values (space-separated):
").split()))
right_boundary = list(map(float, input(f"Enter {grid_size} right boundary values (space-separated):
").split()))

# Generate interior points (ignoring boundary points)

interior_points = []
boundary_value = 0 # Default Dirichlet condition for unused edges (if any)

for i in range(1, grid_size - 1): # Rows (excluding boundary rows)

for j in range(1, grid_size - 1): # Columns (excluding boundary columns)
interior_points.append((i, j)) # Store grid indices as (row, col)

num_points = len(interior_points) # Total number of interior points

# Map indices for interior points

point_indices = {point: idx for idx, point in enumerate(interior_points)}

# Dynamically construct the coefficient matrix A

A = np.zeros((num_points, num_points))
# Compute the source term vector b
b = np.zeros(num_points)

# Populate the coefficient matrix A and vector b

for point, idx in point_indices.items():
i, j = point
A[idx, idx] = -4 # Diagonal (u_ij coefficient)

# Left neighbor
if (i, j - 1) in point_indices:

A[idx, point_indices[(i, j - 1)]] = 1

else:
b[idx] -= left_boundary[i] # Contribution from left boundary

# Right neighbor
if (i, j + 1) in point_indices:
A[idx, point_indices[(i, j + 1)]] = 1

else:
b[idx] -= right_boundary[i] # Contribution from right boundary

# Top neighbor
if (i - 1, j) in point_indices:
A[idx, point_indices[(i - 1, j)]] = 1
else:
b[idx] -= top_boundary[j] # Contribution from top boundary

# Bottom neighbor
if (i + 1, j) in point_indices:
A[idx, point_indices[(i + 1, j)]] = 1
else:
b[idx] -= bottom_boundary[j] # Contribution from bottom boundary

# Compute source term

x, y = i * h, j * h # Convert grid indices to coordinates
b[idx] += source_function(x, y)

# Solve the linear system A * U = b

interior_values = np.linalg.solve(A, b)

# Map results back to the grid

grid = np.zeros((grid_size, grid_size)) # Initialize with zeros
grid[0, :] = top_boundary
grid[-1, :] = bottom_boundary
grid[:, 0] = left_boundary
grid[:, -1] = right_boundary

for (i, j), idx in point_indices.items():

grid[i, j] = interior_values[idx]

# Display the coefficient matrix A

print("\nCoefficient Matrix A:")
print(A)
 # Display the right-hand side vector b
print("\nSource Term Vector b:")
print(b)

# Display the updated grid with solved values

print("\nUpdated Grid with Solved Interior Values:")
print(grid)

 Output:
Code 2 Output: