From 00da7f87251c2a78701d8342846fc753f357c11e Mon Sep 17 00:00:00 2001 From: Lena Castillo Date: Mon, 13 Jan 2025 15:19:43 -0500 Subject: [PATCH] Nuevas funciones --- marlonpy/marlonpy/Interpolation/Lagrange.py | 122 +++++++++++++++ marlonpy/marlonpy/SystemEquations/Cholesky.py | 145 ++++++++++++++++++ marlonpy/marlonpy/SystemEquations/__init__.py | 3 +- marlonpy/marlonpy/__init__.py | 3 + marlonpy/marlonpy/requirements.txt | 4 +- 5 files changed, 275 insertions(+), 2 deletions(-) create mode 100644 marlonpy/marlonpy/Interpolation/Lagrange.py create mode 100644 marlonpy/marlonpy/SystemEquations/Cholesky.py diff --git a/marlonpy/marlonpy/Interpolation/Lagrange.py b/marlonpy/marlonpy/Interpolation/Lagrange.py new file mode 100644 index 0000000..fc5d7ed --- /dev/null +++ b/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() +""" diff --git a/marlonpy/marlonpy/SystemEquations/Cholesky.py b/marlonpy/marlonpy/SystemEquations/Cholesky.py new file mode 100644 index 0000000..4295209 --- /dev/null +++ b/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() +""" \ No newline at end of file diff --git a/marlonpy/marlonpy/SystemEquations/__init__.py b/marlonpy/marlonpy/SystemEquations/__init__.py index 88ecd43..3015fd0 100644 --- a/marlonpy/marlonpy/SystemEquations/__init__.py +++ b/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 \ No newline at end of file +#from Jacobi import jacobi +#from Cholesky import Cholesky \ No newline at end of file diff --git a/marlonpy/marlonpy/__init__.py b/marlonpy/marlonpy/__init__.py index 3d20afe..1218de3 100644 --- a/marlonpy/marlonpy/__init__.py +++ b/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 + diff --git a/marlonpy/marlonpy/requirements.txt b/marlonpy/marlonpy/requirements.txt index 820f99c..b6f326d 100644 --- a/marlonpy/marlonpy/requirements.txt +++ b/marlonpy/marlonpy/requirements.txt @@ -1,3 +1,5 @@ sympy numpy -tabulate \ No newline at end of file +tabulate +fractions +ipython \ No newline at end of file