In [1]: import numpy as np
import [Link] as plt
from [Link] import interp1d
x1=[Link](1,5,10)
print(x1)
y1=x1**3
print(y1)
#y_f=interp1d(x1,y1, 'cubic')
#print(y_f(2.6))
[Link](x1,y1)
[Link]('X_Data')
[Link]('Y_Data')
#[Link]()
[1. 1.44444444 1.88888889 2.33333333 2.77777778 3.22222222
3.66666667 4.11111111 4.55555556 5. ]
[ 1. 3.01371742 6.739369 12.7037037 21.43347051
33.45541838 49.2962963 69.48285322 94.54183813 125. ]
Out[1]: Text(0, 0.5, 'Y_Data')
In [2]: import numpy as np
import [Link] as plt
from [Link] import interp1d
x2=[Link](1,5,10)
y2=x**3
[Link](x2,y2,'o')
[Link]()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 5
3 from [Link] import interp1d
4 x2=[Link](1,5,10)
----> 5 y2=x**3
7 [Link](x2,y2,'o')
8 [Link]()
NameError: name 'x' is not defined
In [144… import numpy as np
import [Link] as plt
from [Link] import interp1d
x=[Link]([1,6,7,9,12,20])
print("x=",x)
y=[Link]([2,8,6,12,14,41])
print("y=",y)
#[Link](x,y,"o", label="Data points")
#[Link]()
[Link](x,y)
[Link]()
x= [ 1 6 7 9 12 20]
y= [ 2 8 6 12 14 41]
In [148… import numpy as np
import [Link] as plt
from [Link] import interp1d
x=[Link]([1,6,7,9,12,20])
y=[Link]([2,8,6,12,14,41])
[Link](x,y,"o", label="Data points")
#print(x_interp)
x_interp=[Link]([Link](x),[Link](x),50)
#1D line spline interpolation
y_linear=interp1d(x,y)
#print(y_linear(3))
[Link](x_interp, y_linear(x_interp), "red", label="Linear Spline" )
[Link]()
[Link]()
In [150… import numpy as np
import [Link] as plt
from [Link] import interp1d
x=[Link]([1,6,7,9,12,20])
y=[Link]([2,8,6,12,14,41])
[Link](x,y,"o", label="Data points")
x_interp=[Link]([Link](x),[Link](x),50)
#1D quadratic spline interpolation
y_quadratic=interp1d(x,y,kind="quadratic")
[Link](x_interp, y_quadratic(x_interp), "green", label="Quadratic Spline" )
[Link]()
[Link]()
In [152… import numpy as np
import [Link] as plt
from [Link] import interp1d
x=[Link]([1,6,7,9,12,20])
y=[Link]([2,8,6,12,14,41])
x_interp=[Link]([Link](x),[Link](x),50)
#1D Cubic spline interpolation
y_cubic=interp1d(x, y, kind="cubic")
[Link](x,y,"o", label="Data points")
[Link](x_interp, y_cubic(x_interp), "green", label="Cubic Spline" )
[Link]()
[Link]()
In [57]: import numpy as np
import [Link] as plt
from [Link] import interp1d, CubicSpline
x=[Link]([1,6,7,9,12,20])
y=[Link]([2,8,6,12,14,41])
x_interp=[Link]([Link](x),[Link](x),50)
#1D Cubic spline interpolation Boundary Condition
y_cubicBC=CubicSpline(x, y, bc_type="natural")
[Link](x,y,"o", label="Data points")
[Link](x_interp, y_cubicBC(x_interp), "black", label="Cubic Spline BC" )
[Link]()
[Link]()
In [4]: # Newton-Raphson Method in Python
# Definition-based implementation
def newton_raphson_method():
# Define the function f(x)
def f(x):
return x**3 - 5*x + 3 # Example: f(x) = x^3 - 5x + 3
# Define the derivative f'(x)
def df(x):
return 3*x**2 - 5 # Derivative: f'(x) = 3x^2 - 5
# Initial guess (you can change this)
x0 = 1.0
# Tolerance (how close the result should be)
tol = 1e-6
# Maximum number of iterations
max_iter = 100
print("Newton–Raphson Method")
print("-----------------------")
print(f"Initial guess: x0 = {x0}")
print("Iteration process:")
for i in range(max_iter):
fx = f(x0)
dfx = df(x0)
# If derivative is zero, stop (to avoid division by zero)
if dfx == 0:
print("Error: Derivative is zero. Cannot continue.")
return None
# Newton-Raphson formula
x1 = x0 - fx / dfx
print(f"Iteration {i+1}: x = {x1:.6f}, f(x) = {f(x1):.6f}")
# Check convergence
if abs(x1 - x0) < tol:
print("\nConverged!")
print(f"Approximate root = {x1:.6f}")
return x1
# Update value
x0 = x1
print("\nDid not converge within the maximum iterations.")
return None
# Run the method
newton_raphson_method()
Newton–Raphson Method
-----------------------
Initial guess: x0 = 1.0
Iteration process:
Iteration 1: x = 0.500000, f(x) = 0.625000
Iteration 2: x = 0.647059, f(x) = 0.035620
Iteration 3: x = 0.656573, f(x) = 0.000177
Iteration 4: x = 0.656620, f(x) = 0.000000
Iteration 5: x = 0.656620, f(x) = 0.000000
Converged!
Approximate root = 0.656620
Out[4]: 0.6566204310471104
In [6]: # Bisection Method (step-by-step)
# Step 1: Define the function
def f(x):
return x**3 - 4*x - 9 # Example function; change as needed
# Step 2: Set the initial interval and tolerance
a = 2
b = 3
tolerance = 0.001
# Step 3: Check if a root exists between a and b
if f(a) * f(b) > 0:
print("Root may not exist in this interval. f(a) and f(b) should have opposite sig
else:
iteration = 1
print("Iter\t a\t\t b\t\t c\t\t f(c)")
# Step 4: Apply the bisection formula repeatedly
while abs(b - a) >= tolerance:
c = (a + b) / 2 # Formula: midpoint of interval
print(f"{iteration}\t {a:.6f}\t {b:.6f}\t {c:.6f}\t {f(c):.6f}")
# Step 5: Check which subinterval contains the root
if f(c) == 0:
break
elif f(a) * f(c) < 0:
b = c
else:
a = c
iteration += 1
# Step 6: Display the final result
print(f"\nApproximate root after {iteration-1} iterations = {c:.6f}")
Iter a b c f(c)
1 2.000000 3.000000 2.500000 -3.375000
2 2.500000 3.000000 2.750000 0.796875
3 2.500000 2.750000 2.625000 -1.412109
4 2.625000 2.750000 2.687500 -0.339111
5 2.687500 2.750000 2.718750 0.220917
6 2.687500 2.718750 2.703125 -0.061077
7 2.703125 2.718750 2.710938 0.079423
8 2.703125 2.710938 2.707031 0.009049
9 2.703125 2.707031 2.705078 -0.026045
10 2.705078 2.707031 2.706055 -0.008506
Approximate root after 10 iterations = 2.706055
In [8]: import math
def jacobi(A, b, x0=None, tol=1e-8, maxiter=1000):
n = len(A)
if x0 is None:
x = [0.0]*n
else:
x = x0[:]
x_new = [0.0]*n
for k in range(1, maxiter+1):
for i in range(n):
s = 0.0
for j in range(n):
if j != i:
s += A[i][j] * x[j]
x_new[i] = (b[i] - s) / A[i][i]
# compute infinity-norm of difference
diff = max(abs(x_new[i] - x[i]) for i in range(n))
# copy new -> x
x[:] = x_new[:]
if diff < tol:
return x, k, diff
return x, maxiter, diff
# Example
A = [[4.0, -1.0, 0.0],
[-1.0, 4.0, -1.0],
[0.0, -1.0, 3.0]]
b = [15.0, 10.0, 10.0]
x, iters, final_err = jacobi(A, b, tol=1e-10, maxiter=200)
print("solution:", x)
print("iterations:", iters, "final_error:", final_err)
solution: [4.999999999983133, 4.999999999966265, 4.99999999997751]
iterations: 27 final_error: 5.46185319194592e-11
In [10]: import math
def gauss_seidel(A, b, x0=None, tol=1e-8, maxiter=1000):
n = len(A)
if x0 is None:
x = [0.0]*n
else:
x = x0[:]
for k in range(1, maxiter + 1):
x_old = x[:]
for i in range(n):
s1 = sum(A[i][j] * x[j] for j in range(i)) # new values
s2 = sum(A[i][j] * x_old[j] for j in range(i+1, n)) # old values
x[i] = (b[i] - s1 - s2) / A[i][i]
# check convergence
diff = max(abs(x[i] - x_old[i]) for i in range(n))
if diff < tol:
return x, k, diff
return x, maxiter, diff
# Example system
A = [[4.0, -1.0, 0.0],
[-1.0, 4.0, -1.0],
[0.0, -1.0, 3.0]]
b = [15.0, 10.0, 10.0]
x, iters, final_err = gauss_seidel(A, b, t#ol=1e-10, maxiter=200)
print("solution:", x)
print("iterations:", iters, "final_error:", final_err)
solution: [4.999999999994729, 4.999999999996925, 4.999999999998975]
iterations: 15 final_error: 3.087308186877635e-11
In [12]: # SV Deco Code
import numpy as np
from [Link] import svd
In [14]: #This function calculates the Eigenvectors corresponding for U matrice
def calculU(M):
B = [Link](M, M.T)
eigenvalues, eigenvectors = [Link](B)
ncols = [Link](eigenvalues)[::-1]
return eigenvectors[:,ncols]
In [16]: #This function calculates the Eigenvectors corresponding for V matrice
def calculVt(M):
B = [Link](M.T, M)
eigenvalues, eigenvectors = [Link](B)
ncols = [Link](eigenvalues)[::-1]
return eigenvectors[:,ncols].T
In [18]: #Function that calculates Eigenvalues corresponding to the Sigma Matrix
def calculSigma(M):
if ([Link]([Link](M, M.T)) > [Link]([Link](M.T, M))):
newM = [Link](M.T, M)
else:
newM = [Link](M, M.T)
eigenvalues, eigenvectors = [Link](newM)
eigenvalues = [Link](eigenvalues)
#Sorting in descending order as the svd function does
return eigenvalues[::-1]
In [20]: #Creating our matrix
A = [Link]([[1,1],[[Link](3),0],[0,[Link](3)]])
In [36]: print("-------------------U-------------------")
print(U)
print("\n--------------Sigma----------------")
print(Sigma)
print("\n-------------V transpose---------------")
print(Vt)
-------------------U-------------------
[[ 0.28978415 -0.95709203]
[ 0.95709203 0.28978415]]
--------------Sigma----------------
[8.13872588 3.97003036]
-------------V transpose---------------
[[ 0.26001965 0.65919758 0.70558368]
[-0.89132415 -0.11719389 0.43795758]
[ 0.37139068 -0.74278135 0.55708601]]
In [ ]: newSigma = [Link]((2, 3)) # define acoor. as matrix A to check [Link]
newSigma[:2, :2] = [Link](Sigma[:2])
print(A,"\n")
A_remake = (U @ newSigma @ Vt)
print(A_remake)
In [22]: #Calling the corresponding Fuctions and saving the values in variables
U = calculU(A)
Sigma = calculSigma(A)
Vt = calculVt(A)
In [38]: #Using Direct Pack method
In [30]: #Creating our matrix
A = [Link]([[4,2,0],[1,5,6]])
In [32]: U, Sigma, Vt = [Link](A)
In [40]: # Convert Sigma (1D array) → diagonal matrix
Sigma1 = [Link](([Link][0], [Link][0]))
np.fill_diagonal(Sigma1, Sigma)
In [42]: # LD Deco Method code
In [44]: #step by step run the code
In [ ]: def lu_decomposition(A, n):
"""Performs LU Decomposition of matrix A (n x n) into L and U"""
L = [[0.0] * n for _ in range(n)]
U = [[0.0] * n for _ in range(n)]
print(L,U)
for i in range(n):
# Upper Triangular
for k in range(i, n):
sum_val = 0
for j in range(i):
sum_val += (L[i][j] * U[j][k])
U[i][k] = A[i][k] - sum_val
# Lower Triangular
L[i][i] = 1 # diagonal as 1
for k in range(i+1, n):
sum_val = 0
for j in range(i):
sum_val += (L[k][j] * U[j][i])
L[k][i] = (A[k][i] - sum_val) / U[i][i]
return L, U
In [46]: A = [[1, 1, 1],
[4, 3, -1],
[3, 5, 3]]
In [ ]: # Now call the function To calculate the L and U Matrix
lu_decomposition(A,3)
In [50]: # Solving system of equation using SVD decomposition
In [ ]: # Calculating variable values x
In [ ]: def forward_substitution(L, b, n):
"""Solve Ly = b"""
y = [0.0] * n
for i in range(n):
sum_val = 0
for j in range(i):
sum_val += L[i][j] * y[j]
y[i] = b[i] - sum_val
return y
def backward_substitution(U, y, n):
"""Solve Ux = y"""
x = [0.0] * n
for i in range(n-1, -1, -1):
sum_val = 0
for j in range(i+1, n):
sum_val += U[i][j] * x[j]
x[i] = (y[i] - sum_val) / U[i][i]
return x
In [ ]: #Insert matrix and print the LU matrix and Values of variable X
In [ ]: # System definition
A = [[1, 1, 1],
[4, 3, -1],
[3, 5, 3]]
b = [1, 6, 4]
n = 3
# LU Decomposition
L, U = lu_decomposition(A, n)
# Solve Ly = b
y = forward_substitution(L, b, n)
# Solve Ux = y
x = backward_substitution(U, y, n)
print("Solution: x1 = {:.2f}, x2 = {:.2f}, x3 = {:.2f}".format(x[0], x[1], x[2]))
In [ ]: