diff --git a/arithmetic_analysis/gaussian_elimination.py b/arithmetic_analysis/gaussian_elimination.py
index 2dada4fbf..89ed3b323 100644
--- a/arithmetic_analysis/gaussian_elimination.py
+++ b/arithmetic_analysis/gaussian_elimination.py
@@ -5,9 +5,13 @@ Gaussian elimination - https://en.wikipedia.org/wiki/Gaussian_elimination
 
 
 import numpy as np
+from numpy import float64
+from numpy.typing import NDArray
 
 
-def retroactive_resolution(coefficients: np.matrix, vector: np.ndarray) -> np.ndarray:
+def retroactive_resolution(
+    coefficients: NDArray[float64], vector: NDArray[float64]
+) -> NDArray[float64]:
     """
     This function performs a retroactive linear system resolution
         for triangular matrix
@@ -27,7 +31,7 @@ def retroactive_resolution(coefficients: np.matrix, vector: np.ndarray) -> np.nd
 
     rows, columns = np.shape(coefficients)
 
-    x = np.zeros((rows, 1), dtype=float)
+    x: NDArray[float64] = np.zeros((rows, 1), dtype=float)
     for row in reversed(range(rows)):
         sum = 0
         for col in range(row + 1, columns):
@@ -38,7 +42,9 @@ def retroactive_resolution(coefficients: np.matrix, vector: np.ndarray) -> np.nd
     return x
 
 
-def gaussian_elimination(coefficients: np.matrix, vector: np.ndarray) -> np.ndarray:
+def gaussian_elimination(
+    coefficients: NDArray[float64], vector: NDArray[float64]
+) -> NDArray[float64]:
     """
     This function performs Gaussian elimination method
 
@@ -60,7 +66,7 @@ def gaussian_elimination(coefficients: np.matrix, vector: np.ndarray) -> np.ndar
         return np.array((), dtype=float)
 
     # augmented matrix
-    augmented_mat = np.concatenate((coefficients, vector), axis=1)
+    augmented_mat: NDArray[float64] = np.concatenate((coefficients, vector), axis=1)
     augmented_mat = augmented_mat.astype("float64")
 
     # scale the matrix leaving it triangular
diff --git a/arithmetic_analysis/in_static_equilibrium.py b/arithmetic_analysis/in_static_equilibrium.py
index ed0d1eb98..d762a376f 100644
--- a/arithmetic_analysis/in_static_equilibrium.py
+++ b/arithmetic_analysis/in_static_equilibrium.py
@@ -3,7 +3,8 @@ Checks if a system of forces is in static equilibrium.
 """
 from __future__ import annotations
 
-from numpy import array, cos, cross, ndarray, radians, sin
+from numpy import array, cos, cross, float64, radians, sin
+from numpy.typing import NDArray
 
 
 def polar_force(
@@ -27,7 +28,7 @@ def polar_force(
 
 
 def in_static_equilibrium(
-    forces: ndarray, location: ndarray, eps: float = 10**-1
+    forces: NDArray[float64], location: NDArray[float64], eps: float = 10**-1
 ) -> bool:
     """
     Check if a system is in equilibrium.
@@ -46,7 +47,7 @@ def in_static_equilibrium(
     False
     """
     # summation of moments is zero
-    moments: ndarray = cross(location, forces)
+    moments: NDArray[float64] = cross(location, forces)
     sum_moments: float = sum(moments)
     return abs(sum_moments) < eps
 
@@ -61,7 +62,7 @@ if __name__ == "__main__":
         ]
     )
 
-    location = array([[0, 0], [0, 0], [0, 0]])
+    location: NDArray[float64] = array([[0, 0], [0, 0], [0, 0]])
 
     assert in_static_equilibrium(forces, location)
 
diff --git a/arithmetic_analysis/jacobi_iteration_method.py b/arithmetic_analysis/jacobi_iteration_method.py
index 667482425..4336aaa91 100644
--- a/arithmetic_analysis/jacobi_iteration_method.py
+++ b/arithmetic_analysis/jacobi_iteration_method.py
@@ -4,13 +4,15 @@ Jacobi Iteration Method - https://en.wikipedia.org/wiki/Jacobi_method
 from __future__ import annotations
 
 import numpy as np
+from numpy import float64
+from numpy.typing import NDArray
 
 
 # Method to find solution of system of linear equations
 def jacobi_iteration_method(
-    coefficient_matrix: np.ndarray,
-    constant_matrix: np.ndarray,
-    init_val: list,
+    coefficient_matrix: NDArray[float64],
+    constant_matrix: NDArray[float64],
+    init_val: list[int],
     iterations: int,
 ) -> list[float]:
     """
@@ -99,7 +101,9 @@ def jacobi_iteration_method(
     if iterations <= 0:
         raise ValueError("Iterations must be at least 1")
 
-    table = np.concatenate((coefficient_matrix, constant_matrix), axis=1)
+    table: NDArray[float64] = np.concatenate(
+        (coefficient_matrix, constant_matrix), axis=1
+    )
 
     rows, cols = table.shape
 
@@ -125,7 +129,7 @@ def jacobi_iteration_method(
 
 
 # Checks if the given matrix is strictly diagonally dominant
-def strictly_diagonally_dominant(table: np.ndarray) -> bool:
+def strictly_diagonally_dominant(table: NDArray[float64]) -> bool:
     """
     >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 4, -4]])
     >>> strictly_diagonally_dominant(table)
diff --git a/arithmetic_analysis/lu_decomposition.py b/arithmetic_analysis/lu_decomposition.py
index b488b1bb3..371f7b166 100644
--- a/arithmetic_analysis/lu_decomposition.py
+++ b/arithmetic_analysis/lu_decomposition.py
@@ -6,9 +6,13 @@ Reference:
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as NDArray
+from numpy import float64
 
 
-def lower_upper_decomposition(table: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+def lower_upper_decomposition(
+    table: NDArray[float64],
+) -> tuple[NDArray[float64], NDArray[float64]]:
     """Lower-Upper (LU) Decomposition
 
     Example: