User avatar
@anonymous · 4h
1.Gauss elimination
def printMatrix(A):
for row in A:
print("[", end="")
for elem in row:
print(round(elem, 5), " ", end="")
print("]")
def GaussElimination(n, A, b):
for k in range(n - 1):
for i in range(k + 1, n):
factor = A[i][k] / A[k][k]
for j in range(k, n):
A[i][j] = A[i][j] - factor * A[k][j]
b[i] = b[i] - factor * b[k]
print(f"\nRow operation in row: {i+1}")
printMatrix(A)
x = [0] * n
x[n - 1] = b[n - 1] / A[n - 1][n - 1]
for i in range(n - 2, -1, -1):
total = b[i]
for j in range(i + 1, n):
total -= A[i][j] * x[j]
x[i] = total / A[i][i]
return x
N = int(input("Enter the size of the system: "))
print("Enter the coefficient matrix A:")
matrix = []
for i in range(N):
row = list(map(float, input().split()))
matrix.append(row)
print("Enter the RHS vector b:")
b = list(map(float, input().split()))
solution = GaussElimination(N, matrix, b)
print("\nSolution Vector:")
for i in range(N):
print(f"x{i+1} = {solution[i]:.2f}")
OUTPUT
Enter the size of the matrix: 3
Enter the matrix A:
2 1 -1
-3 -1 2
-2 1 2
Enter the answer vector b:
8 -11 -3
2.GAUSS ELI AND ROW PIVOTING
from fractions import Fraction
def printMatrix(A):
for row in A:
print("[",end="")
for elem in row:
print(round(Fraction(elem),5)," ",end="")
print("]",end="")
print()
def gauss_elimination(n, A, b):
TOL = 1e-10
for k in range(n - 1):
if abs(A[k][k]) <= TOL:
print("The first row first column in 0")
max_row = k
max_value = abs(A[k][k])
for p in range(k + 1, n):
if abs(A[p][k]) > max_value:
max_value = abs(A[p][k])
max_row = p
if max_value <= TOL:
raise ValueError("No unique solution exists.")
A[k], A[max_row] = A[max_row], A[k]
b[k], b[max_row] = b[max_row], b[k]
print(f"After {k} swap")
for i in range(k + 1, n):
factor = A[i][k] / A[k][k]
for j in range(k + 1, n):
A[i][j] = A[i][j] - factor * A[k][j]
A[i][k] = 0
b[i] = b[i] - factor * b[k]
print(f" -- Row operation in row: {i+1}")
printMatrix(A)
x = [0] * n
x[n - 1] = b[n - 1] / A[n - 1][n - 1]
for k in range(n - 2, -1, -1):
sum_ax = 0
for j in range(k + 1, n):
sum_ax += A[k][j] * x[j]
x[k] = (b[k] - sum_ax) / A[k][k]
return x
n = int(input("Enter the size of the system: "))
print("Enter the coefficient matrix A:")
A = []
for i in range(n):
row = list(map(float, input().split()))
A.append(row)
print("Enter the RHS vector b:")
b = list(map(float, input().split()))
solution = gauss_elimination(n, A, b)
print("\nSolution Vector:")
for i in range(n):
print(f"x{i+1} = {solution[i]:.2f}")
OUTPUT
Enter the size of the system: 3
Enter the coefficient matrix A:
2 1 -1
-3 -1 2
-2 1 2
Enter the RHS vector b:
8 -11 -3
3.JACOBI
def gauss_jacobi(A,b,x_old,TOL,KMAX):
n=len(A)
for i in range(n):
print(f"Old x[{i}] = {x_old[i]}")
TAG = True
for k in range(1,KMAX+1):
if TAG==False:
break
x_new = [0]*n
for i in range(n):
summ=0
for j in range(n):
if j!=i:
summ+=A[i][j]*x_old[j]
x_new[i]=(b[i]-summ)/A[i][i]
print(f"x_new[{i}] = {x_new[i]}")
error = max(abs(x_new[i]-x_old[i]) for i in range(n))
print(f"Iteration {k}: {x_new}")
print(f"Error: {error}")
if error<=TOL:
TAG=False
else:
for i in range(n):
x_old[i]=x_new[i]
return x_new
n = int(input("Enter no. of equ: "))
A=[]
print("Enter the coeff. matrix: ")
for i in range(n):
row = list(map(float,input(f"Row {i+1}: ").split()))
A.append(row)
b = list(map(float,input("Enter RHS value: ").split()))
x0=list(map(float,input("Enter inital guess values: ").split()))
TOL=float(input("Enter tolerance value: "))
KMAX=int(input("Etner maximum no. of iteration: "))
ans = gauss_jacobi(A,b,x0,TOL,KMAX)
print("\nFinal Answer: ")
for i in range(n):
print(f"x{i+1} = {ans[i]}")
OUTPUT
Enter no. of equ: 3
Enter the coeff. matrix:
Row 1: 10 1 1
Row 2: 2 10 1
Row 3: 2 2 10
Enter RHS value: 12 13 14
Enter inital guess values: 0 0 0
Enter tolerance value: 0.001
Etner maximum no. of iteration: 20
4.SEIDAL
def gauss_seidal(A,b,x_old,TOL,KMAX):
n=len(b)
for k in range(1,KMAX+1):
x_new = x_old.copy()
for i in range(n):
summ=0
for j in range(i):
summ+=A[i][j]*x_new[j]
for j in range(i+1,n):
summ+=A[i][j]* x_old[j]
x_new[i]=(b[i]-summ)/A[i][i]
error = max(abs(x_new[i]-x_old[i]) for i in range(n))
print(f"Iteration {k}: {x_new}")
print(f"Error: {error}")
if error<=TOL:
break
x_old=x_new.copy()
return x_new
n=int(input("Enter no.of equ: "))
A=[]
print("Enter coeff. matrix: ")
for i in range(n):
row = list(map(float,input(f"Row {i+1}: ").split()))
A.append(row)
b = list(map(float,input("Enter RHS value: ").split()))
x0=list(map(float,input("Enter inital guess values: ").split()))
TOL=float(input("Enter tolerance value: "))
KMAX=int(input("Enter maximum no. of iteration: "))
ans = gauss_seidal(A,b,x0,TOL,KMAX)
print("\nFinal Answer: ")
for i in range(n):
print(f"x{i+1} = {ans[i]}")
OUTPUT
Enter no.of equ: 3
Enter coeff. matrix:
Row 1: 10 1 1
Row 2: 2 10 1
Row 3: 2 2 10
Enter RHS value: 12 13 14
Enter inital guess values: 0 0 0
Enter tolerance value: 0.001
Enter maximum no. of iteration: 20
5.POWER METHOD
import math
MAX = 10
def main():
n = int(input("Enter order of matrix: "))
print("Enter matrix:")
A = []
for i in range(n):
row = list(map(float, input().split()))
A.append(row)
print("Enter initial vector:")
V = list(map(float, input().split()))
tol = float(input("Enter tolerance: "))
kmax = int(input("Enter maximum iterations: "))
lambda_val = 0.0
tag = True
W = [0.0] * n
for k in range(1, kmax + 1):
if not tag:
break
for i in range(n):
W[i] = 0.0
for j in range(n):
W[i] = W[i] + A[i][j] * V[j]
c1 = W[0]
for i in range(n):
V[i] = W[i] / c1
lambda_error = lambda_val - c1
if lambda_error < 0:
lambda_error = -lambda_error
total_sum = 0.0
print(f"\nIteration {k}")
print(f"Eigen Value = {lambda_val}")
print("Eigen Vector: ", end=" ")
for i in range(n):
print(f"{V[i]}", end=" ")
print()
if lambda_error <= tol:
tag = False
lambda_val = c1
print(f"\nDominant Eigen Value = {lambda_val}")
print("Dominant Eigen Vector = ", end=" ")
for i in range(n):
print(f"{V[i]}", end=" ")
print()
if __name__ == "__main__":
main()
OUTPUT
Enter order of matrix: 2
Enter matrix:
4 1
2 3
Enter initial vector:
1 1
Enter tolerance: 0.001
Enter maximum iterations: 20
6.JORDAN INVERSE
n = int(input("Enter the size of matrix: "))
A = []
for i in range(n):
row = list(map(float, input().split()))
A.append(row)
aug = []
for i in range(n):
identity = [1 if i == j else 0 for j in range(n)]
aug.append(A[i] + identity)
print("\nInitial Augmented Matrix:")
for row in aug:
print(" ".join(f"{x:.3f}" for x in row))
for i in range(n):
if aug[i][i] == 0:
for k in range(i + 1, n):
if aug[k][i] != 0:
aug[i], aug[k] = aug[k], aug[i]
break
else:
print("Inverse does not exist")
exit()
pivot = aug[i][i]
for j in range(2 * n):
aug[i][j] /= pivot
for k in range(n):
if k != i:
factor = aug[k][i]
for j in range(2 * n):
aug[k][j] -= factor * aug[i][j]
print(f"\nIteration {i+1}:")
for row in aug:
print(" ".join(f"{x:.3f}" for x in row))
print("\nInverse matrix:")
for i in range(n):
print(" ".join(f"{aug[i][j]:.3f}" for j in range(n, 2 * n)))
print("\nA * A-inverse:")
result = []
for i in range(n):
row = []
for j in range(n):
total = 0
for k in range(n):
total += A[i][k] * aug[k][n+j]
row.append(total)
result.append(row)
for row in result:
print(" ".join(f"{x:.3f}" for x in row))
OUTPUT
Enter the size of matrix: 3
1 2 3
0 1 4
5 6 0
7.LAGRANGE
import sympy as sp
n = int(input("Enter the number of points: "))
l = []
for i in range(n):
row = list(map(float, input().split()))
l.append(row)
x_value = float(input("Enter x value: "))
x = sp.Symbol('x')
value_count=0
sum = 0
term2 = 0
for i in range(n):
p = 1
term1 = l[i][1]
for j in range(n):
if j!=i:
term1*=(x-l[j][0])/(l[i][0]-l[j][0])
p*=(x_value-l[j][0])/(l[i][0]-l[j][0])
value_count+=4
sum+=p*l[i][1]
value_count+=2
term2+=term1*l[i][1]
poly_form = sp.expand(term2)
degree = sp.degree(poly_form, x)
print(x, x_value)
print("\nLagrange form : ",term2)
print("Polynomial form: ",poly_form)
print(f"Interpolation of {x_value} :",sum)
print("No of operations to get the value: ",value_count)
print("No of operations to get the polynomial eqn: ",sp.count_ops(poly_form))
print("No of operations to interpolate at x: ",2*degree)
OUTPUT
Enter the number of points: 3
1 1
2 4
3 9
Enter x value: 2.5
8.NEWTON DIVIDED DIFF
n = int(input("Enter number of data points: "))
x = []
y = []
print("Enter x and y values:")
for i in range(n):
xi = float(input(f"x[{i}] = "))
yi = float(input(f"y[{i}] = "))
x.append(xi)
y.append(yi)
d = [[0.0 for j in range(n)] for i in range(n)]
for j in range(n):
d[0][j] = y[j]
for i in range(1, n):
for j in range(n - i):
d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (x[j+i] - x[j])
print("\nDivided Difference Table:")
header = f"{'x':<10}{'y':<10}"
for i in range(1, n):
header += f"{f'dd{i}':<10}"
print(header)
print("-" * (10 * (n + 1)))
for j in range(n):
row_str = f"{x[j]:<10.4f}{d[0][j]:<10.4f}"
for i in range(1, n - j):
row_str += f"{d[i][j]:<10.4f}"
print(row_str)
equation = f"{d[0][0]}"
for i in range(1, n):
equation += f" + ({d[i][0]})"
for j in range(i):
equation += f"(x - {x[j]})"
print("\nNewton Interpolation Equation:")
print("P(x) =", equation)
x_value = float(input("\nEnter x value: "))
y_interp = d[0][0]
p = 1
for i in range(n - 1):
p = p * (x_value - x[i])
y_interp += p * d[i+1][0]
print(f"\nValue at x = {x_value}: {y_interp}")
OUTPUT
Enter number of data points: 3
Enter x and y values:
x[0] = 1
y[0] = 1
x[1] = 2
y[1] = 4
x[2] = 3
y[2] = 9
Enter x value: 2.5
9.FORWARD DIFF
n = int(input("Enter number of data points: "))
x = []
y = []
print("Enter x and y values:")
for i in range(n):
xi = float(input(f"x[{i}] = "))
yi = float(input(f"y[{i}] = "))
x.append(xi)
y.append(yi)
h = x[1] - x[0]
for i in range(2, n):
if abs((x[i] - x[i-1]) - h) > 1e-10:
print("Error: x values must be equally spaced.")
exit()
d = [[0.0 for j in range(n)] for i in range(n)]
for j in range(n):
d[0][j] = y[j]
for i in range(1, n):
for j in range(n - i):
d[i][j] = d[i-1][j+1] - d[i-1][j]
print("\nForward Difference Table:")
header = f"{'x':<10}{'y':<10}"
for i in range(1, n):
header += f"{f'dd{i}':<10}"
print(header)
print("-" * (10 * (n + 1)))
for j in range(n):
row_str = f"{x[j]:<10.4f}{y[j]:<10.4f}"
for i in range(1, n):
if j >= i:
row_str += f"{d[i][j-i]:<10.4f}"
else:
row_str += f"{'':<10}"
print(row_str)
equation = f"{d[0][0]}"
for i in range(1, n):
equation += f" + ({d[i][0]})"
for j in range(i):
equation += f"(u - {j})"
equation += f"/{i}!"
print("\nNewton Forward Interpolation Equation:")
print("P(x) =", equation)
l = int(input("\nHow many x values do you want to interpolate? "))
for k in range(l):
x_value = float(input(f"\nEnter x value {k+1}: "))
u = (x_value - x[0]) / h
y_interp = d[0][0]
u_product = 1
factorial = 1
for i in range(1, n):
u_product = u_product * (u - (i - 1))
factorial = factorial * i
y_interp += (u_product / factorial) * d[i][0]
print(f"u = (x - x0) / h = ({x_value} - {x[0]}) / {h} = {u}")
print(f"Value at x = {x_value}: {y_interp}")
OUTPUT
Enter number of data points: 3
Enter x and y values:
x[0] = 1
y[0] = 1
x[1] = 2
y[1] = 4
x[2] = 3
y[2] = 9
Enter x value: 2.5
10.BACKWARD DIFF
import matplotlib.pyplot as plt
n = int(input("Enter number of data points: "))
x = []
y = []
print("Enter x and y values:")
for i in range(n):
xi = float(input(f"x[{i}] = "))
yi = float(input(f"y[{i}] = "))
x.append(xi)
y.append(yi)
h = x[1] - x[0]
for i in range(2, n):
if abs((x[i] - x[i-1]) - h) > 1e-10:
print("Error: x values must be equally spaced.")
exit()
d = [[0.0 for j in range(n)] for i in range(n)]
for j in range(n):
d[0][j] = y[j]
for i in range(1, n):
for j in range(n - i):
d[i][j] = d[i-1][j+1] - d[i-1][j]
print("\nBackward Difference Table:")
header = f"{'x':<10}{'y':<10}"
for i in range(1, n):
header += f"{f'dd{i}':<10}"
print(header)
print("-" * (10 * (n + 1)))
for j in range(n):
row_str = f"{x[j]:<10.4f}{y[j]:<10.4f}"
for i in range(1, n):
if j >= i:
row_str += f"{d[i][j-i]:<10.4f}"
else:
row_str += f"{'':<10}"
print(row_str)
equation = f"{d[0][n-1]}"
for i in range(1, n):
equation += f" + ({d[i][n-i-1]})"
for j in range(i):
equation += f"(u + {j})"
equation += f"/{i}!"
print("\nNewton Backward Interpolation Equation:")
print("P(x) =", equation)
l = int(input("\nHow many x values do you want to interpolate? "))
x_inter = []
y_inter = []
for k in range(l):
x_value = float(input(f"\nEnter x value {k+1}: "))
u = (x_value - x[n-1]) / h
y_interp = d[0][n-1]
u_product = 1
factorial = 1
for i in range(1, n):
u_product = u_product * (u + (i - 1))
factorial = factorial * i
y_interp += (u_product / factorial) * d[i][n-i-1]
x_inter.append(x_value)
y_inter.append(y_interp)
print(f"u = (x - xn) / h = ({x_value} - {x[n-1]}) / {h} = {u}")
print(f"Value at x = {x_value}: {y_interp}")
plt.plot(x, y, 'bo-', label='Given Data')
plt.plot(x_inter, y_inter, 'go', label='Interpolated Points')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Newton Backward Interpolation")
plt.legend()
plt.grid(True)
plt.show()
OUTPUT
Enter number of data points: 3
Enter x and y values:
x[0] = 1
y[0] = 1
x[1] = 2
y[1] = 4
x[2] = 3
y[2] = 9
How many x values do you want to interpolate? 2
Enter x value 1: 2.5
u = (x - xn) / h = (2.5 - 3.0) / 1.0 = -0.5
Value at x = 2.5: 6.25
Enter x value 2: 2.8
u = (x - xn) / h = (2.8 - 3.0) / 1.0 = -0.20000000000000018
Value at x = 2.8: 7.84
11.CUBIC SPLINE
import sympy as sp
n = int(input("Enter number of data points: "))
x_data = []
y_data = []
for i in range(n):
x_val, y_val = map(float, input(f"Enter x[{i}] and y[{i}]: ").split())
x_data.append(x_val)
y_data.append(y_val)
l = int(input("Enter number of interpolated values: "))
x_values = []
for k in range(l):
while True:
x_val = float(input(f"Enter interpolation x[{k}]: "))
if x_data[0] <= x_val <= x_data[-1]:
x_values.append(x_val)
break
else:
print("Incorrect interpolated value. Enter again.")
x = sp.Symbol('x')
h=[0]
for i in range (1,n):
h.append(x_data[i]-x_data[i-1])
m = [0] * n
for i in range(1, n - 1):
m[i] = sp.Symbol(f'm{i}')
equation = []
for i in range(1, n - 1):
eq = (
m[i - 1]
+ 4 * m[i]
+ m[i + 1]
- (6 / h[i]**2) *
(y_data[i + 1] - 2 * y_data[i] + y_data[i - 1]))
equation.append(eq)
unknowns = m[1:n-1]
values = sp.solve(equation, unknowns)
for i in range(1, n - 1):
m[i] = values[m[i]]
fx = []
for i in range(1, n):
eq = (
((x_data[i] - x)**3 * m[i - 1]) / (6 * h[i])
+ ((x - x_data[i - 1])**3 * m[i]) / (6 * h[i])
+ ((x_data[i] - x) / h[i]) *
(y_data[i - 1] - (h[i]**2 * m[i - 1]) / 6)
+ ((x - x_data[i - 1]) / h[i]) *
(y_data[i] - (h[i]**2 * m[i]) / 6))
fx.append(sp.expand(eq))
print("\nCubic Spline Functions:")
for i in range(n - 1):
print(f"F{i + 1}(x) =", fx[i])
print("\nInterpolated Values:")
for x_val in x_values:
interval = None
for i in range(n - 1):
if x_data[i] <= x_val <= x_data[i + 1]:
interval = i
break
if interval is not None:
result = fx[interval].subs(x, x_val)
result = float(result)
print(f"S({x_val}) = {result}")
plot=[]
for i in range (n-1):
p=sp.plot(fx[i],(x,x_data[i],x_data[i+1]),show=False,label=f'F{i+1}X')
plot.append(p)
final_plot=plot[0]
for p in plot[1:]:
final_plot.extend(p)
final_plot.show()
OUTPUT
Enter number of data points: 3
Enter x[0] and y[0]: 1 1
Enter x[1] and y[1]: 2 4
Enter x[2] and y[2]: 3 9
Enter number of interpolated values: 2
Enter interpolation x[0]: 1.5
Enter interpolation x[1]: 2.5
12.TRAPEZOIDAL
import sympy as sp
n = int(input("Enter number of subintervals: "))
a = float(input("Enter lower limit a: "))
b = float(input("Enter upper limit b: "))
function_input = input("Enter the function f(x): ")
x = sp.Symbol('x')
f_expr = sp.sympify(function_input)
f = sp.lambdify(x, f_expr, 'math')
h = (b - a) / n
sum_value = f(a) + f(b)
for i in range(1, n):
x_value = a + i * h
sum_value += 2 * f(x_value)
integral = (h / 2) * sum_value
exact_integral = sp.integrate(f_expr, (x, a, b))
exact_value = float(exact_integral)
error = abs(exact_value - integral)
print("\nFunction f(x) =", f_expr)
print("Step size h =", h)
print("Approximate integral =", integral)
print("Exact integral =", exact_value)
print("Error =", error)
OUTPUT
Enter number of subintervals: 4
Enter lower limit a: 0
Enter upper limit b: 2
Enter the function f(x): x**2
13.SIMPSON
import sympy as sp
n=int(input("Number of Sub-Intervals:"))
a=int(input("Start Interval: "))
b=int(input("End Interval: "))
h= (b-a)/n
x = sp.Symbol('x')
f = sp.sympify(input("Enter f(x): "))
area = f.subs(x,a) + f.subs(x,b)
for i in range (1,n):
m= a +(i*h)
if i%2==1:
area += 4*(f.subs(x,m))
else:
area += 2*(f.subs(x,m))
area = area * (h/3)
print("The Area under the Curve (Using Composite Simson's 1/3 Rule): ",area)
integral_area=sp.integrate(f,(x,a,b))
print("The Area under the Curve (Using Integration): ",float(integral_area))
print("Error : ",(float(integral_area)-area))
OUTPUT
Number of Sub-Intervals: 4
Start Interval: 0
End Interval: 2
Enter f(x): x**2
14.QUADRATURE
import sympy as sp
n=int(input("Number of Sub-Intervals: "))
a=int(input("Start Interval: "))
b=int(input("End Interval: "))
h= (b-a)/n
x = sp.Symbol('x')
f = sp.sympify(input("Enter f(x): "))
area = 0
d = h/(2*(3**0.5))
for i in range(n):
xi = a + (i*h)
xj = xi + h
midpoint = (xi+xj)/2
x1 = midpoint - d
x2 = midpoint + d
area += f.subs(x,x1) + f.subs(x,x2)
area = area * (h/2)
print("The Area under the Curve (Using Composite Two Point Gauss Quadrature Rule): ",area)
integral_area=sp.integrate(f,(x,a,b))
print("The Area under the Curve (Using Integration): ",floa