Scribd
Upload
6 views15 pages

Codess

The document provides various numerical methods for solving linear equations, interpolation, and differential equations. It includes implementations of the Gauss elimination method, Jacobi method, Gauss-Seidel method, LU decomposition, Cholesky method, and several root-finding algorithms. Additionally, it covers numerical integration techniques such as the trapezoidal rule and Simpson's rule, as well as interpolation methods like Lagrange and Newton's divided difference.

Uploaded by

princesharma200503
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
6 views15 pages

Codess

The document provides various numerical methods for solving linear equations, interpolation, and differential equations. It includes implementations of the Gauss elimination method, Jacobi method, Gauss-Seidel method, LU decomposition, Cholesky method, and several root-finding algorithms. Additionally, it covers numerical integration techniques such as the trapezoidal rule and Simpson's rule, as well as interpolation methods like Lagrange and Newton's divided difference.

Uploaded by

princesharma200503
Copyright
© All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

📌 Gauss Elimination Method :

def gauss_elimination(A, b):


n = len(b)

# Forward Elimination
for i in range(n):
for j in range(i + 1, n):
ratio = A[j][i] / A[i][i]

for k in range(n):
A[j][k] = A[j][k] - ratio * A[i][k]

b[j] = b[j] - ratio * b[i]

# Back Substitution
x = [0] * n

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


x[i] = b[i]

for j in range(i + 1, n):


x[i] = x[i] - A[i][j] * x[j]

x[i] = x[i] / A[i][i]

return x

# Example
A=[
[2, 1, -1],
[-3, -1, 2],
[-2, 1, 2]
]

b = [8, -11, -3]

result = gauss_elimination(A, b)
print("Solution:", result)
📌 Jacobi Method
def jacobi(A, b, x, iterations):
n = len(A)

for it in range(iterations):
x_new = [0] * n

for i in range(n):
s=0
for j in range(n):
if j != i:
s += A[i][j] * x[j]

x_new[i] = (b[i] - s) / A[i][i]

x = x_new

return x

# Example
A=[
[10, -1, 2],
[-1, 11, -1],
[2, -1, 10]
]

b = [6, 25, -11]


x0 = [0, 0, 0]

result = jacobi(A, b, x0, 10)


print("Solution:", result)

📌 Gauss-Seidel Method
def gauss_seidel(A, b, x, iterations):
n = len(A)

for it in range(iterations):
for i in range(n):
s=0

for j in range(n):
if j != i:
s += A[i][j] * x[j]

x[i] = (b[i] - s) / A[i][i]

return x
📌 1. LU Decomposition (Doolittle Method)

def doolittle(A):

n = len(A)

L = [[0]*n for _ in range(n)]

U = [[0]*n for _ in range(n)]

for i in range(n):

L[i][i] = 1 # diagonal = 1

for i in range(n):

# U matrix

for j in range(i, n):

s=0

for k in range(i):

s += L[i][k] * U[k][j]

U[i][j] = A[i][j] - s

# L matrix

for j in range(i+1, n):

s=0

for k in range(i):

s += L[j][k] * U[k][i]

L[j][i] = (A[j][i] - s) / U[i][i]

return L, U

# Example

Take matrix A : L , U = doolittle(A)

print("L:", L) , print(“U:”, U)
📌 2. LU Decomposition (Crout Method)
def crout(A):

n = len(A)

L = [[0]*n for _ in range(n)]

U = [[0]*n for _ in range(n)]

for j in range(n):

U[j][j] = 1 # diagonal = 1

for j in range(n):

# L matrix

for i in range(j, n):

s=0

for k in range(j):

s += L[i][k] * U[k][j]

L[i][j] = A[i][j] - s

# U matrix

for i in range(j+1, n):

s=0

for k in range(j):

s += L[j][k] * U[k][i]

U[j][i] = (A[j][i] - s) / L[j][j]

return L, U
📌 3. Cholesky Method

import math

def cholesky(A):

n = len(A)

L = [[0]*n for _ in range(n)]

for i in range(n):

for j in range(i+1):

s=0

for k in range(j):

s += L[i][k] * L[j][k]

if i == j:

L[i][j] = [Link](A[i][i] - s)

else:

L[i][j] = (A[i][j] - s) / L[j][j]

return L

📌 SVD (Singular Value Decomposition)

def svd_decomposition(A):

U, S, VT = [Link](A)

return U, S, VT

# Example
Take A ;

U, S, VT = svd_decomposition(A)

print("U =\n", U) , print(“Singular Values =\n", S) , print(“VT =\n", VT)


📌 1. Incremental Search Method
def incremental_search(f, a, b, h):
x=a

while x < b:
if f(x) * f(x + h) < 0:
print("Root between", x, "and", x + h)
x=x+h

# Example
f = lambda x: x**3 - x - 2
incremental_search(f, 1, 2, 0.1)

📌 2. Bisection Method
def bisection(f, a, b, tol):

if f(a) * f(b) >= 0:

print("Invalid interval")

return

while (b - a) > tol:

c = (a + b) / 2

if f(c) == 0:

return c

elif f(a) * f(c) < 0:

b=c

else:

a=c

return (a + b) / 2

# Example

f = lambda x: x**3 - x - 2

root = bisection(f, 1, 2, 0.0001)

print("Root:", root)
📌 3. Newton–Raphson Method
def newton_raphson(f, df, x0, tol):

x = x0

while True:

x_new = x - f(x) / df(x)

if abs(x_new - x) < tol:

return x_new

x = x_new

📌 Regular Falsi Method

def regula_falsi(f, a, b, tol):

if f(a) * f(b) >= 0:

print("Invalid interval")

return

while True:

c = (a*f(b) - b*f(a)) / (f(b) - f(a))

if abs(f(c)) < tol:

return c

elif f(a) * f(c) < 0:

b=c

else:

a=c

# Example

f = lambda x: x**3 - x - 2 , root = regula_falsi(f, 1, 2, 0.0001) , print(“Root:", root)


📌 1. EULER
import numpy as np

def f(x, y):


return (same expression)

x=0
y=1
h = 0.1

for i in range(4):
y = y + h * f(x, y)
x=x+h

print("Euler:", y)

📌 2. MODIFIED EULER
import numpy as np

def f(x, y):


return (same expression)

x=0
y=1
h = 0.1

for i in range(4):
y_star = y + h * f(x, y)
y = y + (h/2) * (f(x, y) + f(x+h, y_star))
x=x+h

print("Modi ed Euler:", y)

📌 3. RK2
import numpy as np

def f(x, y):


return (same expression)

x=0
y=1
h = 0.1

for i in range(4):
k1 = h * f(x, y)
k2 = h * f(x+h, y+k1)
y = y + 0.5 * (k1 + k2)
x=x+h

print("RK2:", y)
fi
📌 4. RK4
import numpy as np

def f(x, y):


return (same expression)

x=0
y=1
h = 0.1

for i in range(4):
k1 = h * f(x, y)
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)

y = y + (k1 + 2*k2 + 2*k3 + k4)/6


x=x+h

print("RK4:", y)
👉 Common function
import numpy as np

def f(x):
return [Link](1 + x**2) # change this

a=0
b=3
n=6

h = (b - a)/n
x = [Link](a, b, n+1)
y = f(x)

📌 TRAPEZOIDAL
I = h/2 * (y[0] + 2*sum(y[1:n]) + y[n])
print("Trapezoidal:", I)

📌 SIMPSON 1/3
I = y[0] + y[n]

for i in range(1, n):


if i % 2 == 0:
I += 2*y[i]
else:
I += 4*y[I]

I = (h/3)*I
print("Simpson 1/3:", I)

📌 SIMPSON 3/8
I = y[0] + y[n]

for i in range(1, n):


if i % 3 == 0:
I += 2*y[i]
else:
I += 3*y[i]

I = (3*h/8)*I
print("Simpson 3/8:", I)
📌 Lagrange Interpolation

def lagrange(x_points, y_points, x):

n = len(x_points)

result = 0

for i in range(n):
term = y_points[i]

for j in range(n):
if i != j:
term *= (x - x_points[j]) / (x_points[i] - x_points[j])

result += term

return result

# Example

x_points = [1, 2, 3]
y_points = [2, 3, 5]

x = 2.5

print("Interpolated value:", lagrange(x_points, y_points, x))


📌 Newton divided di erence Interpolation
def newton_divided(x, y, value):
n = len(x)

# Create divided di erence table


table = [[0 for _ in range(n)] for _ in range(n)]

# First column = y values


for i in range(n):
table[i][0] = y[i]

# Fill table
for j in range(1, n):
for i in range(n - j):
table[i][j] = (table[i+1][j-1] - table[i][j-1]) / (x[i+j] - x[i])

# Evaluate polynomial
result = table[0][0]
product = 1

for i in range(1, n):


product *= (value - x[i-1])
result += table[0][i] * product

return result

# Example
x = [1, 2, 3]
y = [2, 3, 5]

print("Interpolated value:", newton_divided(x, y, 2.5))


ff
ff
📌 Forward Interpolation Code
def newton_forward(x, y, value):
n = len(x)
h = x[1] - x[0]

# Di erence table
di = [[Link]()]
for i in range(1, n):
temp = []
for j in range(n - i):
[Link](di [i-1][j+1] - di [i-1][j])
di .append(temp)

u = (value - x[0]) / h
result = y[0]
u_term = 1

for i in range(1, n):


u_term *= (u - (i-1))
result += (u_term * di [i][0]) / factorial(i)

return result

from math import factorial

# Example
x = [1, 2, 3, 4]
y = [1, 4, 9, 16]

print(newton_forward(x, y, 2.5))
ff
ff
ff
ff
ff
ff
📌 Backward Interpolation Code
def newton_backward(x, y, value):

n = len(x)

h = x[1] - x[0]

# Di erence table

di = [[Link]()]

for i in range(1, n):

temp = []

for j in range(n - i):

[Link](di [i-1][j+1] - di [i-1][j])

di .append(temp)

u = (value - x[-1]) / h

result = y[-1]

u_term = 1

for i in range(1, n):

u_term *= (u + (i-1))

result += (u_term * di [i][-1]) / factorial(i)

return result

from math import factorial

# Example

x = [1, 2, 3, 4]

y = [1, 4, 9, 16]

print(newton_backward(x, y, 3.5))
ff
ff
ff
ff
ff
ff
📌 Cubic spline
import numpy as np
from [Link] import CubicSpline

# Given data points


x = [Link]([1, 2, 3, 4])
y = [Link]([1, 4, 9, 16])

# Create spline
cs = CubicSpline(x, y, bc_type='natural')

# Find value
print(cs(2.5))

You might also like