Nuevas funciones

marlonpy/marlonpy/Interpolation/Lagrange.py

@@ -0,0 +1,122 @@
+import sympy as sp
+from sympy import latex
+from fractions import Fraction
+
+"""
+Performs Lagrange polynomial interpolation for a given set of data points.
+
+Args:
+    x_values (list or numpy.ndarray): Known x data points.
+    y_values (list or numpy.ndarray): Corresponding y data points.
+
+Attributes:
+    x_values (list or numpy.ndarray): Known x data points.
+    y_values (list or numpy.ndarray): Corresponding y data points.
+    n (int): Number of data points.
+    poly (sympy.Expr): Simplified Lagrange interpolating polynomial.
+
+Methods:
+    build_polynomial():
+        Constructs the Lagrange interpolating polynomial using the given x and y data points.
+        Simplifies the polynomial and displays it in LaTeX format.
+
+    evaluate_roots():
+        Evaluates the Lagrange polynomial at each x data point.
+        Displays the calculation steps, including the contribution of each term, 
+        and confirms if the evaluated values match the corresponding y data points.
+
+    evaluate_at_point(x_sub):
+        Evaluates the Lagrange polynomial at a specific point x_sub.
+        Args:
+            x_sub (float): The x-value where the polynomial is evaluated.
+        Returns:
+            str: A formatted string with the evaluated value in both decimal 
+                 and fractional forms. Example: "f(2) = 3.3333 = 10/3".
+
+    interpolate():
+        Performs the complete interpolation process:
+        - Constructs the Lagrange polynomial.
+        - Evaluates the polynomial at the given x data points.
+        - Confirms the interpolation accuracy.
+
+Returns:
+    str: A formatted string representing the evaluation process and results,
+         showing the polynomial, intermediate calculations, and final results
+         with both decimal values and their fractional equivalents.
+"""
+
+class Lagrange:
+    def __init__(self, x_values, y_values):
+        self.x_values = x_values  # Known x data points
+        self.y_values = y_values  # Corresponding y data points
+        self.n = len(x_values)    # Number of data points
+        self.poly = None          # Lagrange polynomial
+
+    def build_polynomial(self):
+        x = sp.Symbol("x")
+        polynomial_sum = 0
+
+        # Construct the Lagrange interpolating polynomial
+        for i in range(self.n):  # Summation
+            term = 1
+            for j in range(self.n):  # Product
+                if i != j:
+                    term *= (x - self.x_values[j]) / (self.x_values[i] - self.x_values[j])
+            polynomial_term = self.y_values[i] * term
+            polynomial_sum += polynomial_term
+
+        self.poly = sp.simplify(polynomial_sum)
+        print("\nInterpolating Polynomial:")
+        print(f"f(x) = {latex(self.poly)}")  # Display the polynomial in LaTeX format
+
+    def evaluate_roots(self):
+        print("\nNow, we evaluate the polynomial at the given roots:")
+        x = sp.Symbol("x")
+        for xi in self.x_values:
+            terms = []  # List to store the terms in string format
+            term_values = []  # List to store the evaluated values of each term
+            result = 0  # Initialize result to 0
+
+            # Extract each term of the polynomial and replace the 'x' with each value of x_data
+            for coef in self.poly.as_ordered_terms():
+                if coef.has(x):
+                    # Format the term as 'a*(x)^n'
+                    term_str = str(coef).replace("x", f"({xi})")  # Replace x with xi
+                    term_value = coef.subs(x, xi)  # Calculate the value for each term
+                    terms.append(term_str)
+                    term_values.append(term_value)
+                else:
+                    terms.append(str(coef))
+                    term_values.append(coef)
+
+            # Join all terms into a string representation
+            expression = " + ".join(terms).replace("+ -", "- ")
+            simplified_terms = " + ".join(
+                [f"{term_value:.2f}" if term_value >= 0 else f"- {abs(term_value):.2f}" for term_value in term_values]
+            ).replace("+ -", "- ")
+
+            # Final result after evaluating each term
+            result = sum(term_values)
+            print(f"f({xi}) = {expression} = {simplified_terms} = {result:.2f}")
+
+            # Confirm if the result matches the corresponding y_value
+            if result == self.y_values[self.x_values.index(xi)]:
+                print(f"Thus, f({xi}) = {result:.2f} matches y{self.x_values.index(xi)}")
+
+    def evaluate_at_point(self, x_sub):
+        y_sub = self.poly.subs(sp.Symbol('x'), x_sub)
+        fractional_result = Fraction(float(y_sub)).limit_denominator()
+        print(f"f({x_sub}) = {y_sub:.2f} = {fractional_result}")
+
+    def interpolate(self):
+        self.build_polynomial()
+        self.evaluate_roots()
+
+"""
+# Example usage
+x_data = [1, 2, 3]
+y_data = [3, 1, 5]
+
+lagrange = Lagrange(x_data, y_data)
+lagrange.interpolate()
+"""

marlonpy/marlonpy/SystemEquations/Cholesky.py

@@ -0,0 +1,145 @@
+import numpy as np
+from fractions import Fraction
+
+class Cholesky:
+    """
+Performs Cholesky decomposition of a symmetric, positive definite matrix A and solves a linear system Ax = sol.
+
+Args:
+    A (numpy.ndarray): Coefficient matrix of the linear system.
+    sol (numpy.ndarray): Vector of constants in the linear system.
+
+Attributes:
+    L (numpy.ndarray): Lower triangular matrix obtained from Cholesky decomposition.
+    x (numpy.ndarray): Solution vector containing the values of the variables that satisfy Ax = sol.
+    y (numpy.ndarray): Intermediate solution vector for the system Ly = sol.
+
+Returns:
+    str: A formatted string representing the solution vector x with both decimal values and their fractional equivalents.
+         Example: "Thus, x = 28.5833 = 343/12, y = -7.6667 = -23/3, z = 1.3333 = 4/3".
+"""
+
+    def __init__(self, A, sol):
+        self.A = A  # Coefficient matrix
+        self.sol = sol  # Solution vector
+        self.n = A.shape[0]  # Matrix size
+        self.L = np.zeros((self.n, self.n))  # Lower triangular matrix
+        self.y = np.zeros(self.n)  # Intermediate vector y
+        self.x = np.zeros(self.n)  # Solution vector x
+
+    def is_symmetric(self):
+        return np.allclose(self.A, self.transpose(self.A))
+
+    def transpose(self, matrix):
+        return matrix.T
+
+    def det(self, submatrix):
+        submatrix = np.array(submatrix)  # Ensure it's a NumPy array
+        if submatrix.shape == (1, 1):  # Handle single-element matrices
+            D = submatrix[0, 0]
+        elif submatrix.shape == (2, 2):  # Handle 2x2 matrices
+            D = submatrix[0, 0] * submatrix[1, 1] - submatrix[0, 1] * submatrix[1, 0]
+        elif submatrix.shape == (3, 3):  # Handle 3x3 matrices
+            D = (
+                submatrix[0, 0] * submatrix[1, 1] * submatrix[2, 2]
+                + submatrix[1, 0] * submatrix[2, 1] * submatrix[0, 2]
+                + submatrix[2, 0] * submatrix[0, 1] * submatrix[1, 2]
+                - submatrix[0, 2] * submatrix[1, 1] * submatrix[2, 0]
+                - submatrix[1, 2] * submatrix[2, 1] * submatrix[0, 0]
+                - submatrix[2, 2] * submatrix[0, 1] * submatrix[1, 0]
+            )
+        else:  # For larger matrices
+            D = 0
+            size = len(submatrix)
+            for i in range(size):
+                for j in range(size):
+                    sign = (-1) ** j
+                    sub = [
+                        [submatrix[row][col] for col in range(size) if col != j]
+                        for row in range(size) if row != i
+                    ]
+                    sub = np.array(sub)  # Ensure submatrix is a NumPy array
+                    D += sign * submatrix[0][j] * self.det(sub)
+        return D
+
+    def is_positive_definite(self):
+        for k in range(1, self.n + 1):
+            submatrix = self.A[:k, :k]
+            D = self.det(submatrix)
+            if D <= 0:
+                return False
+
+        return True
+
+    def decompose(self):
+        print("\033[1mMatrix A:\033[0m")
+        print(self.A)
+        print()
+        if self.is_symmetric():
+            if self.is_positive_definite():
+                print("\033[1mThe matrix is symmetric and positive definite.\033[0m\n")
+            else:
+                print("\033[1mThe matrix is symmetric but not positive definite.\033[0m\n")
+        else:
+            if self.is_positive_definite():
+                print("\033[1mThe matrix is not symmetric but it is positive definite.\033[0m\n")
+            else:
+                print("\033[1mThe matrix is neither symmetric nor positive definite.\033[0m\n")
+
+        if self.is_symmetric():
+            if self.is_positive_definite():
+                # Step 2: Calculate L
+                for i in range(self.n):
+                    for j in range(i + 1):
+                        sum_val = self.A[i, j]
+                        for k in range(j):
+                            sum_val -= self.L[i, k] * self.L[j, k]
+                        if i == j:
+                            self.L[i, j] = np.sqrt(sum_val)
+                        else:
+                            self.L[i, j] = sum_val / self.L[j, j]
+                    print(f"\033[1mIteration {i + 1}:\033[0m")
+                    print("Matrix L (Lower Triangular):")
+                    print(self.L)
+                    print()
+
+                # Step 2.1: Solve Ly = b
+                for i in range(self.n):
+                    self.y[i] = (self.sol[i] - np.dot(self.L[i, :i], self.y[:i])) / self.L[i, i]
+                print("\033[1my = \033[0m", self.y, "\n")
+
+                # Step 3: Calculate U
+                U = self.transpose(self.L)
+
+                # Print upper triangular matrix U
+                print("Matrix U (Upper Triangular):")
+                print(U, "\n")
+
+                # Step 3.1: Solve Ux = y
+                for i in range(self.n - 1, -1, -1):
+                    self.x[i] = (self.y[i] - np.dot(self.L[i + 1:self.n, i], self.x[i + 1:self.n])) / self.L[i, i]
+
+                # Show the solutions x
+                print("\033[1mx = \033[0m", self.x)
+
+                self.show_results()
+
+    def show_results(self):
+        sol_x = float(self.x[0].item())
+        sol_y = float(self.x[1].item())
+        sol_z = float(self.x[2].item())
+
+        # Convert solutions to fractions
+        frac_x = Fraction(sol_x).limit_denominator()
+        frac_y = Fraction(sol_y).limit_denominator()
+        frac_z = Fraction(sol_z).limit_denominator()
+
+        print(f"\n\033[1mThus, x = {self.x[0]:.4f} = {frac_x}, y = {self.x[1]:.4f} = {frac_y}, z = {self.x[2]:.4f} = {frac_z}\033[0m")
+
+""""
+A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]])
+sol = np.array([1, 2, 3, 4])
+
+cholesky = Cholesky(A, sol)
+cholesky.decompose()
+"""

marlonpy/marlonpy/SystemEquations/__init__.py

@@ -2,4 +2,5 @@
 #from DDM import DDM
 #from Doolittle import doolittle
 #from GaussSeidel import gauss_seidel
-#from Jacobi import jacobi
+#from Jacobi import jacobi
+#from Cholesky import Cholesky

marlonpy/marlonpy/__init__.py

@@ -17,5 +17,8 @@
 from .NumericalIntegration.TrapezoidalRule import trapezoidal
 from .NumericalIntegration.GaussLegendre import gauss_legendre
 from .DifferentialEquations.RungeKutta4th import runge_kutta
+from .SystemEquations.Cholesky import Cholesky
+from .Interpolation import Lagrange
+
 
 

marlonpy/marlonpy/requirements.txt

@@ -1,3 +1,5 @@
 sympy
 numpy
-tabulate
+tabulate
+fractions 
+ipython