diff --git a/DIRECTORY.md b/DIRECTORY.md
index c07e1550..1e3711fe 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -1,17 +1,4 @@
 
-## Arithmetic Analysis
-  * [Bisection](arithmetic_analysis/bisection.py)
-  * [Gaussian Elimination](arithmetic_analysis/gaussian_elimination.py)
-  * [In Static Equilibrium](arithmetic_analysis/in_static_equilibrium.py)
-  * [Intersection](arithmetic_analysis/intersection.py)
-  * [Jacobi Iteration Method](arithmetic_analysis/jacobi_iteration_method.py)
-  * [Lu Decomposition](arithmetic_analysis/lu_decomposition.py)
-  * [Newton Forward Interpolation](arithmetic_analysis/newton_forward_interpolation.py)
-  * [Newton Method](arithmetic_analysis/newton_method.py)
-  * [Newton Raphson](arithmetic_analysis/newton_raphson.py)
-  * [Newton Raphson New](arithmetic_analysis/newton_raphson_new.py)
-  * [Secant Method](arithmetic_analysis/secant_method.py)
-
 ## Audio Filters
   * [Butterworth Filter](audio_filters/butterworth_filter.py)
   * [Iir Filter](audio_filters/iir_filter.py)
@@ -520,6 +507,9 @@
     * [Test Knapsack](knapsack/tests/test_knapsack.py)
 
 ## Linear Algebra
+  * [Gaussian Elimination](linear_algebra/gaussian_elimination.py)
+  * [Jacobi Iteration Method](linear_algebra/jacobi_iteration_method.py)
+  * [Lu Decomposition](linear_algebra/lu_decomposition.py)
   * Src
     * [Conjugate Gradient](linear_algebra/src/conjugate_gradient.py)
     * [Lib](linear_algebra/src/lib.py)
@@ -583,7 +573,6 @@
   * [Binary Multiplication](maths/binary_multiplication.py)
   * [Binomial Coefficient](maths/binomial_coefficient.py)
   * [Binomial Distribution](maths/binomial_distribution.py)
-  * [Bisection](maths/bisection.py)
   * [Ceil](maths/ceil.py)
   * [Chebyshev Distance](maths/chebyshev_distance.py)
   * [Check Polygon](maths/check_polygon.py)
@@ -617,7 +606,6 @@
   * [Germain Primes](maths/germain_primes.py)
   * [Greatest Common Divisor](maths/greatest_common_divisor.py)
   * [Hardy Ramanujanalgo](maths/hardy_ramanujanalgo.py)
-  * [Integration By Simpson Approx](maths/integration_by_simpson_approx.py)
   * [Interquartile Range](maths/interquartile_range.py)
   * [Is Int Palindrome](maths/is_int_palindrome.py)
   * [Is Ip V4 Address Valid](maths/is_ip_v4_address_valid.py)
@@ -644,10 +632,24 @@
   * [Modular Exponential](maths/modular_exponential.py)
   * [Monte Carlo](maths/monte_carlo.py)
   * [Monte Carlo Dice](maths/monte_carlo_dice.py)
-  * [Nevilles Method](maths/nevilles_method.py)
-  * [Newton Raphson](maths/newton_raphson.py)
   * [Number Of Digits](maths/number_of_digits.py)
-  * [Numerical Integration](maths/numerical_integration.py)
+  * Numerical Analysis
+    * [Bisection](maths/numerical_analysis/bisection.py)
+    * [Bisection 2](maths/numerical_analysis/bisection_2.py)
+    * [Integration By Simpson Approx](maths/numerical_analysis/integration_by_simpson_approx.py)
+    * [Intersection](maths/numerical_analysis/intersection.py)
+    * [Nevilles Method](maths/numerical_analysis/nevilles_method.py)
+    * [Newton Forward Interpolation](maths/numerical_analysis/newton_forward_interpolation.py)
+    * [Newton Method](maths/numerical_analysis/newton_method.py)
+    * [Newton Raphson](maths/numerical_analysis/newton_raphson.py)
+    * [Newton Raphson 2](maths/numerical_analysis/newton_raphson_2.py)
+    * [Newton Raphson New](maths/numerical_analysis/newton_raphson_new.py)
+    * [Numerical Integration](maths/numerical_analysis/numerical_integration.py)
+    * [Runge Kutta](maths/numerical_analysis/runge_kutta.py)
+    * [Runge Kutta Fehlberg 45](maths/numerical_analysis/runge_kutta_fehlberg_45.py)
+    * [Secant Method](maths/numerical_analysis/secant_method.py)
+    * [Simpson Rule](maths/numerical_analysis/simpson_rule.py)
+    * [Square Root](maths/numerical_analysis/square_root.py)
   * [Odd Sieve](maths/odd_sieve.py)
   * [Perfect Cube](maths/perfect_cube.py)
   * [Perfect Number](maths/perfect_number.py)
@@ -673,8 +675,6 @@
   * [Radians](maths/radians.py)
   * [Radix2 Fft](maths/radix2_fft.py)
   * [Remove Digit](maths/remove_digit.py)
-  * [Runge Kutta](maths/runge_kutta.py)
-  * [Runge Kutta Fehlberg 45](maths/runge_kutta_fehlberg_45.py)
   * [Segmented Sieve](maths/segmented_sieve.py)
   * Series
     * [Arithmetic](maths/series/arithmetic.py)
@@ -687,7 +687,6 @@
   * [Sieve Of Eratosthenes](maths/sieve_of_eratosthenes.py)
   * [Sigmoid](maths/sigmoid.py)
   * [Signum](maths/signum.py)
-  * [Simpson Rule](maths/simpson_rule.py)
   * [Simultaneous Linear Equation Solver](maths/simultaneous_linear_equation_solver.py)
   * [Sin](maths/sin.py)
   * [Sock Merchant](maths/sock_merchant.py)
@@ -709,7 +708,6 @@
     * [Proth Number](maths/special_numbers/proth_number.py)
     * [Ugly Numbers](maths/special_numbers/ugly_numbers.py)
     * [Weird Number](maths/special_numbers/weird_number.py)
-  * [Square Root](maths/square_root.py)
   * [Sum Of Arithmetic Series](maths/sum_of_arithmetic_series.py)
   * [Sum Of Digits](maths/sum_of_digits.py)
   * [Sum Of Geometric Progression](maths/sum_of_geometric_progression.py)
@@ -812,6 +810,7 @@
   * [Horizontal Projectile Motion](physics/horizontal_projectile_motion.py)
   * [Hubble Parameter](physics/hubble_parameter.py)
   * [Ideal Gas Law](physics/ideal_gas_law.py)
+  * [In Static Equilibrium](physics/in_static_equilibrium.py)
   * [Kinetic Energy](physics/kinetic_energy.py)
   * [Lorentz Transformation Four Vector](physics/lorentz_transformation_four_vector.py)
   * [Malus Law](physics/malus_law.py)
diff --git a/arithmetic_analysis/README.md b/arithmetic_analysis/README.md
deleted file mode 100644
index 45cf321e..00000000
--- a/arithmetic_analysis/README.md
+++ /dev/null
@@ -1,7 +0,0 @@
-# Arithmetic analysis
-
-Arithmetic analysis is a branch of mathematics that deals with solving linear equations.
-
-* <https://en.wikipedia.org/wiki/System_of_linear_equations>
-* <https://en.wikipedia.org/wiki/Gaussian_elimination>
-* <https://en.wikipedia.org/wiki/Root-finding_algorithms>
diff --git a/arithmetic_analysis/image_data/__init__.py b/arithmetic_analysis/image_data/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/arithmetic_analysis/gaussian_elimination.py b/linear_algebra/gaussian_elimination.py
similarity index 100%
rename from arithmetic_analysis/gaussian_elimination.py
rename to linear_algebra/gaussian_elimination.py
diff --git a/arithmetic_analysis/jacobi_iteration_method.py b/linear_algebra/jacobi_iteration_method.py
similarity index 96%
rename from arithmetic_analysis/jacobi_iteration_method.py
rename to linear_algebra/jacobi_iteration_method.py
index 44c52dd4..8c91a19e 100644
--- a/arithmetic_analysis/jacobi_iteration_method.py
+++ b/linear_algebra/jacobi_iteration_method.py
@@ -1,203 +1,203 @@
-"""
-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: NDArray[float64],
-    constant_matrix: NDArray[float64],
-    init_val: list[float],
-    iterations: int,
-) -> list[float]:
-    """
-    Jacobi Iteration Method:
-    An iterative algorithm to determine the solutions of strictly diagonally dominant
-    system of linear equations
-
-    4x1 +  x2 +  x3 =  2
-     x1 + 5x2 + 2x3 = -6
-     x1 + 2x2 + 4x3 = -4
-
-    x_init = [0.5, -0.5 , -0.5]
-
-    Examples:
-
-    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
-    >>> constant = np.array([[2], [-6], [-4]])
-    >>> init_val = [0.5, -0.5, -0.5]
-    >>> iterations = 3
-    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
-    [0.909375, -1.14375, -0.7484375]
-
-
-    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2]])
-    >>> constant = np.array([[2], [-6], [-4]])
-    >>> init_val = [0.5, -0.5, -0.5]
-    >>> iterations = 3
-    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
-    Traceback (most recent call last):
-        ...
-    ValueError: Coefficient matrix dimensions must be nxn but received 2x3
-
-    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
-    >>> constant = np.array([[2], [-6]])
-    >>> init_val = [0.5, -0.5, -0.5]
-    >>> iterations = 3
-    >>> jacobi_iteration_method(
-    ...     coefficient, constant, init_val, iterations
-    ... )  # doctest: +NORMALIZE_WHITESPACE
-    Traceback (most recent call last):
-        ...
-    ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but
-                received 3x3 and 2x1
-
-    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
-    >>> constant = np.array([[2], [-6], [-4]])
-    >>> init_val = [0.5, -0.5]
-    >>> iterations = 3
-    >>> jacobi_iteration_method(
-    ...     coefficient, constant, init_val, iterations
-    ... )  # doctest: +NORMALIZE_WHITESPACE
-    Traceback (most recent call last):
-        ...
-    ValueError: Number of initial values must be equal to number of rows in coefficient
-                matrix but received 2 and 3
-
-    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
-    >>> constant = np.array([[2], [-6], [-4]])
-    >>> init_val = [0.5, -0.5, -0.5]
-    >>> iterations = 0
-    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
-    Traceback (most recent call last):
-        ...
-    ValueError: Iterations must be at least 1
-    """
-
-    rows1, cols1 = coefficient_matrix.shape
-    rows2, cols2 = constant_matrix.shape
-
-    if rows1 != cols1:
-        msg = f"Coefficient matrix dimensions must be nxn but received {rows1}x{cols1}"
-        raise ValueError(msg)
-
-    if cols2 != 1:
-        msg = f"Constant matrix must be nx1 but received {rows2}x{cols2}"
-        raise ValueError(msg)
-
-    if rows1 != rows2:
-        msg = (
-            "Coefficient and constant matrices dimensions must be nxn and nx1 but "
-            f"received {rows1}x{cols1} and {rows2}x{cols2}"
-        )
-        raise ValueError(msg)
-
-    if len(init_val) != rows1:
-        msg = (
-            "Number of initial values must be equal to number of rows in coefficient "
-            f"matrix but received {len(init_val)} and {rows1}"
-        )
-        raise ValueError(msg)
-
-    if iterations <= 0:
-        raise ValueError("Iterations must be at least 1")
-
-    table: NDArray[float64] = np.concatenate(
-        (coefficient_matrix, constant_matrix), axis=1
-    )
-
-    rows, cols = table.shape
-
-    strictly_diagonally_dominant(table)
-
-    """
-    # Iterates the whole matrix for given number of times
-    for _ in range(iterations):
-        new_val = []
-        for row in range(rows):
-            temp = 0
-            for col in range(cols):
-                if col == row:
-                    denom = table[row][col]
-                elif col == cols - 1:
-                    val = table[row][col]
-                else:
-                    temp += (-1) * table[row][col] * init_val[col]
-            temp = (temp + val) / denom
-            new_val.append(temp)
-        init_val = new_val
-    """
-
-    # denominator - a list of values along the diagonal
-    denominator = np.diag(coefficient_matrix)
-
-    # val_last - values of the last column of the table array
-    val_last = table[:, -1]
-
-    # masks - boolean mask of all strings without diagonal
-    # elements array coefficient_matrix
-    masks = ~np.eye(coefficient_matrix.shape[0], dtype=bool)
-
-    # no_diagonals - coefficient_matrix array values without diagonal elements
-    no_diagonals = coefficient_matrix[masks].reshape(-1, rows - 1)
-
-    # Here we get 'i_col' - these are the column numbers, for each row
-    # without diagonal elements, except for the last column.
-    i_row, i_col = np.where(masks)
-    ind = i_col.reshape(-1, rows - 1)
-
-    #'i_col' is converted to a two-dimensional list 'ind', which will be
-    # used to make selections from 'init_val' ('arr' array see below).
-
-    # Iterates the whole matrix for given number of times
-    for _ in range(iterations):
-        arr = np.take(init_val, ind)
-        sum_product_rows = np.sum((-1) * no_diagonals * arr, axis=1)
-        new_val = (sum_product_rows + val_last) / denominator
-        init_val = new_val
-
-    return new_val.tolist()
-
-
-# Checks if the given matrix is strictly diagonally dominant
-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)
-    True
-
-    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 3, -4]])
-    >>> strictly_diagonally_dominant(table)
-    Traceback (most recent call last):
-        ...
-    ValueError: Coefficient matrix is not strictly diagonally dominant
-    """
-
-    rows, cols = table.shape
-
-    is_diagonally_dominant = True
-
-    for i in range(rows):
-        total = 0
-        for j in range(cols - 1):
-            if i == j:
-                continue
-            else:
-                total += table[i][j]
-
-        if table[i][i] <= total:
-            raise ValueError("Coefficient matrix is not strictly diagonally dominant")
-
-    return is_diagonally_dominant
-
-
-# Test Cases
-if __name__ == "__main__":
-    import doctest
-
-    doctest.testmod()
+"""
+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: NDArray[float64],
+    constant_matrix: NDArray[float64],
+    init_val: list[float],
+    iterations: int,
+) -> list[float]:
+    """
+    Jacobi Iteration Method:
+    An iterative algorithm to determine the solutions of strictly diagonally dominant
+    system of linear equations
+
+    4x1 +  x2 +  x3 =  2
+     x1 + 5x2 + 2x3 = -6
+     x1 + 2x2 + 4x3 = -4
+
+    x_init = [0.5, -0.5 , -0.5]
+
+    Examples:
+
+    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
+    >>> constant = np.array([[2], [-6], [-4]])
+    >>> init_val = [0.5, -0.5, -0.5]
+    >>> iterations = 3
+    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
+    [0.909375, -1.14375, -0.7484375]
+
+
+    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2]])
+    >>> constant = np.array([[2], [-6], [-4]])
+    >>> init_val = [0.5, -0.5, -0.5]
+    >>> iterations = 3
+    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
+    Traceback (most recent call last):
+        ...
+    ValueError: Coefficient matrix dimensions must be nxn but received 2x3
+
+    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
+    >>> constant = np.array([[2], [-6]])
+    >>> init_val = [0.5, -0.5, -0.5]
+    >>> iterations = 3
+    >>> jacobi_iteration_method(
+    ...     coefficient, constant, init_val, iterations
+    ... )  # doctest: +NORMALIZE_WHITESPACE
+    Traceback (most recent call last):
+        ...
+    ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but
+                received 3x3 and 2x1
+
+    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
+    >>> constant = np.array([[2], [-6], [-4]])
+    >>> init_val = [0.5, -0.5]
+    >>> iterations = 3
+    >>> jacobi_iteration_method(
+    ...     coefficient, constant, init_val, iterations
+    ... )  # doctest: +NORMALIZE_WHITESPACE
+    Traceback (most recent call last):
+        ...
+    ValueError: Number of initial values must be equal to number of rows in coefficient
+                matrix but received 2 and 3
+
+    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
+    >>> constant = np.array([[2], [-6], [-4]])
+    >>> init_val = [0.5, -0.5, -0.5]
+    >>> iterations = 0
+    >>> jacobi_iteration_method(coefficient, constant, init_val, iterations)
+    Traceback (most recent call last):
+        ...
+    ValueError: Iterations must be at least 1
+    """
+
+    rows1, cols1 = coefficient_matrix.shape
+    rows2, cols2 = constant_matrix.shape
+
+    if rows1 != cols1:
+        msg = f"Coefficient matrix dimensions must be nxn but received {rows1}x{cols1}"
+        raise ValueError(msg)
+
+    if cols2 != 1:
+        msg = f"Constant matrix must be nx1 but received {rows2}x{cols2}"
+        raise ValueError(msg)
+
+    if rows1 != rows2:
+        msg = (
+            "Coefficient and constant matrices dimensions must be nxn and nx1 but "
+            f"received {rows1}x{cols1} and {rows2}x{cols2}"
+        )
+        raise ValueError(msg)
+
+    if len(init_val) != rows1:
+        msg = (
+            "Number of initial values must be equal to number of rows in coefficient "
+            f"matrix but received {len(init_val)} and {rows1}"
+        )
+        raise ValueError(msg)
+
+    if iterations <= 0:
+        raise ValueError("Iterations must be at least 1")
+
+    table: NDArray[float64] = np.concatenate(
+        (coefficient_matrix, constant_matrix), axis=1
+    )
+
+    rows, cols = table.shape
+
+    strictly_diagonally_dominant(table)
+
+    """
+    # Iterates the whole matrix for given number of times
+    for _ in range(iterations):
+        new_val = []
+        for row in range(rows):
+            temp = 0
+            for col in range(cols):
+                if col == row:
+                    denom = table[row][col]
+                elif col == cols - 1:
+                    val = table[row][col]
+                else:
+                    temp += (-1) * table[row][col] * init_val[col]
+            temp = (temp + val) / denom
+            new_val.append(temp)
+        init_val = new_val
+    """
+
+    # denominator - a list of values along the diagonal
+    denominator = np.diag(coefficient_matrix)
+
+    # val_last - values of the last column of the table array
+    val_last = table[:, -1]
+
+    # masks - boolean mask of all strings without diagonal
+    # elements array coefficient_matrix
+    masks = ~np.eye(coefficient_matrix.shape[0], dtype=bool)
+
+    # no_diagonals - coefficient_matrix array values without diagonal elements
+    no_diagonals = coefficient_matrix[masks].reshape(-1, rows - 1)
+
+    # Here we get 'i_col' - these are the column numbers, for each row
+    # without diagonal elements, except for the last column.
+    i_row, i_col = np.where(masks)
+    ind = i_col.reshape(-1, rows - 1)
+
+    #'i_col' is converted to a two-dimensional list 'ind', which will be
+    # used to make selections from 'init_val' ('arr' array see below).
+
+    # Iterates the whole matrix for given number of times
+    for _ in range(iterations):
+        arr = np.take(init_val, ind)
+        sum_product_rows = np.sum((-1) * no_diagonals * arr, axis=1)
+        new_val = (sum_product_rows + val_last) / denominator
+        init_val = new_val
+
+    return new_val.tolist()
+
+
+# Checks if the given matrix is strictly diagonally dominant
+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)
+    True
+
+    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 3, -4]])
+    >>> strictly_diagonally_dominant(table)
+    Traceback (most recent call last):
+        ...
+    ValueError: Coefficient matrix is not strictly diagonally dominant
+    """
+
+    rows, cols = table.shape
+
+    is_diagonally_dominant = True
+
+    for i in range(rows):
+        total = 0
+        for j in range(cols - 1):
+            if i == j:
+                continue
+            else:
+                total += table[i][j]
+
+        if table[i][i] <= total:
+            raise ValueError("Coefficient matrix is not strictly diagonally dominant")
+
+    return is_diagonally_dominant
+
+
+# Test Cases
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/arithmetic_analysis/lu_decomposition.py b/linear_algebra/lu_decomposition.py
similarity index 100%
rename from arithmetic_analysis/lu_decomposition.py
rename to linear_algebra/lu_decomposition.py
diff --git a/arithmetic_analysis/bisection.py b/maths/numerical_analysis/bisection.py
similarity index 100%
rename from arithmetic_analysis/bisection.py
rename to maths/numerical_analysis/bisection.py
diff --git a/maths/bisection.py b/maths/numerical_analysis/bisection_2.py
similarity index 100%
rename from maths/bisection.py
rename to maths/numerical_analysis/bisection_2.py
diff --git a/maths/integration_by_simpson_approx.py b/maths/numerical_analysis/integration_by_simpson_approx.py
similarity index 100%
rename from maths/integration_by_simpson_approx.py
rename to maths/numerical_analysis/integration_by_simpson_approx.py
diff --git a/arithmetic_analysis/intersection.py b/maths/numerical_analysis/intersection.py
similarity index 100%
rename from arithmetic_analysis/intersection.py
rename to maths/numerical_analysis/intersection.py
diff --git a/maths/nevilles_method.py b/maths/numerical_analysis/nevilles_method.py
similarity index 100%
rename from maths/nevilles_method.py
rename to maths/numerical_analysis/nevilles_method.py
diff --git a/arithmetic_analysis/newton_forward_interpolation.py b/maths/numerical_analysis/newton_forward_interpolation.py
similarity index 100%
rename from arithmetic_analysis/newton_forward_interpolation.py
rename to maths/numerical_analysis/newton_forward_interpolation.py
diff --git a/arithmetic_analysis/newton_method.py b/maths/numerical_analysis/newton_method.py
similarity index 100%
rename from arithmetic_analysis/newton_method.py
rename to maths/numerical_analysis/newton_method.py
diff --git a/arithmetic_analysis/newton_raphson.py b/maths/numerical_analysis/newton_raphson.py
similarity index 61%
rename from arithmetic_analysis/newton_raphson.py
rename to maths/numerical_analysis/newton_raphson.py
index 1b90ad41..8491ca80 100644
--- a/arithmetic_analysis/newton_raphson.py
+++ b/maths/numerical_analysis/newton_raphson.py
@@ -5,42 +5,41 @@
 from __future__ import annotations
 
 from decimal import Decimal
-from math import *  # noqa: F403
 
-from sympy import diff
+from sympy import diff, lambdify, symbols
 
 
-def newton_raphson(
-    func: str, a: float | Decimal, precision: float = 10**-10
-) -> float:
+def newton_raphson(func: str, a: float | Decimal, precision: float = 1e-10) -> float:
     """Finds root from the point 'a' onwards by Newton-Raphson method
     >>> newton_raphson("sin(x)", 2)
     3.1415926536808043
-    >>> newton_raphson("x**2 - 5*x +2", 0.4)
+    >>> newton_raphson("x**2 - 5*x + 2", 0.4)
     0.4384471871911695
     >>> newton_raphson("x**2 - 5", 0.1)
     2.23606797749979
-    >>> newton_raphson("log(x)- 1", 2)
+    >>> newton_raphson("log(x) - 1", 2)
     2.718281828458938
     """
-    x = a
+    x = symbols("x")
+    f = lambdify(x, func, "math")
+    f_derivative = lambdify(x, diff(func), "math")
+    x_curr = a
     while True:
-        x = Decimal(x) - (
-            Decimal(eval(func)) / Decimal(eval(str(diff(func))))  # noqa: S307
-        )
-        # This number dictates the accuracy of the answer
-        if abs(eval(func)) < precision:  # noqa: S307
-            return float(x)
+        x_curr = Decimal(x_curr) - Decimal(f(x_curr)) / Decimal(f_derivative(x_curr))
+        if abs(f(x_curr)) < precision:
+            return float(x_curr)
 
 
-# Let's Execute
 if __name__ == "__main__":
-    # Find root of trigonometric function
+    import doctest
+
+    doctest.testmod()
+
     # Find value of pi
     print(f"The root of sin(x) = 0 is {newton_raphson('sin(x)', 2)}")
     # Find root of polynomial
     print(f"The root of x**2 - 5*x + 2 = 0 is {newton_raphson('x**2 - 5*x + 2', 0.4)}")
-    # Find Square Root of 5
+    # Find value of e
     print(f"The root of log(x) - 1 = 0 is {newton_raphson('log(x) - 1', 2)}")
-    # Exponential Roots
+    # Find root of exponential function
     print(f"The root of exp(x) - 1 = 0 is {newton_raphson('exp(x) - 1', 0)}")
diff --git a/maths/newton_raphson.py b/maths/numerical_analysis/newton_raphson_2.py
similarity index 100%
rename from maths/newton_raphson.py
rename to maths/numerical_analysis/newton_raphson_2.py
diff --git a/arithmetic_analysis/newton_raphson_new.py b/maths/numerical_analysis/newton_raphson_new.py
similarity index 100%
rename from arithmetic_analysis/newton_raphson_new.py
rename to maths/numerical_analysis/newton_raphson_new.py
diff --git a/maths/numerical_integration.py b/maths/numerical_analysis/numerical_integration.py
similarity index 100%
rename from maths/numerical_integration.py
rename to maths/numerical_analysis/numerical_integration.py
diff --git a/maths/runge_kutta.py b/maths/numerical_analysis/runge_kutta.py
similarity index 100%
rename from maths/runge_kutta.py
rename to maths/numerical_analysis/runge_kutta.py
diff --git a/maths/runge_kutta_fehlberg_45.py b/maths/numerical_analysis/runge_kutta_fehlberg_45.py
similarity index 100%
rename from maths/runge_kutta_fehlberg_45.py
rename to maths/numerical_analysis/runge_kutta_fehlberg_45.py
diff --git a/arithmetic_analysis/secant_method.py b/maths/numerical_analysis/secant_method.py
similarity index 100%
rename from arithmetic_analysis/secant_method.py
rename to maths/numerical_analysis/secant_method.py
diff --git a/maths/simpson_rule.py b/maths/numerical_analysis/simpson_rule.py
similarity index 100%
rename from maths/simpson_rule.py
rename to maths/numerical_analysis/simpson_rule.py
diff --git a/maths/square_root.py b/maths/numerical_analysis/square_root.py
similarity index 100%
rename from maths/square_root.py
rename to maths/numerical_analysis/square_root.py
diff --git a/arithmetic_analysis/image_data/2D_problems.jpg b/physics/image_data/2D_problems.jpg
similarity index 100%
rename from arithmetic_analysis/image_data/2D_problems.jpg
rename to physics/image_data/2D_problems.jpg
diff --git a/arithmetic_analysis/image_data/2D_problems_1.jpg b/physics/image_data/2D_problems_1.jpg
similarity index 100%
rename from arithmetic_analysis/image_data/2D_problems_1.jpg
rename to physics/image_data/2D_problems_1.jpg
diff --git a/arithmetic_analysis/__init__.py b/physics/image_data/__init__.py
similarity index 100%
rename from arithmetic_analysis/__init__.py
rename to physics/image_data/__init__.py
diff --git a/arithmetic_analysis/in_static_equilibrium.py b/physics/in_static_equilibrium.py
similarity index 96%
rename from arithmetic_analysis/in_static_equilibrium.py
rename to physics/in_static_equilibrium.py
index 7aaecf17..d56299f6 100644
--- a/arithmetic_analysis/in_static_equilibrium.py
+++ b/physics/in_static_equilibrium.py
@@ -1,94 +1,94 @@
-"""
-Checks if a system of forces is in static equilibrium.
-"""
-from __future__ import annotations
-
-from numpy import array, cos, cross, float64, radians, sin
-from numpy.typing import NDArray
-
-
-def polar_force(
-    magnitude: float, angle: float, radian_mode: bool = False
-) -> list[float]:
-    """
-    Resolves force along rectangular components.
-    (force, angle) => (force_x, force_y)
-    >>> import math
-    >>> force = polar_force(10, 45)
-    >>> math.isclose(force[0], 7.071067811865477)
-    True
-    >>> math.isclose(force[1], 7.0710678118654755)
-    True
-    >>> force = polar_force(10, 3.14, radian_mode=True)
-    >>> math.isclose(force[0], -9.999987317275396)
-    True
-    >>> math.isclose(force[1], 0.01592652916486828)
-    True
-    """
-    if radian_mode:
-        return [magnitude * cos(angle), magnitude * sin(angle)]
-    return [magnitude * cos(radians(angle)), magnitude * sin(radians(angle))]
-
-
-def in_static_equilibrium(
-    forces: NDArray[float64], location: NDArray[float64], eps: float = 10**-1
-) -> bool:
-    """
-    Check if a system is in equilibrium.
-    It takes two numpy.array objects.
-    forces ==>  [
-                        [force1_x, force1_y],
-                        [force2_x, force2_y],
-                        ....]
-    location ==>  [
-                        [x1, y1],
-                        [x2, y2],
-                        ....]
-    >>> force = array([[1, 1], [-1, 2]])
-    >>> location = array([[1, 0], [10, 0]])
-    >>> in_static_equilibrium(force, location)
-    False
-    """
-    # summation of moments is zero
-    moments: NDArray[float64] = cross(location, forces)
-    sum_moments: float = sum(moments)
-    return abs(sum_moments) < eps
-
-
-if __name__ == "__main__":
-    # Test to check if it works
-    forces = array(
-        [
-            polar_force(718.4, 180 - 30),
-            polar_force(879.54, 45),
-            polar_force(100, -90),
-        ]
-    )
-
-    location: NDArray[float64] = array([[0, 0], [0, 0], [0, 0]])
-
-    assert in_static_equilibrium(forces, location)
-
-    # Problem 1 in image_data/2D_problems.jpg
-    forces = array(
-        [
-            polar_force(30 * 9.81, 15),
-            polar_force(215, 180 - 45),
-            polar_force(264, 90 - 30),
-        ]
-    )
-
-    location = array([[0, 0], [0, 0], [0, 0]])
-
-    assert in_static_equilibrium(forces, location)
-
-    # Problem in image_data/2D_problems_1.jpg
-    forces = array([[0, -2000], [0, -1200], [0, 15600], [0, -12400]])
-
-    location = array([[0, 0], [6, 0], [10, 0], [12, 0]])
-
-    assert in_static_equilibrium(forces, location)
-
-    import doctest
-
-    doctest.testmod()
+"""
+Checks if a system of forces is in static equilibrium.
+"""
+from __future__ import annotations
+
+from numpy import array, cos, cross, float64, radians, sin
+from numpy.typing import NDArray
+
+
+def polar_force(
+    magnitude: float, angle: float, radian_mode: bool = False
+) -> list[float]:
+    """
+    Resolves force along rectangular components.
+    (force, angle) => (force_x, force_y)
+    >>> import math
+    >>> force = polar_force(10, 45)
+    >>> math.isclose(force[0], 7.071067811865477)
+    True
+    >>> math.isclose(force[1], 7.0710678118654755)
+    True
+    >>> force = polar_force(10, 3.14, radian_mode=True)
+    >>> math.isclose(force[0], -9.999987317275396)
+    True
+    >>> math.isclose(force[1], 0.01592652916486828)
+    True
+    """
+    if radian_mode:
+        return [magnitude * cos(angle), magnitude * sin(angle)]
+    return [magnitude * cos(radians(angle)), magnitude * sin(radians(angle))]
+
+
+def in_static_equilibrium(
+    forces: NDArray[float64], location: NDArray[float64], eps: float = 10**-1
+) -> bool:
+    """
+    Check if a system is in equilibrium.
+    It takes two numpy.array objects.
+    forces ==>  [
+                        [force1_x, force1_y],
+                        [force2_x, force2_y],
+                        ....]
+    location ==>  [
+                        [x1, y1],
+                        [x2, y2],
+                        ....]
+    >>> force = array([[1, 1], [-1, 2]])
+    >>> location = array([[1, 0], [10, 0]])
+    >>> in_static_equilibrium(force, location)
+    False
+    """
+    # summation of moments is zero
+    moments: NDArray[float64] = cross(location, forces)
+    sum_moments: float = sum(moments)
+    return abs(sum_moments) < eps
+
+
+if __name__ == "__main__":
+    # Test to check if it works
+    forces = array(
+        [
+            polar_force(718.4, 180 - 30),
+            polar_force(879.54, 45),
+            polar_force(100, -90),
+        ]
+    )
+
+    location: NDArray[float64] = array([[0, 0], [0, 0], [0, 0]])
+
+    assert in_static_equilibrium(forces, location)
+
+    # Problem 1 in image_data/2D_problems.jpg
+    forces = array(
+        [
+            polar_force(30 * 9.81, 15),
+            polar_force(215, 180 - 45),
+            polar_force(264, 90 - 30),
+        ]
+    )
+
+    location = array([[0, 0], [0, 0], [0, 0]])
+
+    assert in_static_equilibrium(forces, location)
+
+    # Problem in image_data/2D_problems_1.jpg
+    forces = array([[0, -2000], [0, -1200], [0, 15600], [0, -12400]])
+
+    location = array([[0, 0], [6, 0], [10, 0], [12, 0]])
+
+    assert in_static_equilibrium(forces, location)
+
+    import doctest
+
+    doctest.testmod()