diff --git a/.gitignore b/.gitignore index b282c0c..c28a45b 100644 --- a/.gitignore +++ b/.gitignore @@ -161,3 +161,4 @@ cython_debug/ .idea/ /src/error_config.py /main.py +/double_integrators.py diff --git a/.pylintrc b/.pylintrc index 8e0a6db..1188020 100644 --- a/.pylintrc +++ b/.pylintrc @@ -183,12 +183,13 @@ good-names=i, j, k, ex, + dx, Run, _ # Good variable names regexes, separated by a comma. If names match any regex, # they will always be accepted -good-names-rgxs= +good-names-rgxs=(f|k|x|y|z)[0-9]?(_?)* # Include a hint for the correct naming format with invalid-name. include-naming-hint=no diff --git a/README.md b/README.md index e5487c7..452a21f 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Solvers +# num_solvers A library for some basic numerical solvers. diff --git a/setup.cfg b/setup.cfg index bc92209..d7da790 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,5 +1,5 @@ [metadata] -name = Solvers +name = num_solvers version = 0.1.0 url = https://github.com/syedalimohsinbukhari/Solvers author = Syed Ali Mohsin Bukhari, Astrophysics and Python diff --git a/setup.py b/setup.py index 4ac5866..36f6d3f 100644 --- a/setup.py +++ b/setup.py @@ -5,14 +5,14 @@ readme = f.read() setup( - name='Solvers', + name='num_solvers', version='0.1.0', packages=find_packages(where="src"), url='https://github.com/syedalimohsinbukhari/Solvers', license='MIT', author='Syed Ali Mohsin Bukhari, Astrophysics and Python', author_email='syedali.b@outlook.com, astrophysicsandpython@gmail.com', - description='A library for some basic numerical solvers.', + description='A library for some basic numerical num_solvers.', long_description=readme, long_description_content_type="text/markdown", python_requires=">=3.9", diff --git a/src/Solvers/__init__.py b/src/Solvers/__init__.py deleted file mode 100644 index 4ce4377..0000000 --- a/src/Solvers/__init__.py +++ /dev/null @@ -1,19 +0,0 @@ -"""Created on Dec 20 14:10:28 2023""" - -from typing import Callable, List, Optional, Union - -TOLERANCE = 1e-10 - -DOC_STYLE = 'numpy_napoleon_with_merge' - -IFloat = Union[float, int] - -FList = List[IFloat] -LList = List[FList] -LLList = List[LList] -OptList = Optional[FList] - -FListOrLList = Union[FList, LList] -IFloatOrFList = Union[IFloat, FList] - -OptFunc = Optional[Callable] diff --git a/src/Solvers/integrators/NewtonCotesIntegrators.py b/src/Solvers/integrators/NewtonCotesIntegrators.py deleted file mode 100644 index ecfffe4..0000000 --- a/src/Solvers/integrators/NewtonCotesIntegrators.py +++ /dev/null @@ -1,198 +0,0 @@ -"""Created on Nov 29 14:48:53 2023""" - -from typing import Callable, List, Union - -from ._backend.NEWTON_COTES import SingleVariableNewtonCotesIntegrators as SvNCI - -INTEGRATOR_OUTPUT = Union[List[float], float] - - -def newton_raphson_solver(function: callable, derivative_of_function: callable, initial_guess: float, - tolerance: float = 1e-14) -> float: - """ - Find the root of a function using the Newton-Raphson method. - - Parameters - ---------- - function : callable - The function for which the root is being sought. - derivative_of_function : callable - The derivative of the target function. - initial_guess : float - The initial guess for the root. - tolerance : float, optional - Tolerance for convergence. Default is 1e-14. - - Returns - ------- - float - The approximate root of the function. - """ - f, df, x_0 = function, derivative_of_function, initial_guess - - if abs(f(x_0)) < tolerance: - return x_0 - else: - return newton_raphson_solver(f, df, x_0 - f(x_0) / df(x_0), tolerance) - - -def trapezoid_rule(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using the Trapezoidal Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'trapezoid') - return solver.value_table() if get_table else solver.solve() - - -def simpson_13(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using Simpson's 1/3 Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'simpson_13') - return solver.value_table() if get_table else solver.solve() - - -def simpson_38(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using Simpson's 3/8 Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'simpson_38') - return solver.value_table() if get_table else solver.solve() - - -def boole(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using Boole's Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'boole') - return solver.value_table() if get_table else solver.solve() - - -def weddle_3h10(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using Weddle's 3/10 Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'weddle_3h10') - return solver.value_table() if get_table else solver.solve() - - -def weddle_41h140(function: Callable, x_0: float, x_n: float, composites: int = 1, - get_table: bool = False) -> INTEGRATOR_OUTPUT: - """ - Perform numerical integration using Weddle's 41/140 Rule. - - Parameters - ---------- - function : callable - The function to be integrated. - x_0 : float - The lower limit of integration. - x_n : float - The upper limit of integration. - composites : int, optional - Number of sub-intervals for composite method. Default is 1. - get_table : bool, optional - Whether to return the table of values. Default is False. - - Returns - ------- - INTEGRAL_OUTPUT - The numerical result of integration or the table of values if ``get_table`` is True. - """ - solver = SvNCI(function, x_0, x_n, composites, 'weddle_41h140') - return solver.value_table() if get_table else solver.solve() diff --git a/src/Solvers/integrators/__init__.py b/src/Solvers/integrators/__init__.py deleted file mode 100644 index fd9848c..0000000 --- a/src/Solvers/integrators/__init__.py +++ /dev/null @@ -1,7 +0,0 @@ -"""Created on Dec 20 22:59:44 2023""" - -__all__ = ['newton_raphson_solver', 'trapezoid_rule', 'simpson_13', 'simpson_38', 'boole', 'weddle_3h10', - 'weddle_41h140'] - -from .NewtonCotesIntegrators import (boole, newton_raphson_solver, simpson_13, simpson_38, trapezoid_rule, weddle_3h10, - weddle_41h140) diff --git a/src/Solvers/integrators/_backend/NEWTON_COTES.py b/src/Solvers/integrators/_backend/NEWTON_COTES.py deleted file mode 100644 index 36f8719..0000000 --- a/src/Solvers/integrators/_backend/NEWTON_COTES.py +++ /dev/null @@ -1,117 +0,0 @@ -"""Created on Nov 30 16:37:17 2023""" - -import enum - - -class DataTable: - def __init__(self, x, y): - self.x = x - self.y = y - - def print_table(self, precision=100): - # Find the maximum length of strings in x and y - n_len = 6 - max_len_x = max(len(str(element)) for element in self.x) - max_len_y = max(len(str(element)) for element in self.y) - - # Add extra space for padding - max_len_x += 4 - max_len_y += 4 - n_len += 2 - - print('-' * (n_len + max_len_x + max_len_y + 4)) - # Print the table header - print(f"|{'n':^{n_len}}|{'x':^{max_len_x}}|{'y':^{max_len_y}}|") - - # Print the table border - print('-' * (n_len + max_len_x + max_len_y + 4)) - - # Print the table rows - for i in range(len(self.x)): - print(f"|{str(i + 1):^{n_len}}|{str(self.x[i]):^{max_len_x}}|" - f"{str(round(self.y[i], precision)):^{max_len_y}}|") - - # Print the table border - print('-' * (n_len + max_len_x + max_len_y + 4)) - - -@enum.unique -class Iterators(enum.Enum): - TRAPEZOID = 1 - SIMPSON_13 = 2 - SIMPSON_38 = 3 - BOOLE = 4 - WEDDLE_3H10 = 5.1 - WEDDLE_41H140 = 5.2 - - -class SingleVariableNewtonCotesIntegrators: - - _I = Iterators - - str_iters = {'trapezoid': _I.TRAPEZOID, - 'simpson_13': _I.SIMPSON_13, - 'simpson_38': _I.SIMPSON_38, - 'boole': _I.BOOLE, - 'weddle_3h10': _I.WEDDLE_3H10, - 'weddle_41h140': _I.WEDDLE_41H140} - - nodes = {_I.TRAPEZOID: 2, _I.SIMPSON_13: 2, _I.SIMPSON_38: 3, _I.BOOLE: 4, _I.WEDDLE_3H10: 6, _I.WEDDLE_41H140: 6} - - norm = {_I.TRAPEZOID: lambda x: x / 2, - _I.SIMPSON_13: lambda x: x / 3, - _I.SIMPSON_38: lambda x: (3 * x) / 8, - _I.BOOLE: lambda x: (2 * x) / 45, - _I.WEDDLE_3H10: lambda x: (3 * x) / 10, - _I.WEDDLE_41H140: lambda x: x / 140} - - outer_multipliers = {_I.TRAPEZOID: 1, _I.SIMPSON_13: 1, _I.SIMPSON_38: 1, _I.BOOLE: 7, _I.WEDDLE_3H10: 1, - _I.WEDDLE_41H140: 41} - - inner_multipliers = {_I.TRAPEZOID: [2], - _I.SIMPSON_13: [4, 2], - _I.SIMPSON_38: [3, 3, 2], - _I.BOOLE: [32, 12, 32, 2 * 7], - _I.WEDDLE_3H10: [5, 1, 6, 1, 5, 2 * 1], - _I.WEDDLE_41H140: [216, 27, 272, 27, 216, 2 * 41]} - - def __init__(self, function, x_0, x_n, composites=1, integration_method='weddle'): - self.function = function - self.x_0 = x_0 - self.x_n = x_n - self.solver = self.str_iters[integration_method] if isinstance(integration_method, str) else integration_method - - self.composites = composites * self.nodes[self.solver] - self.dx = (x_n - x_0) / self.composites - - self.n_points = self.nodes[self.solver] - self.norm_value = self.norm[self.solver](self.dx) - self.outer_multiplier_value = self.outer_multipliers[self.solver] - self.inner_multipliers_values = self.inner_multipliers[self.solver] - - def _value_table(self): - x_values = [self.x_0 + i * self.dx for i in range(self.composites + 1)] - return {'x': x_values, 'y': [self.function(i) for i in x_values]} - - def value_table(self, precision=100): - vals = self._value_table() - DataTable(vals['x'], vals['y']).print_table(precision=precision) - - def solve(self, rounded=False, decimal_places=3): - fs = self._value_table()['y'] - - if rounded: - fs = [round(i, decimal_places) for i in fs] - - f_0 = [self.outer_multiplier_value * i for i in [fs[0], fs[-1]]] - - f_list = [] - for multiplier, start in zip(self.inner_multipliers_values, range(1, self.n_points + 1)): - if self.solver == self._I.TRAPEZOID: - f_list.extend(multiplier * i for i in fs[start:-1]) - else: - f_list.extend(multiplier * i for i in fs[start:-1:self.n_points]) - - f_list[0:0] = f_0 - - return self.norm_value * sum(f_list) diff --git a/src/Solvers/integrators/_backend/__init__.py b/src/Solvers/integrators/_backend/__init__.py deleted file mode 100644 index ed70602..0000000 --- a/src/Solvers/integrators/_backend/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Created on Dec 20 23:00:12 2023""" diff --git a/src/Solvers/interpolators/LagrangeInterpolation.py b/src/Solvers/interpolators/LagrangeInterpolation.py deleted file mode 100644 index 74df8c7..0000000 --- a/src/Solvers/interpolators/LagrangeInterpolation.py +++ /dev/null @@ -1,25 +0,0 @@ -"""Created on Nov 01 16:47:03 2023""" - -from ._backend.INTERPOLATION_ import INTERPOLATION - - -class LagrangeInterpolation(INTERPOLATION): - - def __init__(self, given_values, value_to_approximate, function=None, function_values=None): - super().__init__(given_values, value_to_approximate, function, function_values) - - def interpolate(self, give_values=False): - answer = [] - given = self.given_values - for i in range(len(given)): - modified_array = given[:i] + given[i + 1:] - numerator = 1 - denominator = 1 - for value in modified_array: - numerator *= (self.value_to_approximate - value) - denominator *= (given[i] - value) - answer.append((numerator / denominator) * self.function_values[i]) - - result = sum(answer) - - return (result, answer) if give_values else result diff --git a/src/Solvers/interpolators/NewtonianInterpolationMethods.py b/src/Solvers/interpolators/NewtonianInterpolationMethods.py deleted file mode 100644 index eb4b608..0000000 --- a/src/Solvers/interpolators/NewtonianInterpolationMethods.py +++ /dev/null @@ -1,115 +0,0 @@ -"""Created on Nov 02 00:25:14 2023""" - -from typing import Callable, List, Optional - -from ._backend.INTERPOLATION_ import BkwInterpolation, DividedInterpolation, FwdInterpolation, get_result - -L_FLOAT = List[float] - - -def forward_interpolation(given_values: L_FLOAT, value_to_approximate: float, function: Optional[Callable] = None, - function_values: Optional[L_FLOAT] = None, get_difference_table: bool = False): - """ - Perform forward Newton Interpolation. - - Parameters - ---------- - given_values: L_FLOAT - List of x values. - value_to_approximate: float - The x value for which to perform interpolation. - function: Optional[Callable] - The function to use for interpolation. - function_values: Optional[L_FLOAT] - List of function values corresponding to given_values. - get_difference_table: bool - Flag indicating whether to return the difference table. - - Returns - ------- - dict - A dictionary containing - - step_values - - results - Additionally, with ``get_difference_table`` flag, - - difference_table - - Notes - ------ - Either ``function`` or ``function_values`` must be specified. - """ - - interpolation = FwdInterpolation(given_values, value_to_approximate, function, function_values) - return get_result(interpolation, get_difference_table) - - -def backward_interpolation(given_values: L_FLOAT, value_to_approximate: float, function: Optional[Callable] = None, - function_values: Optional[L_FLOAT] = None, get_difference_table: bool = False): - """ - Perform backwards Newton Interpolation. - - Parameters - ---------- - given_values: L_FLOAT - List of x values. - value_to_approximate: float - The x value for which to perform interpolation. - function: Optional[Callable] - The function to use for interpolation. - function_values: Optional[L_FLOAT] - List of function values corresponding to given_values. - get_difference_table: bool - Flag indicating whether to return the difference table. - - Returns - ------- - dict - A dictionary containing - - step_values - - results - Additionally, with ``get_difference_table`` flag, - - difference_table - - Notes - ------ - Either ``function`` or ``function_values`` must be specified. - """ - - interpolation = BkwInterpolation(given_values, value_to_approximate, function, function_values) - return get_result(interpolation, get_difference_table) - - -def divided_interpolation(given_values: L_FLOAT, value_to_approximate: float, function: Optional[Callable] = None, - function_values: Optional[L_FLOAT] = None, get_difference_table: bool = False): - """ - Perform divided Newton Interpolation. - - Parameters - ---------- - given_values: L_FLOAT - List of x values. - value_to_approximate: float - The x value for which to perform interpolation. - function: Optional[Callable] - The function to use for interpolation. - function_values: Optional[L_FLOAT] - List of function values corresponding to given_values. - get_difference_table: bool - Flag indicating whether to return the difference table. - - Returns - ------- - dict - A dictionary containing - - step_values - - results - Additionally, with ``get_difference_table`` flag, - - difference_table - - Notes - ------ - Either ``function`` or ``function_values`` must be specified. - """ - - interpolation = DividedInterpolation(given_values, value_to_approximate, function, function_values) - return get_result(interpolation, get_difference_table) diff --git a/src/Solvers/interpolators/_backend/ERRORS_.py b/src/Solvers/interpolators/_backend/ERRORS_.py deleted file mode 100644 index 4a6036c..0000000 --- a/src/Solvers/interpolators/_backend/ERRORS_.py +++ /dev/null @@ -1,24 +0,0 @@ -"""ERRORS_ - -This module contains custom error classes for interpolation methods. The error classes include, -- InterpolationException: The base exception class for interpolation methods. -- AtLeastOneParameterRequired: This error is raised when neither ``function`` not ``function_values`` parameter is - supplied in an interpolation class. -- WrongBoundaryEquation: This error is raised when the user provides wrong string for boundary equation. - -Created on Oct 29 18:21:27 2023 -""" - -__all__ = ['AtLeastOneParameterRequired', 'WrongBoundaryEquation'] - - -class InterpolationException(BaseException): - pass - - -class AtLeastOneParameterRequired(InterpolationException): - pass - - -class WrongBoundaryEquation(InterpolationException): - pass diff --git a/src/Solvers/interpolators/_backend/INTERPOLATION_.py b/src/Solvers/interpolators/_backend/INTERPOLATION_.py deleted file mode 100644 index 7407554..0000000 --- a/src/Solvers/interpolators/_backend/INTERPOLATION_.py +++ /dev/null @@ -1,237 +0,0 @@ -"""Created on Oct 19 03:46:05 2023""" - -from typing import List, Optional - -from . import ERRORS_ - - -class INTERPOLATION: - - def __init__(self, given_values: List[float], value_to_approximate: float, function: Optional[callable] = None, - function_values: Optional[List[float]] = None): - """ - Initialize the INTERPOLATION object. - - Parameters - ---------- - given_values : List[float] - List of x-values. - value_to_approximate : float - The x-value to interpolate the corresponding y-value. - function : callable, optional - The function to use for calculating y-values if ``function_values`` is not provided. - function_values : List[float], optional - List of y-values corresponding to ``given_values``. If not provided, ``function`` will be used to calculate - these values. - - Raises - ------ - ERRORS_.AtLeastOneParameterRequired - If neither ``function`` nor ``function_values`` is provided. - """ - self.given_values = given_values - self.value_to_approximate = value_to_approximate - - if function is None and function_values is None: - raise ERRORS_.AtLeastOneParameterRequired("One of `function` or `function_values` parameter is required.") - - self.function_values = function_values if function_values else [function(value) for value in given_values] - - def difference_table(self): - """ - Calculate the difference table. - """ - pass - - def interpolate(self): - """ - Interpolate the y-value corresponding to the given x-value. - """ - pass - - -class BaseInterpolation(INTERPOLATION): - - def __init__(self, given_values, value_to_approximate, function=None, function_values=None, n_decimal: int = 15): - """ - Initialize the BaseInterpolation object. - - Parameters - ---------- - given_values : List[float] - List of x-values. - value_to_approximate : float - The x-value to interpolate the corresponding y-value. - function : callable, optional - The function to use for calculating y-values if `function_values` is not provided. - function_values : List[float], optional - List of y-values corresponding to `given_values`. If not provided, `function` will be used to calculate - these values. - n_decimal: int - The number of decimal places to round off to. Default is 15 - """ - - super().__init__(given_values, value_to_approximate, function, function_values) - self.h = given_values[1] - given_values[0] - self.round_ = n_decimal - - def _class_check(self) -> str: - """ - Check the class type. - - Returns - ------- - str - 'Fwd' if the class is FwdInterpolation, 'Bkw' if the class is BkwInterpolation. - """ - return 'Fwd' if self.__class__.__name__ == 'FwdInterpolation' else 'Bkw' - - def _get_index(self) -> int: - """ - Get the index of the value to approximate in the sorted given values. - - Returns - ------- - int - The index of the value to approximate in the sorted given values. - """ - temp = self.given_values - temp.insert(0, self.value_to_approximate) - temp.sort() - - temp_idx = temp.index(self.value_to_approximate) - del temp[temp_idx] - return temp_idx - - def difference_table(self) -> List[float]: - """ - Calculate the divided difference table. - - Returns - ------- - List[float] - The divided difference table. - """ - idx_ = self._get_index() - - table_limit = len(self.given_values) - 1 - - difference_table = [self.function_values] - - for i in range(table_limit): - temp_ = [] - for j, k in zip(difference_table[-1][:-1], difference_table[-1][1:]): - temp_.append(round(k - j, self.round_)) - - difference_table.append(temp_) - - if self._class_check() == 'Fwd': - if idx_ - 1 != 0: - reduced_index = idx_ - 1 - reduced_table = difference_table[:-reduced_index] - d_table = [i[reduced_index] for i in reduced_table] - else: - d_table = difference_table - d_table = [i[0] for i in d_table] - else: - d_table = [v[idx_ - i - 1] for i, v in enumerate(difference_table[:idx_])] - - return d_table - - def interpolate(self): - """ - Interpolate the y-value corresponding to the given x-value. - - Returns - ------- - dict - Dictionary with ``step_values`` (list of interpolated values for each step) and ``result`` - (sum of the interpolated values). - """ - - def find_p(): - return (self.value_to_approximate - self.given_values[self._get_index() - 1]) / self.h - - def factorial(number): - return 1 if number == 0 else number * factorial(number - 1) - - difference_table = self.difference_table()[1:] - initial_value = self.function_values[self._get_index() - 1] - - result, iter_condition = [initial_value], len(difference_table) - - _, p_value = find_p(), [find_p()] - - for i in range(1, iter_condition): - _ *= (find_p() - i) if self._class_check() == 'Fwd' else find_p() + i - p_value.append(_) - - result = [(difference_table[i] * p_value[i]) / factorial(i + 1) for i in range(iter_condition)] - result.insert(0, initial_value) - result = list(map(lambda x: round(x, self.round_), result)) - - return result - - -class FwdInterpolation(BaseInterpolation): - """ - Forward Interpolation class, inheriting from BaseInterpolation. - - This class specializes in implementing the forward interpolation method. - It utilizes the BaseInterpolation class for shared functionality. - - Parameters and methods are inherited from BaseInterpolation. - """ - pass - - -class BkwInterpolation(BaseInterpolation): - """ - Backward Interpolation class, inheriting from BaseInterpolation. - - This class specializes in implementing the backward interpolation method. - It utilizes the BaseInterpolation class for shared functionality. - - Parameters and methods are inherited from BaseInterpolation. - """ - pass - - -class DividedInterpolation(BaseInterpolation): - - def difference_table(self): - - n = len(self.function_values) - difference_table = [[0] * n for _ in range(n)] - - for i in range(n): - difference_table[i][0] = self.function_values[i] - - for j in range(1, n): - for i in range(n - j): - numerator = difference_table[i + 1][j - 1] - difference_table[i][j - 1] - denominator = self.given_values[i + j] - self.given_values[i] - difference_table[i][j] = numerator / denominator - - return difference_table[0] - - def interpolate(self): - n = len(self.function_values) - - product, all_products = 1, [1] - for i in range(1, n): - product *= self.value_to_approximate - self.given_values[i - 1] - all_products.append(product) - - difference_table = self.difference_table() - result = [i * j for i, j in zip(difference_table, all_products)] - - return result - - -def get_result(interpolation, difference_table): - result = interpolation.interpolate() - if difference_table: - result['difference_table'] = interpolation.difference_table() - - return result diff --git a/src/Solvers/interpolators/_backend/__init__.py b/src/Solvers/interpolators/_backend/__init__.py deleted file mode 100644 index 96307e9..0000000 --- a/src/Solvers/interpolators/_backend/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -"""Created on Dec 20 13:35:59 2023""" - -TOLERANCE = 1e-10 diff --git a/src/Solvers/iterative_solvers/non_polynomial.py b/src/Solvers/iterative_solvers/non_polynomial.py deleted file mode 100644 index 07219ce..0000000 --- a/src/Solvers/iterative_solvers/non_polynomial.py +++ /dev/null @@ -1,248 +0,0 @@ -"""Created on Dec 24 12:56:39 2023""" - -from math import sqrt -from typing import Callable, List, Union - -from . import FLOAT_LIST, TOLERANCE - - -def bisection_method(function: Callable, x_0: float, x_n: float, tolerance: float = TOLERANCE, - get_all_values: bool = False) -> FLOAT_LIST: - f, root = function, [] - - def is_sane(): - return f(x_0) < 0 < f(x_n) - - while is_sane(): - x_new = (x_0 + x_n) / 2 - root.append(x_new) - f_x_new = f(x_new) - - if abs(f_x_new) < tolerance: - break - - if f_x_new < 0: - x_0 = x_new - else: - x_n = x_new - - return root if get_all_values else root[-1] - - -def false_position_method(function: Callable, x_0: float, x_n: float, tolerance: float = TOLERANCE, - get_all_values: bool = False) -> FLOAT_LIST: - f, root = function, [] - - def is_sane(): - return f(x_0) < 0 < f(x_n) - - while is_sane(): - f_x0, f_x1 = f(x_0), f(x_n) - - x_new = (x_0 * f_x1 - x_n * f_x0) / (f_x1 - f_x0) - root.append(x_new) - - f_x_new = f(x_new) - - if abs(f_x_new) < tolerance: - break - - if f_x_new == 0: - return root if get_all_values else root[-1] - - if f_x0 * f_x_new < 0: - x_n = x_new - else: - x_0 = x_new - - return root if get_all_values else root[-1] - - -def regula_falsi_method(function: Callable, x_0: float, x_n: float, tolerance: float = TOLERANCE, - get_all_values: bool = False) -> FLOAT_LIST: - return false_position_method(function, x_0, x_n, tolerance, get_all_values) - - -def secant_method(function: Callable, x_0: float, x_n: float, tolerance: float = TOLERANCE, - get_all_values: bool = False) -> FLOAT_LIST: - f, root = function, [] - - while True: - x_new = x_n - f(x_n) * (x_n - x_0) / (f(x_n) - f(x_0)) - - root.append(x_new) - x_0, x_n = x_n, x_new - - f_x_new = f(x_new) - if abs(f_x_new) < tolerance: - break - - if f_x_new == 0: - return root if get_all_values else root[-1] - - return root if get_all_values else root[-1] - - -# Taken form https://en.wikipedia.org/wiki/Muller%27s_method#Computational_example -# minor tweaking applied on variable namings for consistency -def div_diff(f_: Callable, xs_: list[float]): - """Calculate the divided difference f[x0, x1, ...].""" - if len(xs_) == 2: - a, b = xs_ - return (f_(a) - f_(b)) / (a - b) - else: - return (div_diff(f_, xs_[1:]) - div_diff(f_, xs_[0:-1])) / (xs_[-1] - xs_[0]) - - -def muller_method(function: Callable, x0: float, x1: float, x2: float, iterations: int, - get_full_result: bool = True) -> FLOAT_LIST: - """Return the root calculated using Muller's method.""" - - root_ = [x0, x1, x2] - - for _ in range(iterations): - w = div_diff(function, [x2, x1]) + div_diff(function, [x2, x0]) - div_diff(function, [x2, x1]) - s_delta = sqrt(w**2 - 4 * function(x2) * div_diff(function, [x2, x1, x0])) - denominators = [w + s_delta, w - s_delta] - # Take the higher-magnitude denominator - x3 = x2 - 2 * function(x2) / max(denominators, key=abs) - # Advance - x0, x1, x2 = x1, x2, x3 - - root_.append(x3) - - return root_ if get_full_result else root_[-1] - - -def generalized_secant_method(function: Callable, x_0: float, x_1: float, tolerance: float = TOLERANCE, - get_full_result: bool = True) -> FLOAT_LIST: - x_2 = x_1 - function(x_1) / div_diff(function, [x_0, x_1]) - - root_ = [x_0, x_1, x_2] - - while True: - f0 = function(root_[-1]) - - f1 = div_diff(function, [x_2, x_1]) - f2 = div_diff(function, [x_2, x_1, x_0]) * (x_2 - x_1) - - f_ = f1 + f2 - x_3 = x_2 - (f0 / f_) - - x_0, x_1, x_2 = x_1, x_2, x_3 - root_.append(x_3) - - if abs(f0) < tolerance: - break - - return root_ if get_full_result else root_[-1] - - -def sidi_method(function: Callable, x_0: float, x_1: float, tolerance: float = TOLERANCE, - get_full_result: bool = True) -> FLOAT_LIST: - return generalized_secant_method(function, x_0, x_1, tolerance, get_full_result) - - -def ridder_method(function: Callable, x_0: float, x_1: float, tolerance: float = TOLERANCE, - get_full_result: bool = True) -> Union[float, List]: - """ - Find the root of a function using Ridder method. - - Parameters - ---------- - function : Callable - The target function for root finding. - x_0 : float - The lower bound of the initial bracket. - x_1 : float - The upper bound of the initial bracket. - tolerance : float, optional - Tolerance for the convergence criterion. Default is 1e-6. - get_full_result : bool, optional - If True, returns the list of all intermediate roots. If False, returns only the final root. - Default is True. - - Returns - ------- - Union[float, List] - List of roots found during the iterations if ``get_full_result`` is True, else returns the final root. - """ - - def sign_function(functions_to_evaluate, value): - f1, f2 = functions_to_evaluate - return value if f1 - f2 > 0 else -value if f2 - f1 > 0 else 0 - - def is_sane(): - return True if f(x_0) * f(x_1) < 0 else False - - f, root_ = function, [x_0, x_1, (x_0 + x_1) / 2] - - while is_sane(): - x_mid = (x_0 + x_1) / 2 - - num_ = f(x_mid) * (x_mid - x_0) - den_ = sqrt(f(x_mid)**2 - f(x_0) * f(x_1)) - f_ = num_ / den_ - - x_new = x_mid + sign_function([f(x_0), f(x_1)], f_) - - if abs(f(x_new)) < tolerance or f_ == 0: - break - - if f(x_mid) * f(x_new) < 0: - x_0, x_1 = x_mid, x_new - else: - x_0 = x_new - - root_.append(x_new) - - return root_ if get_full_result else root_[-1] - - -def steffensen_method(function: Callable, x_0: float, x_1: float, tolerance: float = TOLERANCE, - get_full_result: bool = True) -> Union[float, List]: - """ - Find the root of a function using Steffensen method. - - Parameters - ---------- - function : Callable - The target function for root finding. - x_0 : float - The lower bound of the initial bracket. - x_1 : float - The upper bound of the initial bracket. - tolerance : float, optional - Tolerance for the convergence criterion. Default is 1e-6. - get_full_result : bool, optional - If True, returns the list of all intermediate roots. If False, returns only the final root. - Default is True. - - Returns - ------- - Union[float, List] - List of roots found during the iterations if ``get_full_result`` is True, else returns the final root. - """ - - def is_sane(): - return True if f(x_0) * f(x_1) < 0 else False - - f, root_ = function, [x_0, x_1, (x_0 + x_1) / 2] - - solve = True if is_sane() else False - - while solve: - f_mid = f(root_[-1]) - f_mid2 = f(root_[-1] + f_mid) - - num_ = f_mid**2 - den_ = f_mid2 - f_mid - - x_new = root_[-1] - (num_ / den_) - - root_.append(x_new) - - if f(x_new) == 0 or abs(f(x_new)) < tolerance: - break - - return root_ if get_full_result else root_[-1] diff --git a/src/Solvers/iterative_solvers/odes/RungeKutta.py b/src/Solvers/iterative_solvers/odes/RungeKutta.py deleted file mode 100644 index dc3c5bf..0000000 --- a/src/Solvers/iterative_solvers/odes/RungeKutta.py +++ /dev/null @@ -1,177 +0,0 @@ -"""Created on Dec 16 03:28:44 2023""" - -from typing import Callable, List, Tuple - -import numpy as np - -NdArray2 = [np.ndarray, np.ndarray] -NdArrayN = [np.ndarray, ...] -RK2_OUTPUT = Tuple[np.ndarray, np.ndarray] - - -# TODO: Multi RK2, RK3 solvers and classes - - -def rk2_solver(ode: Callable, x_initial: float, y_initial: float, step_size: float = 0.1, - x_max: float = 1.0) -> RK2_OUTPUT: - """ - Solve a first-order ordinary differential equation using the Runge-Kutta (RK2) method. - - Parameters - ---------- - ode : Callable - The ordinary differential equation in the form dy/dx = ode(x, y). - x_initial : float - Initial x-value. - y_initial : float - Initial y-value. - step_size : float, optional - Step size for numerical integration, default is 0.1. - x_max : float, optional - Maximum x-value for integration, default is 1.0. - - Returns - ------- - RK2_OUTPUT - Tuple containing arrays of x and y values at each step. - """ - num_steps = int((x_max - x_initial) / step_size) + 1 - x_values = np.zeros(num_steps) - y_values = np.zeros(num_steps) - - x_values[0] = x_initial - y_values[0] = y_initial - - for i in range(1, num_steps): - x_i, y_i = x_values[i - 1], y_values[i - 1] - - k1 = step_size * ode(x_i, y_i) - k2 = step_size * ode(x_i + step_size, y_i + k1) - - y_values[i] = y_i + 0.5 * (k1 + k2) - x_values[i] = x_i + step_size - - return x_values, y_values - - -def rk3_solver(ode: Callable, x_initial: float, y_initial: float, step_size: float = 0.1, - x_max: float = 1.0) -> NdArray2: - """ - Solve a first-order ordinary differential equation using the Runge-Kutta (RK3) method. - - Parameters - ---------- - ode : Callable - The ordinary differential equation in the form dy/dx = ode(x, y). - x_initial : float - Initial x-value. - y_initial : float - Initial y-value. - step_size : float, optional - Step size for numerical integration, default is 0.1. - x_max : float, optional - Maximum x-value for integration, default is 1.0. - - Returns - ------- - NdArray2 - Tuple containing arrays of x and y values at each step. - """ - - num_steps = int((x_max - x_initial) / step_size) + 1 - x_values = np.zeros(num_steps) - y_values = np.zeros(num_steps) - - x_values[0] = x_initial - y_values[0] = y_initial - - for i in range(1, num_steps): - x_i, y_i = x_values[i - 1], y_values[i - 1] - - k1 = step_size * ode(x_i, y_i) - k2 = step_size * ode(x_i + step_size / 2, y_i + k1 / 2) - k3 = step_size * ode(x_i + step_size, y_i - k1 + 2 * k2) - - y_values[i] = y_i + (1 / 6) * (k1 + 4 * k2 + k3) - x_values[i] = x_i + step_size - - return x_values, y_values - - -def rk4_solver(ode: Callable, x_initial: float, y_initial: float, step_size: float = 0.1, - x_max: float = 1.0) -> NdArray2: - - num_steps = int((x_max - x_initial) / step_size) + 1 - x_values = np.zeros(num_steps) - y_values = np.zeros(num_steps) - - x_values[0] = x_initial - y_values[0] = y_initial - - for i in range(1, num_steps): - x_i, y_i = x_values[i - 1], y_values[i - 1] - - k1 = step_size * ode(x_i, y_i) - k2 = step_size * ode(x_i + step_size / 2, y_i + k1 / 2) - k3 = step_size * ode(x_i + step_size / 2, y_i + k2 / 2) - k4 = step_size * ode(x_i + step_size, y_i + k3) - - y_values[i] = y_i + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4) - x_values[i] = x_i + step_size - - return x_values, y_values - - -def rk4_multi_ode(odes: List[Callable], initial_conditions: List[float], step_size: float = 0.1, x_max: float = 1.0, - n_round: int = 9) -> NdArrayN: - """ - Solve a system of first-order ordinary differential equations using the Runge-Kutta (RK4) method. - - Parameters - ---------- - odes : List[Callable] - List of ordinary differential equations in the form dy_i/dx = odes[i](x, y_1, y_2, ..., y_n). - initial_conditions : List[float] - List of initial values for each quantity (x_0, y_1, y_2, ..., y_n). - step_size : float, optional - Step size for numerical integration, default is 0.1. - x_max : float, optional - Maximum x-value for integration, default is 1.0. - n_round : int, optional - Number of decimal places to round the result, default is 9. - - Returns - ------- - NdArrayN - List of arrays containing the values of each quantity at each step. - """ - - num_odes = len(odes) + 1 - num_steps = int((x_max - initial_conditions[0]) / step_size) + 1 - - result = np.zeros((num_steps, num_odes)) - result[0] = initial_conditions - - for i in range(1, num_steps): - quantities = result[i - 1] - m = np.zeros((4, num_odes - 1)) - - for j in range(num_odes - 1): - m[0][j] = step_size * odes[j](*quantities) - - for k in range(1, 3): - temp = [quantities[0] + step_size / 2] + list(quantities[1:] + m[k - 1, :] / 2) - m[k, :] = step_size * np.array([odes[j](*temp) for j in range(num_odes - 1)]) - - temp = [quantities[0] + step_size] + list(quantities[1:] + m[2, :]) - m[3, :] = step_size * np.array([odes[j](*temp) for j in range(num_odes - 1)]) - - for j in range(num_odes): - if j == 0: - result[i, j] = quantities[j] + step_size - else: - temp = m[:, j - 1] - result[i, j] = quantities[j] + (1 / 6) * (temp[0] + 2 * np.sum(temp[1:3]) + temp[3]) - result[i, j] = np.round(result[i, j], n_round) - - return [result[:, i] for i in range(num_odes)] diff --git a/src/Solvers/iterative_solvers/polynomial.py b/src/Solvers/iterative_solvers/polynomial.py deleted file mode 100644 index 35f0e53..0000000 --- a/src/Solvers/iterative_solvers/polynomial.py +++ /dev/null @@ -1,85 +0,0 @@ -"""Created on Dec 20 23:21:44 2023""" - -from cmath import sqrt -from typing import List - -import numpy as np - -TOLERANCE = 1e-8 - - -def generate_random_polynomial(degree: int, low: int = -10, high: int = 10) -> List[float]: - """ - Generate random coefficients for a polynomial. - - Parameters - ---------- - degree : int - The degree of the polynomial. - low : int, optional - The lower bound for the random coefficients. Default is -10. - high : int, optional - The upper bound for the random coefficients. Default is 10. - - Returns - ------- - List[float] - A list of coefficients representing the random polynomial. - """ - - coefficients = np.random.uniform(low, high, degree + 1) - - while coefficients[0] == 0: - coefficients[0] = np.random.uniform(low, high) - - return list(coefficients) - - -def laguerre_method(polynomial: List, x_0: float, degree_of_polynomial: int, tolerance: float = TOLERANCE, - get_full_result: bool = False): - p_val, p_der = np.polyval, np.polyder - root_ = [x_0] - - while True: - poly = p_val(polynomial, x_0) - - if abs(poly) < tolerance: - break - - poly_derivative = p_val(p_der(polynomial), x_0) - poly_double_derivative = p_val(p_der(polynomial, 2), x_0) - - g_ = poly_derivative / poly - h_ = g_**2 - (poly_double_derivative / poly) - - sqrt_in_denominator = sqrt((degree_of_polynomial - 1) * (degree_of_polynomial * h_ - g_**2)) - p1 = max([g_ + sqrt_in_denominator, g_ - sqrt_in_denominator], key=abs) - - x_0 = x_0 - (degree_of_polynomial / p1) - - root_.append(x_0) - - return root_ if get_full_result else root_[-1] - - -def segmented_roots(polynomial: List, x_start: float, x_end: float, degree_of_polynomial: int, num_segments: int = 100, - tolerance: float = TOLERANCE): - x_values = np.linspace(x_start, x_end, num_segments) - all_roots = [] - - for x_0 in x_values: - root = laguerre_method(polynomial, x_0, degree_of_polynomial, tolerance) - - if all(abs(root - existing_root) > tolerance for existing_root in all_roots): - all_roots.append(root) - - unique_roots = np.fromiter(all_roots, dtype=complex) - - new_roots = [] - for i, root_ in enumerate(unique_roots): - if np.allclose(np.imag(root_), 0): - new_roots.append(np.real(root_)) - else: - new_roots.extend([root_, np.conj(root_)]) - - return new_roots diff --git a/src/Solvers/iterative_solvers/system_of_equations/GaussJacobi.py b/src/Solvers/iterative_solvers/system_of_equations/GaussJacobi.py deleted file mode 100644 index 280224d..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/GaussJacobi.py +++ /dev/null @@ -1,52 +0,0 @@ -"""Created on Oct 13 11:37:27 2023""" - -from typing import List, Tuple, Union - -import numpy as np - -from ._backend.ITERATIVE_SOLVER_ import GaussJacobi - -GJ_OUTPUT = Union[np.ndarray, Tuple[np.ndarray, List[np.ndarray]]] -INITIAL_GUESS = Tuple[float, float, float] - - -def gauss_jacobi(system_of_equations: List[List[float]], - solution: List[float], - n_iter: int = 500, - tol: float = 1e-5, - initial_guess: INITIAL_GUESS = (0, 0, 0), - get_solution_set: bool = False) -> GJ_OUTPUT: - """ - Solve a system of linear equations using the Gauss-Jacobi method. - - Parameters - ---------- - system_of_equations : List[List[float]] - Coefficient matrix of the system of linear equations. - solution : List[float] - Right-hand side of the system of linear equations. - n_iter : int, optional - Maximum number of iterations, default is 500. - tol : float, optional - Tolerance for convergence, default is 1e-5. - initial_guess : INITIAL_GUESS, optional - Initial guess for the solution, default is (0, 0, 0). - get_solution_set : bool, optional - If True, returns both the solution and the solution set for each iteration, default is False. - - Returns - ------- - GJ_Output - If ``get_solution_set`` is ``False``, returns the solution as a NumPy array. - If ``get_solution_set`` is ``True``, returns a tuple containing the solution and a list of solution sets for - each iteration. - - Examples - -------- - >>> equations = [[2, 1, 1], [3, 10, 2], [2, 1, 4]] - >>> solutions = [5, 10, 9] - >>> GJ_Solver = gauss_jacobi(equations, solutions) - """ - gj = GaussJacobi(system_of_equations, solution, n_iter, tol, initial_guess) - - return gj.solve() if not get_solution_set else (gj.solve(), gj.solution_set) diff --git a/src/Solvers/iterative_solvers/system_of_equations/GaussSeidel.py b/src/Solvers/iterative_solvers/system_of_equations/GaussSeidel.py deleted file mode 100644 index 95a7dd9..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/GaussSeidel.py +++ /dev/null @@ -1,52 +0,0 @@ -"""Created on Oct 14 05:46:04 2023""" - -from typing import List, Tuple, Union - -import numpy as np - -from ._backend.ITERATIVE_SOLVER_ import GaussSeidel - -GS_OUTPUT = Union[np.ndarray, Tuple[np.ndarray, List[np.ndarray]]] -INITIAL_GUESS = Tuple[float, float, float] - - -def gauss_seidel(system_of_equations: List[List[float]], - solution: List[float], - n_iter: int = 500, - tol: float = 1e-5, - initial_guess: INITIAL_GUESS = (0, 0, 0), - get_solution_set: bool = False) -> GS_OUTPUT: - """ - Solve a system of linear equations using the Gauss-Seidel method. - - Parameters - ---------- - system_of_equations : List[List[float]] - Coefficient matrix of the system of linear equations. - solution : List[float] - Right-hand side of the system of linear equations. - n_iter : int, optional - Maximum number of iterations, default is 500. - tol : float, optional - Tolerance for convergence, default is 1e-5. - initial_guess : INITIAL_GUESS, optional - Initial guess for the solution, default is (0, 0, 0). - get_solution_set : bool, optional - If True, returns both the solution and the solution set for each iteration, default is False. - - Returns - ------- - GS_OUTPUT - If ``get_solution_set`` is ``False``, returns the solution as a NumPy array. - If ``get_solution_set`` is ``True``, returns a tuple containing the solution and a list of solution sets for - each iteration. - - Examples - -------- - >>> equations = [[2, 1, 1], [3, 10, 2], [2, 1, 4]] - >>> solutions = [5, 10, 9] - >>> GS_Solver = gauss_seidel(equations, solutions) - """ - gs = GaussSeidel(system_of_equations, solution, n_iter, tol, initial_guess) - - return gs.solve() if not get_solution_set else (gs.solve(), gs.solution_set) diff --git a/src/Solvers/iterative_solvers/system_of_equations/PredictorCorrectorMethods.py b/src/Solvers/iterative_solvers/system_of_equations/PredictorCorrectorMethods.py deleted file mode 100644 index 4d3104a..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/PredictorCorrectorMethods.py +++ /dev/null @@ -1,105 +0,0 @@ -"""Created on Dec 20 10:55:46 2023""" - -from typing import Callable, Optional, Tuple, Union - -import numpy as np - -from ._backend.PREDICTOR_CORRECTOR import get_x, trapezoidal_rule - -FLOAT = Union[float, np.float32] -EULER_OUTPUT = Union[Tuple[np.ndarray, np.ndarray], float] -EULER_TRAPEZOIDAL_OUTPUT = Union[Tuple[np.ndarray, np.ndarray], Tuple[np.ndarray, Tuple[np.ndarray, np.ndarray]]] - - -def euler_method(ode: Callable, x_initial: FLOAT, y_initial: FLOAT, step_size: Optional[float] = 0.1, - num_steps: Optional[int] = 1, x_max: Optional[float] = None, n_decimal: int = 6, - full_result: bool = False) -> EULER_OUTPUT: - """ - Solve a first-order ordinary differential equation using the Euler method. - - Parameters - ---------- - ode : Callable - The ordinary differential equation in the form dy/dx = ode(x, y). - x_initial : float - Initial x-value. - y_initial : float - Initial y-value. - step_size : float, optional - Step size for numerical integration. If not provided, it is calculated based on num_steps and x_max. - num_steps : int, optional - Number of steps for numerical integration. Default is 1. - x_max : float, optional - Maximum x-value for integration. If not provided, it is calculated based on x_initial, step_size, and num_steps. - n_decimal : int, optional - Number of decimal places to round the result, default is 6. - full_result : bool, optional - If True, return a tuple of arrays (x_values, y_values). If False, return the final y-value only. - - Returns - ------- - EULER_OUTPUT - Result of the Euler method. If full_result is True, returns a tuple of arrays (x_values, y_values). - If full_result is False, returns the final y-value only. - """ - step_size = (x_max - x_initial) / num_steps if step_size is None else step_size - x_max = x_initial + num_steps * step_size if x_max is None else x_max - - y_values, x_values = [y_initial], [x_initial] - - while x_initial < x_max: - y_values.append(y_values[-1] + step_size * ode(x_values[-1], y_values[-1])) - x_values.append(x_values[-1] + step_size) - - x_initial = np.round(x_initial + step_size, n_decimal) - - return (np.array(x_values), np.array(y_values)) if full_result else y_values[-1] - - -def euler_trapezoidal(ode: Callable, x_initial: float, y_initial: float, step_size: float, num_steps: int, - n_decimal: int = 6, get_predictor: bool = False) -> EULER_TRAPEZOIDAL_OUTPUT: - """ - Solve a first-order ordinary differential equation using the Euler method with trapezoidal correction. - - Parameters - ---------- - ode : Callable - The ordinary differential equation in the form dy/dx = ode(x, y). - x_initial : float - Initial x-value. - y_initial : float - Initial y-value. - step_size : float - Step size for numerical integration. - num_steps : int - Number of steps for numerical integration. - n_decimal : int, optional - Number of decimal places to round the result, default is 6. - get_predictor : bool, optional - If True, return both the corrector values and the predictor-corrector pair. - If False, return only the corrector values. - - Returns - ------- - EULER_TRAPEZOIDAL_OUTPUT - Result of the Euler method with trapezoidal correction. - - - If ``get_predictor`` is True, returns (x_values, (predictor_values, corrector_values)). - - If ``get_predictor`` is False, returns (x_values, corrector_values). - """ - x_values = get_x(x_initial, num_steps, step_size, n_decimal) - y_values = np.zeros_like(x_values) - - predictor_values, corrector_values = np.zeros_like(x_values), np.zeros_like(x_values) - y_values[0], predictor_values[0], corrector_values[0] = [y_initial] * 3 - - for i in range(1, num_steps + 1): - predictor = euler_method(ode, x_values[i - 1], corrector_values[i - 1], step_size, n_decimal) - predictor_values[i] = np.round(predictor, n_decimal) - - corrector = trapezoidal_rule(ode, x_values[i - 1], y_values[i - 1], x_values[i], predictor_values[i], - step_size, n_decimal) - corrector_values[i] = corrector - y_values[i] = corrector - - return (x_values, (predictor_values, corrector_values)) if get_predictor else (x_values, corrector_values) diff --git a/src/Solvers/iterative_solvers/system_of_equations/_backend/ITERATIVE_SOLVER_.py b/src/Solvers/iterative_solvers/system_of_equations/_backend/ITERATIVE_SOLVER_.py deleted file mode 100644 index 3658d24..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/_backend/ITERATIVE_SOLVER_.py +++ /dev/null @@ -1,108 +0,0 @@ -"""Created on Oct 14 06:07:34 2023""" - -from typing import List - - -class SysEqnSolver: - def __init__(self, system_of_equations: List[List], solution: List, n_iter: int = 500, tol: float = 1e-5, - initial_guess: tuple = (0, 0, 0)): - self.tol = tol - self.n_iter = n_iter - self.solution = solution - self.initial_guess = initial_guess - self.system_of_equations = system_of_equations - self._arr1, self._arr2, self._arr3 = [self.initial_guess[0]], [self.initial_guess[1]], [self.initial_guess[2]] - - @property - def solution_set(self): - if not self._arr1 or len(self._arr1) == 1: - self.solve() - - return self._arr1, self._arr2, self._arr3 - - def _evaluate(self): - equations, solution, initial_guess = self.system_of_equations, self.solution, self.initial_guess - - iter1_ = solution[0] - iter1_ -= equations[0][1] * initial_guess[1] - iter1_ -= equations[0][2] * initial_guess[2] - iter1_ = iter1_ * equations[0][0]**-1 - - iter2_ = solution[1] - - if self.__class__.__name__ == 'GaussJacobi': - iter2_ -= equations[1][0] * initial_guess[0] - else: - iter2_ -= equations[1][0] * iter1_ - - iter2_ -= equations[1][2] * initial_guess[2] - iter2_ = iter2_ * equations[1][1]**-1 - - iter3_ = solution[2] - - if self.__class__.__name__ == 'GaussJacobi': - iter3_ -= equations[2][0] * initial_guess[0] - iter3_ -= equations[2][1] * initial_guess[1] - else: - iter3_ -= equations[2][0] * iter1_ - iter3_ -= equations[2][1] * iter2_ - - iter3_ = iter3_ * equations[2][2]**-1 - - self._arr1.append(iter1_) - self._arr2.append(iter2_) - self._arr3.append(iter3_) - - return [iter1_, iter2_, iter3_] - - # def _string(self): - # equations, solution, initial_guess = self.system_of_equations, self.solution, self.initial_guess - # str1_ = f'({solution[0][0]}' - # str1_ += f' - ({equations[0][1]}*{initial_guess[1]})' - # str1_ += f' - ({equations[0][2]}*{initial_guess[2]}))' - # - # str1_ = f'({Fraction(1, equations[0][0])})*{str1_}' - # - # ic(str1_, eval(str1_)) - # - # str2_ = f'({solution[1][0]}' - # - # if self.__class__.__name__ == 'GaussJacobi': - # str2_ += f' - ({equations[1][0]}*{initial_guess[0]})' - # else: - # str2_ += f' - ({equations[1][0]}*{eval(str1_)})' - # - # str2_ += f' - ({equations[1][2]}*{initial_guess[2]}))' - # - # str2_ = f'({Fraction(1, equations[1][1])})*{str2_}' - # - # ic(str2_, eval(str2_)) - # - # str3_ = f'({solution[2][0]}' - # - # if self.__class__.__name__ == 'GaussJacobi': - # str3_ += f' - ({equations[2][0]}*{initial_guess[0]})' - # str3_ += f' - ({equations[2][1]}*{initial_guess[1]}))' - # else: - # str3_ += f' - ({equations[2][0]}*{eval(str1_)})' - # str3_ += f' - ({equations[2][1]}*{eval(str2_)}))' - # - # str3_ = f'({Fraction(1, equations[2][2])})*{str3_}' - # - # ic(str3_, eval(str3_)) - - def solve(self): - for iter_ in range(self.n_iter): - self.initial_guess = self._evaluate() - - # self.initial_guess = [round(i) for i in self.initial_guess] - - return self.initial_guess - - -class GaussJacobi(SysEqnSolver): - pass - - -class GaussSeidel(SysEqnSolver): - pass diff --git a/src/Solvers/iterative_solvers/system_of_equations/_backend/PREDICTOR_CORRECTOR.py b/src/Solvers/iterative_solvers/system_of_equations/_backend/PREDICTOR_CORRECTOR.py deleted file mode 100644 index 54cc44f..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/_backend/PREDICTOR_CORRECTOR.py +++ /dev/null @@ -1,23 +0,0 @@ -"""Created on Dec 20 12:59:39 2023""" - -import numpy as np - - -def _trapezoidal(ode, x_, y_, x_n1, y_n1, h): - f_n = ode(x_, y_) - f_p = ode(x_n1, y_n1) - return y_ + (h / 2) * (f_n + f_p) - - -def trapezoidal_rule(ode, x_, y_, x_n1, y_n1, h, n_decimal): - result = _trapezoidal(ode, x_, y_, x_n1, y_n1, h) - return np.round(result, n_decimal) - - -def get_x(x_0, n, h, n_decimal): - def np_round(x_): - return np.round(x_, n_decimal) - - x_n = map(np_round, np.linspace(x_0, x_0 + n * h, n + 1)) - x_n = np.array(list(x_n)) - return np_round(x_n) diff --git a/src/Solvers/iterative_solvers/system_of_equations/_backend/__init__.py b/src/Solvers/iterative_solvers/system_of_equations/_backend/__init__.py deleted file mode 100644 index 446b833..0000000 --- a/src/Solvers/iterative_solvers/system_of_equations/_backend/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Created on Dec 20 13:36:23 2023""" diff --git a/src/num_solvers/__backend/__init__.py b/src/num_solvers/__backend/__init__.py new file mode 100644 index 0000000..14fe82b --- /dev/null +++ b/src/num_solvers/__backend/__init__.py @@ -0,0 +1 @@ +"""Created on Jan 07 16:05:15 2024""" diff --git a/src/num_solvers/__backend/errors_.py b/src/num_solvers/__backend/errors_.py new file mode 100644 index 0000000..ad2a1bc --- /dev/null +++ b/src/num_solvers/__backend/errors_.py @@ -0,0 +1,75 @@ +"""Custom errors module + +This module contains custom error classes for interpolation methods. The error classes include, + +- NumSolversErrors: + The base exception class for this module. + +# Interpolation Errors +- AtLeastOneParameterRequired: + When neither ``function`` not ``function_values`` parameter is supplied for interpolation. +- WrongBoundaryEquation: + When the user provides wrong string for boundary equation. + +# Matrix Errors +- NonSymmetricMatrix: + When the matrix provided isn't symmetric. +- NotPositiveDefinite: + When the matrix provided isn't positive definite. + +# Polynomial Errors +- DegreeOfPolynomialNotCorrect: + When length of the polynomial and the degree doesn't match. + +# NumSolversMatrix Errors +- IndexCanNotBeSlice: + When the index to be set is slice instead of integer. +- InconsistentDimensions: + When the members of the array have inconsistent length. + +Created on Jan 07 16:05:26 2024 +""" + +__all__ = ['AtLeastOneParameterRequired', 'WrongBoundaryEquation', 'NonSymmetricMatrix', 'NumStepCanNotBeFloat', + 'NotPositiveDefinite', 'DegreeOfPolynomialNotCorrect', 'IndexCanNotBeSlice', 'InconsistentDimensions', + 'XFindNotDefined'] + + +class NumSolversErrors(BaseException): + pass + + +class AtLeastOneParameterRequired(NumSolversErrors): + pass + + +class WrongBoundaryEquation(NumSolversErrors): + pass + + +class NumStepCanNotBeFloat(NumSolversErrors): + pass + + +class XFindNotDefined(NumSolversErrors): + pass + + +class NonSymmetricMatrix(NumSolversErrors): + pass + + +class NotPositiveDefinite(NumSolversErrors): + pass + + +class DegreeOfPolynomialNotCorrect(NumSolversErrors): + pass + + +class IndexCanNotBeSlice(NumSolversErrors): + pass + + +class InconsistentDimensions(NumSolversErrors): + pass diff --git a/src/num_solvers/__backend/extra_.py b/src/num_solvers/__backend/extra_.py new file mode 100644 index 0000000..43d8012 --- /dev/null +++ b/src/num_solvers/__backend/extra_.py @@ -0,0 +1,155 @@ +"""Extra functionalities + +This module holds general functionalities that can be used in the entire package, + +- NumSolversMatrix: + Introduced a list-type array object for num_solvers package +- round_value_: + Rounds the given value. +- round_list_: + Rounds the given list. +- linear_list: + Provides a list of linearly spaced elements from start to stop. + +Created on Jan 10 00:01:13 2024 +""" + +__all__ = ['linear_list', 'round_value_', 'round_list_', 'NumSolversMatrix'] + +from .errors_ import AtLeastOneParameterRequired, InconsistentDimensions, IndexCanNotBeSlice +from .. import FList, IFloat, LList, OptIFloat + + +class NumSolversMatrix: + """Class for NumSolversMatrix object.""" + + def __init__(self, values: LList or FList): + """ + Initializer class for NumSolverArray. + + Parameters + ---------- + values: + Values to transform into a NumSolverArray object. + """ + + self.values = values + self.__sanity_check() + + def __sanity_check(self): + """Performs a sanity check on NumSolversMatrix object.""" + + if isinstance(self.values[0], list): + len1 = len(self.values[0]) + if all(len(x) == len1 for x in self.values[1:]): + pass + else: + raise InconsistentDimensions('The members of NumSolversMatrix are not of same length.') + + def __repr__(self): + vals = self.values + len_max = len(str(max(max(vals)))) if isinstance(vals[0], list) else len(str(max(vals))) + len_cond = len_max + 5 if len_max % 2 == 0 else len_max + 4 + + out = '[' + for i, v in enumerate(vals): + if isinstance(v, (int, float)): + out += f'{v:^{len_cond}.{len_cond - 4}f}' if i + 1 == len(vals) else f'{v:^{len_cond}.{len_cond - 4}f}|' + else: + for k in v: + out += f'{k:^{len_cond}.5f}' if i + 1 == len(v) else f'{k:^{len_cond}.5f}|' + if v != vals[-1]: + out += ']\n[' + out += ']' + + return out + + def __getitem__(self, item): + return self.values[item][0] if isinstance(item, slice) else self.values[item] + + def __setitem__(self, index, value): + if isinstance(self.values[index], list): + raise IndexCanNotBeSlice() + self.values[index] = value + + +def round_value_(value: IFloat, n_decimal: int = 8) -> IFloat: + """ + Rounds off a given value. + + Parameters + ---------- + value: + The value to round off. + n_decimal: + Number of decimal places to round off to. + + Returns + ------- + Rounded off value. + """ + + return round(value, n_decimal) + + +def round_list_(values: FList or LList, n_decimal: int = 8) -> FList or LList: + """ + Maps the round function to a list. + + Parameters + ---------- + values: + List of values to round off. + n_decimal: + Number of decimal places to round off to. + + Returns + ------- + List of rounded floats. + """ + + def _auxiliary_function(list_of_values): + return list(map(lambda x: round_value_(x, n_decimal), list_of_values)) + + if not isinstance(values[0], list): + return _auxiliary_function(values) + else: + return [_auxiliary_function(i) for i in values] + + +def linear_list(start: IFloat, stop: IFloat, n_elements: OptIFloat = None, step_size: OptIFloat = None, + n_decimal: int = 8) -> FList: + """ + Creates a linearly spaced list of floats. + + Parameters + ---------- + start: + Starting value for the list. + stop: + Stopping value for the list. + n_elements: + Number of elements for the list + step_size: + Step size for numerical integration, default is None. + n_decimal: + Number of digits to round off the elements of the output list. + + Returns + ------- + Linearly spaced list of floats between ``start`` and ``stop`` + """ + + if step_size is None and n_elements is None: + raise AtLeastOneParameterRequired('Either `step_size` or `n_elements` is required.') + + if step_size is None: + step_size = (stop - start) / (n_elements - 1) + elif n_elements is None: + n_elements = (stop - start) / step_size + n_elements += 1 + n_elements = int(n_elements) + + value = [start + (i * step_size) for i in range(n_elements)] + + return round_list_(value, n_decimal) diff --git a/src/num_solvers/__backend/interpolation_.py b/src/num_solvers/__backend/interpolation_.py new file mode 100644 index 0000000..ceabde4 --- /dev/null +++ b/src/num_solvers/__backend/interpolation_.py @@ -0,0 +1,252 @@ +"""Interpolation base module + +This module provides the basic class structure for various interpolation methods, such as, + +- FwdInterpolation: A class to perform Newton's forwards interpolation. +- BkwInterpolation: A class to perform Newton's backwards interpolation. +- DividedInterpolation: A class to perform Newton's divided difference interpolation. +- LagrangeInterpolation: A class to perform the Lagrange interpolation method. + +Created on Oct 19 03:46:05 2023 +""" + +__all__ = ['Interpolation', 'FwdInterpolation', 'BkwInterpolation', 'DividedInterpolation', 'LagrangeInterpolation'] + +from custom_inherit import doc_inherit + +from .. import DOC_STYLE, FList, IFloat, LList, N_DECIMAL, OptFunc, OptList +from ..__backend.errors_ import AtLeastOneParameterRequired +from ..__backend.extra_ import round_list_ + + +class Interpolation: + """ + Base interpolation class. + + Methods + ------- + difference_table: + To calculate the difference table for the interpolation method. + interpolate: + To interpolate the y-value corresponding to the given x-value. + """ + + def __init__(self, given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None, n_decimal: int = N_DECIMAL): + """ + Initialize the INTERPOLATION object. + + Parameters + ---------- + given_values: + List of x-values. + value_to_approximate: + The x-value to interpolate the corresponding y-value. + function: + The function to use for calculating y-values if ``function_values`` is not provided. + function_values: + List of y-values corresponding to ``given_values``. If not provided, ``function`` will be used to calculate + these values. + n_decimal: + Number of digits to round off to. Default is 8. + + Raises + ------ + ERRORS_.AtLeastOneParameterRequired + If neither ``function`` nor ``function_values`` is provided. + """ + + self.given_values = given_values + self.value_to_approximate = value_to_approximate + + if function is None and function_values is None: + raise AtLeastOneParameterRequired("One of `function` or `function_values` parameter is required.") + + self.function_values = function_values if function_values else [function(value) for value in given_values] + self.step_size = given_values[1] - given_values[0] + self.round_ = n_decimal + + def _class_check(self) -> str: + """ + Check the class type. + + Returns + ------- + 'Fwd' if the class is FwdInterpolation, 'Bkw' if the class is BkwInterpolation. + """ + + return 'Fwd' if self.__class__.__name__ == 'FwdInterpolation' else 'Bkw' + + def _get_index(self) -> int: + """ + Get the index of the value to approximate in the sorted given values. + + Returns + ------- + The index of the value to approximate in the sorted given values. + """ + + temp = self.given_values + temp.insert(0, self.value_to_approximate) + temp.sort() + + temp_idx = temp.index(self.value_to_approximate) + del temp[temp_idx] + + return temp_idx + + def difference_table(self) -> FList: + """ + Calculate the divided difference table. + + Returns + ------- + The difference table. + """ + + idx_ = self._get_index() + + table_limit = len(self.given_values) - 1 + + difference_table = [self.function_values] + + for i in range(table_limit): + temp_ = [] + for j, k in zip(difference_table[-1][:-1], difference_table[-1][1:]): + temp_.append(round(k - j, self.round_)) + + difference_table.append(temp_) + + if self._class_check() == 'Fwd': + if idx_ - 1 != 0: + reduced_index = idx_ - 1 + reduced_table = difference_table[:-reduced_index] + d_table = [i[reduced_index] for i in reduced_table] + else: + d_table = difference_table + d_table = [i[0] for i in d_table] + else: + d_table = [v[idx_ - i - 1] for i, v in enumerate(difference_table[:idx_])] + + return d_table + + def interpolate(self) -> IFloat: + """ + To interpolate the y-value corresponding to the given x-value. + + Returns + ------- + The interpolated value. + """ + + def find_p(): + return (approx - given[index - 1]) / step_size + + def factorial(number): + return 1 if number == 0 else number * factorial(number - 1) + + approx, index, step_size = self.value_to_approximate, self._get_index(), self.step_size + given, func_values, round_to, class_ = self.given_values, self.function_values, self.round_, self._class_check() + + difference_table = self.difference_table()[1:] + initial_value = func_values[index - 1] + + result, iter_condition = [initial_value], len(difference_table) + + _, p_value = find_p(), [find_p()] + + for i in range(1, iter_condition): + _ *= (find_p() - i) if class_ == 'Fwd' else find_p() + i + p_value.append(_) + + result = [(difference_table[i] * p_value[i]) / factorial(i + 1) for i in range(iter_condition)] + result.insert(0, initial_value) + result = sum(round_list_(result, round_to)) + + return result + + +@doc_inherit(Interpolation, style=DOC_STYLE) +class FwdInterpolation(Interpolation): + """Forward Interpolation class specializes in implementing the forward difference method.""" + + +@doc_inherit(Interpolation, style=DOC_STYLE) +class BkwInterpolation(Interpolation): + """Backward Interpolation class specializes in implementing the backward difference method.""" + + +@doc_inherit(Interpolation, style=DOC_STYLE) +class DividedInterpolation(Interpolation): + """Divided Interpolation class specializes in implementing the divided difference method.""" + + def __init__(self, given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None, n_decimal: int = N_DECIMAL): + super().__init__(given_values, value_to_approximate, function, function_values, n_decimal) + + def difference_table(self) -> FList: + given, func_values = self.given_values, self.function_values + n_values = len(func_values) + + difference_table: LList = [[0] * n_values for _ in range(n_values)] + + for i in range(n_values): + difference_table[i][0] = func_values[i] + + for j in range(1, n_values): + for i in range(n_values - j): + numerator = difference_table[i + 1][j - 1] - difference_table[i][j - 1] + denominator = given[i + j] - given[i] + difference_table[i][j] = numerator / denominator + + return difference_table[0] + + def interpolate(self) -> IFloat: + approx, given, func_values = self.value_to_approximate, self.given_values, self.function_values + n_values, round_to = len(func_values), self.round_ + + product, all_products = 1, [1] + for i in range(1, n_values): + product *= approx - given[i - 1] + all_products.append(product) + + difference_table = self.difference_table() + result = round(sum(i * j for i, j in zip(difference_table, all_products)), round_to) + + return result + + +@doc_inherit(Interpolation, style=DOC_STYLE) +class LagrangeInterpolation(Interpolation): + """Class to implement Lagrange interpolation technique. + + Methods + ------- + + difference_table: + None + """ + + def __init__(self, given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None, n_decimal: int = N_DECIMAL): + super().__init__(given_values, value_to_approximate, function, function_values, n_decimal) + + def difference_table(self) -> None: + """Difference table not implemented for Lagrange Interpolation Method.""" + + def interpolate(self) -> IFloat: + answer, given = [], self.given_values + approx, func_values = self.value_to_approximate, self.function_values + + for index, value_ in enumerate(given): + modified_array = given[:index] + given[index + 1:] + numerator, denominator = 1, 1 + + for value in modified_array: + numerator *= (approx - value) + denominator *= (value_ - value) + + answer_ = (numerator / denominator) * func_values[index] + answer.append(round(answer_, self.round_)) + + return sum(answer) diff --git a/src/num_solvers/__backend/iterative_solvers_.py b/src/num_solvers/__backend/iterative_solvers_.py new file mode 100644 index 0000000..ed8b053 --- /dev/null +++ b/src/num_solvers/__backend/iterative_solvers_.py @@ -0,0 +1,191 @@ +"""Iterative Solvers + +This module provides functionality for iterative matrix solvers, such as, + +- GaussJacobi +- GaussSeidel + +using the base class, + +- SysEqnSolver + +Currently, this module supports only 3x3, diagonally dominant matrices. + +Created on Oct 14 06:07:34 2023 +""" + +__all__ = ['SysEqnSolver', 'GaussJacobi', 'GaussSeidel'] + +from .extra_ import round_list_ +from .. import FList, FTuple, IFloat, LLList, LList, N_DECIMAL, TOLERANCE + + +class SysEqnSolver: + """Base class for system of equation solvers.""" + + def __init__(self, system_of_equations: LList, solution: FList, n_iter: int = 500, tolerance: IFloat = TOLERANCE, + initial_guess: FTuple = (0, 0, 0), n_decimal: int = N_DECIMAL): + """ + Initializer for SysEqnSolver class. + + Parameters + ---------- + system_of_equations: + A list of system of equations. + solution: + The solution to the provided system of equations. + n_iter: + Number of iterations to perform. Defaults to 500. + tolerance: + The tolerance level for solution. Defaults to 1e-8. + initial_guess: + The initial guess for the solutions to the system of equations. Defaults to (0, 0, 0). + n_decimal: + Number of digits to round off to. Default is 8. + """ + + self.n_iter = n_iter + self.solution = solution + self.n_decimal = n_decimal + self.tolerance = tolerance + self.initial_guess = initial_guess + self.system_of_equations = system_of_equations + self._arr1, self._arr2, self._arr3 = [self.initial_guess[0]], [self.initial_guess[1]], [self.initial_guess[2]] + + self.__sanity_check() + + def __sanity_check(self): + def _auxiliary(val1, val2, val3): + return abs(val1) >= abs(val2) + abs(val3) + + sys_eq = self.system_of_equations + cond1 = _auxiliary(sys_eq[0][0], sys_eq[0][1], sys_eq[0][2]) + cond2 = _auxiliary(sys_eq[1][1], sys_eq[1][0], sys_eq[1][2]) + cond3 = _auxiliary(sys_eq[2][2], sys_eq[2][0], sys_eq[2][1]) + + if cond1 and cond2 and cond3: + pass + else: + print('The system of equations is not diagonally dominant. The solutions might not be correct.') + + @property + def solution_set(self) -> LList: + """ + Gives the solution set, e.g., all the iterations for the system of equations. + + Returns + ------- + List of all iterations for the solution. + """ + + if not self._arr1 or len(self._arr1) == 1: + self.solve() + + return [self._arr1, self._arr2, self._arr3] + + def _evaluate(self): + class_name, tolerance = self.__class__.__name__, self.tolerance + equations, solution, initial_guess = self.system_of_equations, self.solution, self.initial_guess + + iter1_ = solution[0] + iter1_ -= equations[0][1] * initial_guess[1] + iter1_ -= equations[0][2] * initial_guess[2] + iter1_ = iter1_ * equations[0][0]**-1 + + iter2_ = solution[1] + + if class_name == 'GaussJacobi': + iter2_ -= equations[1][0] * initial_guess[0] + else: + iter2_ -= equations[1][0] * iter1_ + + iter2_ -= equations[1][2] * initial_guess[2] + iter2_ = iter2_ * equations[1][1]**-1 + + iter3_ = solution[2] + + if class_name == 'GaussJacobi': + iter3_ -= equations[2][0] * initial_guess[0] + iter3_ -= equations[2][1] * initial_guess[1] + else: + iter3_ -= equations[2][0] * iter1_ + iter3_ -= equations[2][1] * iter2_ + + iter3_ = iter3_ * equations[2][2]**-1 + + self._arr1.append(iter1_) + self._arr2.append(iter2_) + self._arr3.append(iter3_) + + return [iter1_, iter2_, iter3_] + + def __break_iteration(self): + cond1 = abs(self._arr1[-1] - self._arr1[-2]) < self.tolerance + cond2 = abs(self._arr2[-1] - self._arr2[-2]) < self.tolerance + cond3 = abs(self._arr3[-1] - self._arr3[-2]) < self.tolerance + + return cond1 and cond2 and cond3 + + # def _string(self): + # equations, solution, initial_guess = self.system_of_equations, self.solution, self.initial_guess + # str1_ = f'({solution[0][0]}' + # str1_ += f' - ({equations[0][1]}*{initial_guess[1]})' + # str1_ += f' - ({equations[0][2]}*{initial_guess[2]}))' + # + # str1_ = f'({Fraction(1, equations[0][0])})*{str1_}' + # + # ic(str1_, eval(str1_)) + # + # str2_ = f'({solution[1][0]}' + # + # if self.__class__.__name__ == 'GaussJacobi': + # str2_ += f' - ({equations[1][0]}*{initial_guess[0]})' + # else: + # str2_ += f' - ({equations[1][0]}*{eval(str1_)})' + # + # str2_ += f' - ({equations[1][2]}*{initial_guess[2]}))' + # + # str2_ = f'({Fraction(1, equations[1][1])})*{str2_}' + # + # ic(str2_, eval(str2_)) + # + # str3_ = f'({solution[2][0]}' + # + # if self.__class__.__name__ == 'GaussJacobi': + # str3_ += f' - ({equations[2][0]}*{initial_guess[0]})' + # str3_ += f' - ({equations[2][1]}*{initial_guess[1]}))' + # else: + # str3_ += f' - ({equations[2][0]}*{eval(str1_)})' + # str3_ += f' - ({equations[2][1]}*{eval(str2_)}))' + # + # str3_ = f'({Fraction(1, equations[2][2])})*{str3_}' + # + # ic(str3_, eval(str3_)) + + def solve(self) -> LLList: + """ + Solves the system of equations. + + Returns + ------- + Solutions for the provided system of equation. + """ + + n_iter, n_decimal = self.n_iter, self.n_decimal + + for iter_ in range(n_iter): + self.initial_guess = self._evaluate() + if self.__break_iteration(): + break + + return round_list_(self.initial_guess, n_decimal) + + +class GaussJacobi(SysEqnSolver): + """This class implements the Gauss-Jacobi method of solving system of equations.""" + pass + + +class GaussSeidel(SysEqnSolver): + """This class implements the Gauss-Seidel method of solving system of equations.""" + pass diff --git a/src/num_solvers/__backend/newton_cotes_.py b/src/num_solvers/__backend/newton_cotes_.py new file mode 100644 index 0000000..f87c256 --- /dev/null +++ b/src/num_solvers/__backend/newton_cotes_.py @@ -0,0 +1,216 @@ +"""Newton-Cotes integration base module. + +This module provides basic integration functionality of various Newton-Cotes integrators. Along with that, the module +provides a class for pretty-printing the integrator tables. The classes used in this module are, + +- DataTable: Provides the pretty printing of integrator data tables. +- Iterators: An enumeration class for the integrators. +- SingleVariableNewtonCotesIntegrators: The base class for NC integrators. + +Along with these classes, this module also provides a helper function to the newton_cotes_integrators front-end module + +- solve_helper: Wraps the solver into a function with solution and data table. + +Created on Nov 30 16:37:17 2023 +""" + +__all__ = ['SingleVariableNewtonCotesIntegrators', 'solve_helper'] + +import enum + +from custom_inherit import doc_inherit + +from .extra_ import round_list_ +from .. import DOC_STYLE, FList, Func, IFloat, N_DECIMAL + + +class DataTable: + """Class for making integrator data tables.""" + + def __init__(self, x: FList, y: FList): + """ + Initializer for DataTable class + + Parameters + ---------- + x: + List of x values. + y: + List of y values. + """ + self.x = x + self.y = y + + def print_table(self): + """Prints the data table.""" + + # Find the maximum length of strings in x and y + n_len = 6 + max_len_x = max(len(str(element)) for element in self.x) + max_len_y = max(len(str(element)) for element in self.y) + + # Add extra space for padding + max_len_x += 10 + max_len_y += 10 + n_len += 2 + + print('-' * (n_len + max_len_x + max_len_y + 4)) + # Print the table header + print(f"|{'n':^{n_len}}|{'x':^{max_len_x}}|{'y':^{max_len_y}}|") + + # Print the table border + print('-' * (n_len + max_len_x + max_len_y + 4)) + + # Print the table rows + for index, value in enumerate(self.x): + print(f"|{str(index + 1):^{n_len}}|{str(value):^{max_len_x}}|" + f"{str(self.y[index]):^{max_len_y}}|") + + # Print the table border + print('-' * (n_len + max_len_x + max_len_y + 4)) + + +@enum.unique +class Iterators(enum.Enum): + """Iterator enum class.""" + + TRAPEZOID = 1 + SIMPSON_13 = 2 + SIMPSON_38 = 3 + BOOLE = 4 + WEDDLE_3H10 = 5.1 + WEDDLE_41H140 = 5.2 + + +class SingleVariableNewtonCotesIntegrators: + """Provides functionality to integrate single variables using Newton-Cotes methods.""" + + _I = Iterators + + str_iters = {'trapezoid': _I.TRAPEZOID, + 'simpson_13': _I.SIMPSON_13, + 'simpson_38': _I.SIMPSON_38, + 'boole': _I.BOOLE, + 'weddle_3h10': _I.WEDDLE_3H10, + 'weddle_41h140': _I.WEDDLE_41H140} + + nodes = {_I.TRAPEZOID: 2, _I.SIMPSON_13: 2, _I.SIMPSON_38: 3, _I.BOOLE: 4, _I.WEDDLE_3H10: 6, _I.WEDDLE_41H140: 6} + + norm = {_I.TRAPEZOID: lambda x: x / 2, + _I.SIMPSON_13: lambda x: x / 3, + _I.SIMPSON_38: lambda x: (3 * x) / 8, + _I.BOOLE: lambda x: (2 * x) / 45, + _I.WEDDLE_3H10: lambda x: (3 * x) / 10, + _I.WEDDLE_41H140: lambda x: x / 140} + + outer_multipliers = {_I.TRAPEZOID: 1, _I.SIMPSON_13: 1, _I.SIMPSON_38: 1, _I.BOOLE: 7, _I.WEDDLE_3H10: 1, + _I.WEDDLE_41H140: 41} + + inner_multipliers = {_I.TRAPEZOID: [2], + _I.SIMPSON_13: [4, 2], + _I.SIMPSON_38: [3, 3, 2], + _I.BOOLE: [32, 12, 32, 2 * 7], + _I.WEDDLE_3H10: [5, 1, 6, 1, 5, 2 * 1], + _I.WEDDLE_41H140: [216, 27, 272, 27, 216, 2 * 41]} + + def __init__(self, function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, + integration_method: str = 'weddle'): + """ + Initializer for Newton-Cotes single variable integrator class. + + Parameters + ---------- + function : + The function to be integrated. + x_0 : + The lower limit of integration. + x_n : + The upper limit of integration. + composites : + Number of sub-intervals for composite method. Default is 1. + integration_method: + The integration method to use for evaluation. + """ + self.function = function + self.x_0 = x_0 + self.x_n = x_n + self.solver = self.str_iters[integration_method] if isinstance(integration_method, str) else integration_method + + _solver = self.solver + + self.composites = composites * self.nodes[_solver] + self.dx = (x_n - x_0) / self.composites + + self.n_points = self.nodes[_solver] + self.norm_value = self.norm[_solver](self.dx) + self.outer_multiplier_value = self.outer_multipliers[_solver] + self.inner_multipliers_values = self.inner_multipliers[_solver] + + def _value_table(self): + x, dx, func, composites = self.x_0, self.dx, self.function, self.composites + + x_values = [x + i * dx for i in range(composites + 1)] + y_values = [func(i) for i in x_values] + + return x_values, y_values + + def value_table(self, n_decimal: int = N_DECIMAL): + """ + Prints the data table for the given values. + + Parameters + ---------- + n_decimal: + The number of decimal places to round up to. Default is 8. + """ + + vals = self._value_table() + + x_values = round_list_(vals[0], n_decimal) + y_values = round_list_(vals[1], n_decimal) + + DataTable(x_values, y_values).print_table() + + def solve(self) -> IFloat: + """ + Solves the Newton-Cotes integrator. + + Returns + ------- + The integration value for the provided function. + """ + + value_table, solver_method, normalizer = self._value_table()[1], self.solver, self.norm_value + inner_mul, outer_mul, n_points = self.inner_multipliers_values, self.outer_multiplier_value, self.n_points + + f_0 = [outer_mul * i for i in [value_table[0], value_table[-1]]] + + f_list = [] + for multiplier, start in zip(inner_mul, range(1, n_points + 1)): + if solver_method == self._I.TRAPEZOID: + f_list.extend(multiplier * i for i in value_table[start:-1]) + else: + f_list.extend(multiplier * i for i in value_table[start:-1:n_points]) + + f_list[0:0] = f_0 + + return normalizer * sum(f_list) + + +@doc_inherit(SingleVariableNewtonCotesIntegrators.__init__, style=DOC_STYLE) +def solve_helper(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False, + method: str = 'trapezoid'): + """ + Helper function for performing numerical integration. + + Returns + ------- + The numerical result of integration. + """ + + solver = SingleVariableNewtonCotesIntegrators(function, x_0, x_n, composites, method) + + if get_table: + solver.value_table() + + return solver.solve() diff --git a/src/num_solvers/__backend/root_finding_algorithms_.py b/src/num_solvers/__backend/root_finding_algorithms_.py new file mode 100644 index 0000000..e16f102 --- /dev/null +++ b/src/num_solvers/__backend/root_finding_algorithms_.py @@ -0,0 +1,26 @@ +"""Root finding algorithm - others + +This module provides the divided difference backend method for root-finding algorithms. + +- div_diff + +Created on Jan 09 01:03:46 2024 +""" + +__all__ = ['div_diff'] + +from .. import FList, Func, IFloatOrFList + + +# Taken form https://en.wikipedia.org/wiki/Muller%27s_method#Computational_example +# minor tweaking applied on variable namings for consistency +def div_diff(function: Func, xs_: FList) -> IFloatOrFList: + """Calculate the divided difference f[x0, x1, ...].""" + if len(xs_) == 2: + a, b = xs_ + return (function(a) - function(b)) / (a - b) + else: + f1 = div_diff(function, xs_[1:]) + f2 = div_diff(function, xs_[0:-1]) + + return (f1 - f2) / (xs_[-1] - xs_[0]) diff --git a/src/Solvers/interpolators/_backend/spline_.py b/src/num_solvers/__backend/spline_.py similarity index 98% rename from src/Solvers/interpolators/_backend/spline_.py rename to src/num_solvers/__backend/spline_.py index 1ac86a3..4aac3a6 100644 --- a/src/Solvers/interpolators/_backend/spline_.py +++ b/src/num_solvers/__backend/spline_.py @@ -1,4 +1,4 @@ -"""SPLINE_ +"""Spline interpolation base module This module provides the basic classes to perform different types of spline interpolations, - SPLINE: The base class for spline interpolation. @@ -24,8 +24,8 @@ import numpy as np from custom_inherit import doc_inherit -from .ERRORS_ import AtLeastOneParameterRequired, WrongBoundaryEquation -from ... import DOC_STYLE, FList, FListOrLList, IFloat, IFloatOrFList, LLList, LList, OptFunc, OptList, TOLERANCE +from .. import DOC_STYLE, FList, FListOrLList, IFloat, IFloatOrFList, LLList, LList, OptFunc, OptList, TOLERANCE +from ..__backend.errors_ import AtLeastOneParameterRequired, WrongBoundaryEquation class SPLINE: diff --git a/src/num_solvers/__init__.py b/src/num_solvers/__init__.py new file mode 100644 index 0000000..7f52aa1 --- /dev/null +++ b/src/num_solvers/__init__.py @@ -0,0 +1,30 @@ +"""Created on Dec 20 14:10:28 2023""" + +from typing import Callable, Optional, Union + +from numpy.typing import NDArray + +TOLERANCE = 1e-10 +N_DECIMAL = 8 + +DOC_STYLE = 'numpy_napoleon_with_merge' + +Func = Callable +NdArray = NDArray +IFloat = Union[float, int] + +FList = list[IFloat] + +LList = list[FList] +FTuple = tuple[IFloat, ...] +LLList = list[LList] +LFunc = list[Func] + +FListOrLList = Union[FList, LList] +IFloatOrFList = Union[IFloat, FList] + +OptFunc = Optional[Func] +OptList = Optional[FList] +OptIFloat = Optional[IFloat] + +NdArray2 = tuple[NdArray, NdArray] diff --git a/src/num_solvers/integrators/__init__.py b/src/num_solvers/integrators/__init__.py new file mode 100644 index 0000000..b2f73d3 --- /dev/null +++ b/src/num_solvers/integrators/__init__.py @@ -0,0 +1 @@ +"""Created on Dec 20 22:59:44 2023""" diff --git a/src/num_solvers/integrators/newton_cotes_integrators.py b/src/num_solvers/integrators/newton_cotes_integrators.py new file mode 100644 index 0000000..8d67fe6 --- /dev/null +++ b/src/num_solvers/integrators/newton_cotes_integrators.py @@ -0,0 +1,85 @@ +"""Newton-Cotes integrators + +This module implements several Newton-Cotes integrators, such as, + +- trapezoidal_rule: + The trapezoidal rule integrator. +- simpson_13: + Simpson 1/3 rule integrator. +- simpson_38: + Simpson 3/8 rule integrator. +- boole: + Bool rule integrator. +- weddle_41h140: + Weddle rule integrator. +- weddle_3h10: + Modified Weddle rule integrator. + +Created on Nov 29 14:48:53 2023 +""" + +__all__ = ['trapezoid_rule', 'simpson_13', 'simpson_38', 'boole', 'weddle_41h140', 'weddle_3h10'] + +from custom_inherit import doc_inherit + +from .. import DOC_STYLE, Func, IFloat +from ..__backend.newton_cotes_ import solve_helper + + +def trapezoid_rule(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using the Trapezoidal Rule. + + Parameters + ---------- + function : + The function to be integrated. + x_0 : + The lower limit of integration. + x_n : + The upper limit of integration. + composites : + Number of sub-intervals for composite method. Default is 1. + get_table : + Whether to print the table of values. Default is False. + + Returns + ------- + The numerical result of integration. + """ + + return solve_helper(function, x_0, x_n, composites, get_table, 'trapezoid') + + +@doc_inherit(trapezoid_rule, style=DOC_STYLE) +def simpson_13(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using Simpson's 1/3 Rule.""" + + return solve_helper(function, x_0, x_n, composites, get_table, 'simpson_13') + + +@doc_inherit(trapezoid_rule, style=DOC_STYLE) +def simpson_38(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using Simpson's 3/8 Rule.""" + + return solve_helper(function, x_0, x_n, composites, get_table, 'simpson_38') + + +@doc_inherit(trapezoid_rule, style=DOC_STYLE) +def boole(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using Boole's Rule.""" + + return solve_helper(function, x_0, x_n, composites, get_table, 'boole') + + +@doc_inherit(trapezoid_rule, style=DOC_STYLE) +def weddle_41h140(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using Weddle 41/140 Rule.""" + + return solve_helper(function, x_0, x_n, composites, get_table, 'weddle_41h140') + + +@doc_inherit(trapezoid_rule, style=DOC_STYLE) +def weddle_3h10(function: Func, x_0: IFloat, x_n: IFloat, composites: int = 1, get_table: bool = False) -> IFloat: + """Perform numerical integration using modified Weddle 3/10 Rule.""" + + return solve_helper(function, x_0, x_n, composites, get_table, 'weddle_3h10') diff --git a/src/Solvers/interpolators/__init__.py b/src/num_solvers/interpolators/__init__.py similarity index 100% rename from src/Solvers/interpolators/__init__.py rename to src/num_solvers/interpolators/__init__.py diff --git a/src/num_solvers/interpolators/lagrange_interpolation_method.py b/src/num_solvers/interpolators/lagrange_interpolation_method.py new file mode 100644 index 0000000..0e52f3c --- /dev/null +++ b/src/num_solvers/interpolators/lagrange_interpolation_method.py @@ -0,0 +1,30 @@ +"""Lagrange's interpolation method + +This module provides the functionality to implement the Lagrange Interpolation method via, + +- lagrange_interpolation: A function to compute the interpolation value using Lagrange method. + +Created on Nov 01 16:47:03 2023 +""" + +__all__ = ['lagrange_interpolation'] + +from custom_inherit import doc_inherit + +from .. import DOC_STYLE, FList, IFloat, N_DECIMAL, OptFunc, OptList +from ..__backend.interpolation_ import Interpolation, LagrangeInterpolation + + +@doc_inherit(Interpolation.__init__, style=DOC_STYLE) +def lagrange_interpolation(given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None, n_decimal: int = N_DECIMAL): + """Performs interpolation via lagrange method. + + Returns + ------- + The interpolated result. + """ + + interpolator = LagrangeInterpolation(given_values, value_to_approximate, function, function_values, n_decimal) + + return interpolator.interpolate() diff --git a/src/num_solvers/interpolators/newton_interpolation_methods.py b/src/num_solvers/interpolators/newton_interpolation_methods.py new file mode 100644 index 0000000..659c074 --- /dev/null +++ b/src/num_solvers/interpolators/newton_interpolation_methods.py @@ -0,0 +1,49 @@ +"""Newton interpolation methods + +This module provides the functionality to implement three Newton interpolation methods, +- forwards_interpolation +- backwards_interpolation +- divided_interpolation + +Created on Nov 02 00:25:14 2023 +""" + +__all__ = ['forward_interpolation', 'backward_interpolation', 'divided_interpolation'] + +from custom_inherit import doc_inherit + +from .. import DOC_STYLE, FList, IFloat, N_DECIMAL, OptFunc, OptList +from ..__backend.interpolation_ import BkwInterpolation, DividedInterpolation, FwdInterpolation, Interpolation + + +@doc_inherit(Interpolation.__init__, style=DOC_STYLE) +def forward_interpolation(given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None, n_decimal: int = N_DECIMAL) -> IFloat: + """ + Perform forward difference interpolation. + + Returns + ------- + The interpolated result. + """ + + interpolation = FwdInterpolation(given_values, value_to_approximate, function, function_values, n_decimal) + return interpolation.interpolate() + + +@doc_inherit(forward_interpolation, style=DOC_STYLE) +def backward_interpolation(given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None): + """Perform backwards difference interpolation.""" + + interpolation = BkwInterpolation(given_values, value_to_approximate, function, function_values) + return interpolation.interpolate() + + +@doc_inherit(forward_interpolation, style=DOC_STYLE) +def divided_interpolation(given_values: FList, value_to_approximate: IFloat, function: OptFunc = None, + function_values: OptList = None): + """Perform divided difference interpolation.""" + + interpolation = DividedInterpolation(given_values, value_to_approximate, function, function_values) + return interpolation.interpolate() diff --git a/src/Solvers/interpolators/spline_interpolation.py b/src/num_solvers/interpolators/spline_interpolation.py similarity index 95% rename from src/Solvers/interpolators/spline_interpolation.py rename to src/num_solvers/interpolators/spline_interpolation.py index cf04f5a..ee99b6f 100644 --- a/src/Solvers/interpolators/spline_interpolation.py +++ b/src/num_solvers/interpolators/spline_interpolation.py @@ -1,12 +1,10 @@ -"""SplineInterpolation +"""Spline interpolation methods This module provides the interpolation methods for linear, quadratic and natural cubic spline, - linear_spline_interpolation - quadratic_spline_interpolation - natural_cubic_spline_interpolation -In almost all cases, this user should use this module for spline interpolation. - Created on Nov 01 23:03:42 2023 """ @@ -14,8 +12,8 @@ from custom_inherit import doc_inherit -from ._backend.spline_ import LinearSpline, NaturalCubicSpline, QuadraticSpline, SPLINE from .. import DOC_STYLE, FList, IFloat, OptFunc, OptList +from ..__backend.spline_ import LinearSpline, NaturalCubicSpline, QuadraticSpline, SPLINE @doc_inherit(SPLINE.__init__, style=DOC_STYLE) diff --git a/src/Solvers/iterative_solvers/__init__.py b/src/num_solvers/iterative_solvers/__init__.py similarity index 100% rename from src/Solvers/iterative_solvers/__init__.py rename to src/num_solvers/iterative_solvers/__init__.py diff --git a/src/num_solvers/iterative_solvers/function_root.py b/src/num_solvers/iterative_solvers/function_root.py new file mode 100644 index 0000000..d276ffd --- /dev/null +++ b/src/num_solvers/iterative_solvers/function_root.py @@ -0,0 +1,322 @@ +"""Root finding algorithms + +This module provides several functions for finding the roots of given uni-variable function. An example of such +function will be + +f(x) = x**2 - cos(x) + +Or + +f(y) = y**2 + 2*y - 1 + +For such uni-variable functions, the root finding methods include, + +- bisection_method +- false_position_method +- regula_falsi_method +- generalized_secant_method +- muller_method +- newton_raphson_method +- ridder_method +- secant_method +- sidi_method +- steffensen_method + +Created on Dec 24 12:56:39 2023 +""" + +__all__ = ['bisection_method', 'false_position_method', 'regula_falsi_method', 'secant_method', 'newton_raphson_solver', + 'generalized_secant_method', 'sidi_method', 'ridder_method', 'steffensen_method', 'muller_method'] + +from cmath import sqrt +from math import sqrt + +from custom_inherit import doc_inherit + +from .. import DOC_STYLE, Func, IFloat, IFloatOrFList, TOLERANCE +from ..__backend.root_finding_algorithms_ import div_diff + + +# TODO: modify newton_raphson method out of recursive function + + +def bisection_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_all_values: bool = False) -> IFloatOrFList: + """ + Use the bisection method to find the root of a given function within a specified interval. + + Parameters + ---------- + function : + The function for which the root is being sought. + x_start : + The starting point of the interval. + x_end : + The ending point of the interval. + tolerance : + Tolerance for convergence. Default is 1e-10. + get_all_values : + If True, returns a list of all intermediate values. If False, returns only the final root. + + Returns + ------- + The approximate root of the function or a list of all intermediate values. + + Notes + ----- + The function assumes that the initial interval [x_start, x_end] contains a root. + """ + + f, root = function, [] + + def is_sane(): + return f(x_start) * f(x_end) < 0 + + while is_sane(): + x_new = (x_start + x_end) / 2 + + root.append(x_new) + f_x_new = f(x_new) + + if abs(f_x_new) < tolerance: + break + + if f_x_new < 0: + x_start = x_new + else: + x_end = x_new + + return root if get_all_values else root[-1] + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def false_position_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_all_values: bool = False) -> IFloatOrFList: + """Use the false position method to find the root of a given function within a specified interval.""" + f, root = function, [] + + def is_sane(): + return f(x_start) * f(x_end) < 0 + + while is_sane(): + f_x0, f_x1 = f(x_start), f(x_end) + + x_new = (x_start * f_x1 - x_end * f_x0) / (f_x1 - f_x0) + root.append(x_new) + + f_x_new = f(x_new) + + if abs(f_x_new) < tolerance: + break + + if f_x_new == 0: + return root if get_all_values else root[-1] + + if f_x0 * f_x_new < 0: + x_end = x_new + else: + x_start = x_new + + return root if get_all_values else root[-1] + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def regula_falsi_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_all_values: bool = False) -> IFloatOrFList: + """Use the regula-falsi method to find the root of a given function within a specified interval.""" + return false_position_method(function, x_start, x_end, tolerance, get_all_values) + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def secant_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_all_values: bool = False) -> IFloatOrFList: + """Apply the secant method to find the root of a given function within a specified interval.""" + + f, root = function, [] + + while True: + x_new = x_end - f(x_end) * (x_end - x_start) / (f(x_end) - f(x_start)) + + root.append(x_new) + x_start, x_end = x_end, x_new + + f_x_new = f(x_new) + if abs(f_x_new) < tolerance: + break + + if f_x_new == 0: + return root if get_all_values else root[-1] + + return root if get_all_values else root[-1] + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def generalized_secant_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_full_result: bool = True) -> IFloatOrFList: + """Use the generalized secant method to find the root of a given function within a specified interval.""" + + x_2 = x_end - function(x_end) / div_diff(function, [x_start, x_end]) + + root_ = [x_start, x_end, x_2] + + while True: + f0 = function(root_[-1]) + + f1 = div_diff(function, [x_2, x_end]) + f2 = div_diff(function, [x_2, x_end, x_start]) * (x_2 - x_end) + + f_ = f1 + f2 + x_3 = x_2 - (f0 / f_) + + x_start, x_end, x_2 = x_end, x_2, x_3 + root_.append(x_3) + + if abs(f0) < tolerance: + break + + return root_ if get_full_result else root_[-1] + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def sidi_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_full_result: bool = True) -> IFloatOrFList: + """Use the Sidi method to find the root of a given function within a specified interval.""" + return generalized_secant_method(function, x_start, x_end, tolerance, get_full_result) + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def ridder_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_full_result: bool = True) -> IFloatOrFList: + """Use the Ridder method to find the root of a given function within a specified interval.""" + + def sign_function(functions_to_evaluate, value): + f1, f2 = functions_to_evaluate + return value if f1 - f2 > 0 else -value if f2 - f1 > 0 else 0 + + def is_sane(): + return True if f(x_start) * f(x_end) < 0 else False + + f, root_ = function, [x_start, x_end, (x_start + x_end) / 2] + + while is_sane(): + x_mid = (x_start + x_end) / 2 + + num_ = f(x_mid) * (x_mid - x_start) + den_ = sqrt(f(x_mid)**2 - f(x_start) * f(x_end)) + f_ = num_ / den_ + + x_new = x_mid + sign_function([f(x_start), f(x_end)], f_) + + if abs(f(x_new)) < tolerance or f_ == 0: + break + + if f(x_mid) * f(x_new) < 0: + x_start, x_end = x_mid, x_new + else: + x_start = x_new + + root_.append(x_new) + + return root_ if get_full_result else root_[-1] + + +@doc_inherit(bisection_method, style=DOC_STYLE) +def steffensen_method(function: Func, x_start: IFloat, x_end: IFloat, tolerance: IFloat = TOLERANCE, + get_full_result: bool = True) -> IFloatOrFList: + """Use the Steffensen method to find the root of a given function within a specified interval.""" + + def is_sane(): + return True if f(x_start) * f(x_end) < 0 else False + + f, root_ = function, [x_start, x_end, (x_start + x_end) / 2] + + solve = True if is_sane() else False + + while solve: + f_mid = f(root_[-1]) + f_mid2 = f(root_[-1] + f_mid) + + num_ = f_mid**2 + den_ = f_mid2 - f_mid + + x_new = root_[-1] - (num_ / den_) + + root_.append(x_new) + + if f(x_new) == 0 or abs(f(x_new)) < tolerance: + break + + return root_ if get_full_result else root_[-1] + + +def newton_raphson_method(function: Func, derivative_of_function: Func, initial_guess: IFloat, + tolerance: IFloat = TOLERANCE) -> IFloat: + """ + Find the root of a function using the Newton-Raphson method. + + Parameters + ---------- + function : callable + The function for which the root is being sought. + derivative_of_function : callable + The derivative of the target function. + initial_guess : float + The initial guess for the root. + tolerance : float, optional + Tolerance for convergence. Default is 1e-14. + + Returns + ------- + float + The approximate root of the function. + """ + f, df, x_0 = function, derivative_of_function, initial_guess + + if abs(f(x_0)) < tolerance: + return x_0 + else: + return newton_raphson_solver(f, df, x_0 - f(x_0) / df(x_0), tolerance) + + +def muller_method(function: Func, x_0: IFloat, x_1: IFloat, x_2: IFloat, iterations: int, + get_full_result: bool = True) -> IFloatOrFList: + """Apply the muller method to find the root of a given function within a specified interval. + + Parameters + ---------- + function: + The function for which the root is being sought. + x_0: + The first initial value. + x_1: + The second initial value. + x_2: + The third initial value. + iterations: + Number of iterations to perform. + get_full_result: + If True, returns a list of all intermediate values. If False, returns only the final root. + + Returns + ------- + The approximate root of the function or a list of all intermediate values. + + Notes + ----- + The function assumes that the initial interval [x_start, x_end] contains a root. + """ + + root_ = [x_0, x_1, x_2] + + for _ in range(iterations): + w = div_diff(function, [x_2, x_1]) + div_diff(function, [x_2, x_0]) - div_diff(function, [x_2, x_1]) + s_delta = sqrt(w**2 - 4 * function(x_2) * div_diff(function, [x_2, x_1, x_0])) + denominators = [w + s_delta, w - s_delta] + # Take the higher-magnitude denominator + x3 = x_2 - 2 * function(x_2) / max(denominators, key=abs) + # Advance + x_0, x_1, x_2 = x_1, x_2, x3 + + root_.append(x3) + + return root_ if get_full_result else root_[-1] diff --git a/src/Solvers/iterative_solvers/odes/__init__.py b/src/num_solvers/iterative_solvers/odes/__init__.py similarity index 100% rename from src/Solvers/iterative_solvers/odes/__init__.py rename to src/num_solvers/iterative_solvers/odes/__init__.py diff --git a/src/num_solvers/iterative_solvers/odes/runge_kutta.py b/src/num_solvers/iterative_solvers/odes/runge_kutta.py new file mode 100644 index 0000000..ddd0921 --- /dev/null +++ b/src/num_solvers/iterative_solvers/odes/runge_kutta.py @@ -0,0 +1,109 @@ +"""Runge Kutta Methods for single ODEs + +This module implements runge-kutta methods for single ODEs. The RK methods implemented are, + +- rk2_solver: + RK2 method for odes +- rk3_solver: + RK3 method for odes +- rk4_solver: + RK4 method for odes + +Created on Dec 16 03:28:44 2023 +""" + +__all__ = ['rk2_solver', 'rk3_solver', 'rk4_solver'] + +from custom_inherit import doc_inherit + +from ... import DOC_STYLE, Func, IFloat, LList, N_DECIMAL +from ...__backend.extra_ import round_value_ + + +def rk2_solver(ode: Func, x_initial: IFloat, y_initial: IFloat, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = N_DECIMAL) -> LList: + """ + Solve a first-order ordinary differential equation using the Runge-Kutta (RK2) method. + + Parameters + ---------- + ode: + The ordinary differential equation in the form dy/dx = ode(x, y). + x_initial: + The x value in the initial condition of the ODE. + y_initial: + The y value in the initial condition of the ODE. + step_size: + Step size for numerical integration, default is 0.1. + x_max: + Maximum x-value for integration, default is 1.0. + n_decimal: + Number of decimal places to round up to. Default is 8. + + Returns + ------- + Tuple containing arrays of x and y values at each step. + """ + + num_steps = int((x_max - x_initial) / step_size) + 1 + + x_values, y_values = [0] * num_steps, [0] * num_steps + x_values[0], y_values[0] = x_initial, y_initial + + for i in range(1, num_steps): + x_i, y_i = x_values[i - 1], y_values[i - 1] + + k1 = step_size * ode(x_i, y_i) + k2 = step_size * ode(x_i + step_size, y_i + k1) + + y_values[i] = round_value_(y_i + 0.5 * (k1 + k2), n_decimal) + x_values[i] = round_value_((x_i + step_size), n_decimal) + + return [x_values, y_values] + + +@doc_inherit(rk2_solver, style=DOC_STYLE) +def rk3_solver(ode: Func, x_initial: IFloat, y_initial: IFloat, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = N_DECIMAL) -> LList: + """Solve a first-order ordinary differential equation using the Runge-Kutta (RK3) method.""" + + num_steps = int((x_max - x_initial) / step_size) + 1 + + x_values, y_values = [0] * num_steps, [0] * num_steps + x_values[0], y_values[0] = x_initial, y_initial + + for i in range(1, num_steps): + x_i, y_i = x_values[i - 1], y_values[i - 1] + + k1 = step_size * ode(x_i, y_i) + k2 = step_size * ode(x_i + step_size / 2, y_i + k1 / 2) + k3 = step_size * ode(x_i + step_size, y_i - k1 + 2 * k2) + + y_values[i] = round_value_(y_i + (1 / 6) * (k1 + 4 * k2 + k3), n_decimal) + x_values[i] = round_value_(x_i + step_size, n_decimal) + + return [x_values, y_values] + + +@doc_inherit(rk2_solver, style=DOC_STYLE) +def rk4_solver(ode: Func, x_initial: IFloat, y_initial: IFloat, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = N_DECIMAL) -> LList: + """Solve a first-order ordinary differential equation using the Runge-Kutta (RK4) method.""" + + num_steps = int((x_max - x_initial) / step_size) + 1 + + x_values, y_values = [0] * num_steps, [0] * num_steps + x_values[0], y_values[0] = x_initial, y_initial + + for i in range(1, num_steps): + x_i, y_i = x_values[i - 1], y_values[i - 1] + + k1 = step_size * ode(x_i, y_i) + k2 = step_size * ode(x_i + step_size / 2, y_i + k1 / 2) + k3 = step_size * ode(x_i + step_size / 2, y_i + k2 / 2) + k4 = step_size * ode(x_i + step_size, y_i + k3) + + y_values[i] = round_value_(y_i + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4), n_decimal) + x_values[i] = round_value_(x_i + step_size, n_decimal) + + return [x_values, y_values] diff --git a/src/num_solvers/iterative_solvers/odes/runge_kutta_multi.py b/src/num_solvers/iterative_solvers/odes/runge_kutta_multi.py new file mode 100644 index 0000000..24a868b --- /dev/null +++ b/src/num_solvers/iterative_solvers/odes/runge_kutta_multi.py @@ -0,0 +1,219 @@ +"""Runge Kutta Methods for Multi-ODE Systems + +This module implements runge-kutta methods for Multi-ODE systems. The RK methods implemented are, + +- rk2_multi_ode: + RK2 method for multi-odes +- rk3_multi_ode: + RK3 method for multi-odes +- rk4_multi_ode: + RK4 method for multi-odes + +The module also generalizes a few RK operations, such as, + +- _k_0: For RK2, Rk3, and RK4 + To compute the k0 factor for RK methods. +- _k_n: For RK2, and RK4 + To compute the kn factor for RK methods. +- _multi_rk_result: For RK2, Rk3, and RK4 + Method to implement various final results of RK method depending upon the method used. + +Created on Jan 05 21:54:09 2024 +""" + +__all__ = ['rk2_multi_ode', 'rk3_multi_ode', 'rk4_multi_ode'] + +import numpy as np +from custom_inherit import doc_inherit + +from ... import DOC_STYLE, FList, IFloat, LFunc, LList, N_DECIMAL, NdArray + + +def rk2_multi_ode(odes: LFunc, initial_conditions: FList, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = N_DECIMAL) -> LList: + """ + Solve a system of first-order ordinary differential equations using the Runge-Kutta (RK2) method. + + Parameters + ---------- + odes : List[Callable] + List of ordinary differential equations in the form dy_i/dx = odes[i](x, y_1, y_2, ..., y_n). + initial_conditions : List[float] + List of initial values for each quantity (x_0, y_1, y_2, ..., y_n). + step_size : float, optional + Step size for numerical integration, default is 0.1. + x_max : float, optional + Maximum x-value for integration, default is 1.0. + n_decimal : int, optional + Number of decimal places to round the result, default is 9. + + Returns + ------- + List of lists containing the values of each quantity at each step. + """ + + num_odes = len(odes) + 1 + num_steps = int((x_max - initial_conditions[0]) / step_size) + 1 + + result = np.zeros((num_steps, num_odes)) + result[0] = initial_conditions + + for i in range(1, num_steps): + quantities = result[i - 1] + + temp_ = _k_0(num_odes, odes, quantities, step_size, 2) + temp_ = _k_n(num_odes, odes, quantities, step_size, temp_, 0) + + result = _multi_rk_result(i, temp_, quantities, step_size, num_odes, result, 2, n_decimal) + + return [result[:, i].tolist() for i in range(num_odes)] + + +@doc_inherit(rk2_multi_ode, style=DOC_STYLE) +def rk3_multi_ode(odes: LFunc, initial_conditions: FList, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = 9) -> LList: + """Solve a system of first-order ordinary differential equations using the Runge-Kutta (RK3) method.""" + + num_odes = len(odes) + 1 + num_steps = int((x_max - initial_conditions[0]) / step_size) + 1 + + result = np.zeros((num_steps, num_odes)) + result[0] = initial_conditions + + for i in range(1, num_steps): + quantities = result[i - 1] + + temp_ = _k_0(num_odes, odes, quantities, step_size, 3) + + temp__ = [quantities[0] + step_size / 2] + list(quantities[1:] + temp_[0, :] / 2) + temp_[1, :] = step_size * np.array([odes[j](*temp__) for j in range(num_odes - 1)]) + + temp__ = [quantities[0] + step_size] + list(quantities[1:] - temp_[0, :] + 2 * temp_[1, :]) + temp_[2, :] = step_size * np.array([odes[j](*temp__) for j in range(num_odes - 1)]) + + result = _multi_rk_result(i, temp_, quantities, step_size, num_odes, result, 3, n_decimal) + + return [result[:, i].tolist() for i in range(num_odes)] + + +@doc_inherit(rk2_multi_ode, style=DOC_STYLE) +def rk4_multi_ode(odes: LFunc, initial_conditions: FList, step_size: IFloat = 0.1, x_max: IFloat = 1.0, + n_decimal: IFloat = 9) -> LList: + """Solve a system of first-order ordinary differential equations using the Runge-Kutta (RK4) method.""" + + num_odes = len(odes) + 1 + num_steps = int((x_max - initial_conditions[0]) / step_size) + 1 + + result = np.zeros((num_steps, num_odes)) + result[0] = initial_conditions + + for i in range(1, num_steps): + quantities = result[i - 1] + + temp_ = _k_0(num_odes, odes, quantities, step_size, 4) + + for k in range(1, 3): + temp__ = [quantities[0] + step_size / 2] + list(quantities[1:] + temp_[k - 1, :] / 2) + temp_[k, :] = step_size * np.array([odes[j](*temp__) for j in range(num_odes - 1)]) + + temp_ = _k_n(num_odes, odes, quantities, step_size, temp_, 2) + + result = _multi_rk_result(i, temp_, quantities, step_size, num_odes, result, 4, n_decimal) + + return [result[:, i].tolist() for i in range(num_odes)] + + +def _k_0(num_odes: int, odes: LFunc, quantities: NdArray, step_size: IFloat, rk_type: int) -> NdArray: + """ + Computes the k0 parameter for RK methods. + + Parameters + ---------- + num_odes: + Number of ODEs. + odes: + A list containing the ODEs. + quantities: + The values used for computing the k0 factor + step_size: + The step size for Rk method + rk_type: + The type of RK method used, 2 for RK2, 3 for RK3, and 4 for RK4 method. + + Returns + ------- + k0 parameter of RK method. + """ + + k0_ = np.zeros((rk_type, num_odes - 1)) + + for j in range(num_odes - 1): + k0_[0][j] = step_size * odes[j](*quantities) + + return k0_ + + +@doc_inherit(_k_0, style=DOC_STYLE) +def _k_n(num_odes: int, odes: LFunc, quantities: NdArray, step_size: IFloat, temp_: NdArray, index: int) -> NdArray: + """ + Compute the kn parameter for RK methods. + + Parameters + ---------- + temp_: + The temporary array holding the data + index: + The index of the values to pick. + + Returns + ------- + Modified temporary array. + """ + temp__ = [quantities[0] + step_size] + list(quantities[1:] + temp_[index, :]) + temp_[index + 1, :] = step_size * np.array([odes[j](*temp__) for j in range(num_odes - 1)]) + + return temp_ + + +def _multi_rk_result(i: int, temp_: NdArray, quantities: NdArray, step_size: IFloat, num_odes: int, result: NdArray, + rk_type: int, n_decimal: int) -> NdArray: + """ + Returns the RK method results. + + Parameters + ---------- + i: + The iterative index for the RK method. + temp_: + The temporary array holding the processed values. + quantities: + The array holding the base values to be worked on. + step_size: + The step size for RK method. + num_odes: + The number of ODEs. + result: + The array holding the result values. + rk_type: + The type of RK method used, 2 for RK2, 3 for RK3, and 4 for RK4 method. + n_decimal: + The number of digits to round off to. + + Returns + ------- + The result of RK method calculations. + """ + for j in range(num_odes): + if j == 0: + result[i, j] = np.round(quantities[j] + step_size, n_decimal) + else: + temp__ = temp_[:, j - 1] + if rk_type == 2: + result[i, j] = quantities[j] + (1 / 2) * (temp__[0] + temp__[1]) + elif rk_type == 3: + result[i, j] = quantities[j] + (1 / 6) * (temp__[0] + 4 * temp__[1] + temp__[2]) + else: + result[i, j] = quantities[j] + (1 / 6) * (temp__[0] + 2 * np.sum(temp__[1:3]) + temp__[3]) + result[i, j] = np.round(result[i, j], n_decimal) + + return result diff --git a/src/num_solvers/iterative_solvers/polynomial_root.py b/src/num_solvers/iterative_solvers/polynomial_root.py new file mode 100644 index 0000000..39b0675 --- /dev/null +++ b/src/num_solvers/iterative_solvers/polynomial_root.py @@ -0,0 +1,163 @@ +"""Polynomial root finding algorithms + +This module provides functionality to find roots of given polynomials. The roots are determined via, + +- laguerre_method + +In addition, the user can specify a range in which all roots are to be determined for the given polynomial via, + +- segmented_roots + +This module also provides a handy function to determine random polynomials, via + +- generate_random_polynomial + +Created on Dec 20 23:21:44 2023 +""" + +from cmath import sqrt +from typing import List + +import numpy as np + +from .. import FList, IFloat, LList, OptList, TOLERANCE +from ..__backend.errors_ import DegreeOfPolynomialNotCorrect +from ..__backend.extra_ import round_list_ + + +# TODO: See if numpy can be removed + + +def laguerre_method(polynomial: FList, degree_of_polynomial: int = -1, x_guess: IFloat = 0, + get_full_result: bool = False, tolerance: IFloat = TOLERANCE): + """ + Use Laguerre method to find the roots of the given polynomial. + + Parameters + ---------- + polynomial: + The polynomial as list of coefficients. + degree_of_polynomial: + Degree of the polynomial. Defaults to -1. If polynomial's degree is -1, it is calculated inside the function. + x_guess: + Initial guess for the root. Defaults to 0. + get_full_result: + If True, gives all the guesses. If false, only gives the root of polynomial. + tolerance: + Tolerance for the result. Default is TOLERANCE. + + Returns + ------- + Root of the given polynomial. + """ + + p_val, p_der = np.polyval, np.polyder + root_ = [x_guess] + + while True: + poly = p_val(polynomial, x_guess) + + if abs(poly) < tolerance: + break + + poly_derivative = p_val(p_der(polynomial), x_guess) + poly_double_derivative = p_val(p_der(polynomial, 2), x_guess) + + g_ = poly_derivative / poly + h_ = g_**2 - (poly_double_derivative / poly) + + sqrt_in_denominator = sqrt((degree_of_polynomial - 1) * (degree_of_polynomial * h_ - g_**2)) + p1 = max([g_ + sqrt_in_denominator, g_ - sqrt_in_denominator], key=abs) + + x_guess = x_guess - (degree_of_polynomial / p1) + + root_.append(x_guess) + + return root_ if get_full_result else root_[-1] + + +def segmented_roots(polynomial: List, degree_of_polynomial: int = -1, x_guess: OptList = None, + num_segments: int = 500, n_decimal: int = 8, tolerance: IFloat = TOLERANCE) -> LList: + """ + Segments the given interval to find all possible roots within that interval for a given polynomial. + + Parameters + ---------- + polynomial: + The polynomial as list of coefficients. + degree_of_polynomial: + Degree of the polynomial. Defaults to -1. If polynomial's degree is -1, it is calculated inside the function. + x_guess: + The range in which the root of polynomial is to be determined. Default is [-100, 100]. + num_segments: + Number of segments to divide the ``x_range`` into. Default is 500. + n_decimal: + Number of digits to round off to. Default is 8 + tolerance: + Tolerance for the result. Default is TOLERANCE. + + Returns + ------- + List of all possible roots within the given interval. + """ + + x_guess = x_guess if x_guess else [-100, 100] + degree_of_polynomial = degree_of_polynomial if degree_of_polynomial > 0 else len(polynomial) - 1 + + if degree_of_polynomial: + if len(polynomial) - 1 != degree_of_polynomial: + raise DegreeOfPolynomialNotCorrect('The provided polynomial and degree of polynomial do not match') + + x_values = np.linspace(x_guess[0], x_guess[1], num_segments) + all_roots = [] + + print(f'Checking roots between {x_guess[0]:+} and {x_guess[1]:+} using {num_segments} intervals.') + + for x_0 in x_values: + root = laguerre_method(polynomial, degree_of_polynomial, x_0, tolerance=tolerance) + + if all(abs(root - existing_root) > tolerance for existing_root in all_roots): + all_roots.append(root) + + unique_roots = np.fromiter(all_roots, dtype=complex) + + new_roots = [] + for i, root_ in enumerate(unique_roots): + if np.allclose(np.imag(root_), 0): + new_roots.append(np.real(root_)) + else: + new_roots.extend([root_, np.conj(root_)]) + + new_roots = round_list_(new_roots, n_decimal) + + for i, v in enumerate(new_roots): + print(f'root_{i + 1:03} = {v}') + + return new_roots + + +def generate_random_polynomial(degree: int, low: int = -10, high: int = 10) -> FList: + """ + Generate random coefficients for a polynomial. + + Parameters + ---------- + degree : int + The degree of the polynomial. + low : int, optional + The lower bound for the random coefficients. Default is -10. + high : int, optional + The upper bound for the random coefficients. Default is 10. + + Returns + ------- + List[float] + A list of coefficients representing the random polynomial. + """ + + coefficients = np.random.uniform(low, high, degree + 1) + + while coefficients[0] == 0: + coefficients[0] = np.random.uniform(low, high) + + return list(coefficients) diff --git a/src/Solvers/iterative_solvers/system_of_equations/__init__.py b/src/num_solvers/iterative_solvers/system_of_equations/__init__.py similarity index 100% rename from src/Solvers/iterative_solvers/system_of_equations/__init__.py rename to src/num_solvers/iterative_solvers/system_of_equations/__init__.py diff --git a/src/num_solvers/iterative_solvers/system_of_equations/gauss_jacobi.py b/src/num_solvers/iterative_solvers/system_of_equations/gauss_jacobi.py new file mode 100644 index 0000000..e352417 --- /dev/null +++ b/src/num_solvers/iterative_solvers/system_of_equations/gauss_jacobi.py @@ -0,0 +1,53 @@ +"""Gauss-Jacobi + +This module provides functionality to solve diagonally dominant system of equations using Gauss-Jacobi method. + +Created on Oct 13 11:37:27 2023 +""" + +__all__ = ['gauss_jacobi'] + +from typing import List, Tuple, Union + +import numpy as np +from custom_inherit import doc_inherit + +from ... import DOC_STYLE, FList, FTuple, IFloat, LLList, LList, TOLERANCE +from ...__backend.iterative_solvers_ import GaussJacobi, SysEqnSolver + +GJ_OUTPUT = Union[np.ndarray, Tuple[np.ndarray, List[np.ndarray]]] +INITIAL_GUESS = Tuple[float, float, float] + + +@doc_inherit(SysEqnSolver.__init__, style=DOC_STYLE) +def gauss_jacobi(system_of_equations: LList, + solution: FList, + n_iter: int = 500, + tolerance: IFloat = TOLERANCE, + initial_guess: FTuple = (0, 0, 0), + n_decimal: int = 8, + get_solution_set: bool = False) -> LLList: + """ + Solve a system of linear equations using the Gauss-Jacobi method. + + Parameters + ---------- + + get_solution_set: + If True, returns both the solution and the solution set for each iteration, default is False. + + Returns + ------- + If ``get_solution_set`` is ``False``, returns the solution as a list. + Otherwise, returns a tuple containing the solution and a list of solution sets for each iteration. + + Examples + -------- + >>> equations = [[2, 1, 1], [3, 10, 2], [2, 1, 4]] + >>> solutions = [5, 10, 9] + >>> GJ_Solver = gauss_jacobi(equations, solutions) + """ + + gj = GaussJacobi(system_of_equations, solution, n_iter, tolerance, initial_guess, n_decimal) + + return gj.solve() if not get_solution_set else (gj.solve(), gj.solution_set) diff --git a/src/num_solvers/iterative_solvers/system_of_equations/gauss_seidel.py b/src/num_solvers/iterative_solvers/system_of_equations/gauss_seidel.py new file mode 100644 index 0000000..b4c9216 --- /dev/null +++ b/src/num_solvers/iterative_solvers/system_of_equations/gauss_seidel.py @@ -0,0 +1,47 @@ +"""Gauss-Seidel + +This module provides functionality to solve diagonally dominant system of equations using Gauss-Seidel method. + +Created on Oct 14 05:46:04 2023 +""" + +__all__ = ['gauss_jacobi'] + +from custom_inherit import doc_inherit + +from ... import DOC_STYLE, FList, FTuple, IFloat, LLList, LList, TOLERANCE +from ...__backend.iterative_solvers_ import GaussSeidel, SysEqnSolver + + +@doc_inherit(SysEqnSolver.__init__, style=DOC_STYLE) +def gauss_seidel(system_of_equations: LList, + solution: FList, + n_iter: int = 500, + tolerance: IFloat = TOLERANCE, + initial_guess: FTuple = (0, 0, 0), + n_decimal: int = 8, + get_solution_set: bool = False) -> LLList: + """ + Solve a system of linear equations using the Gauss-Seidel method. + + Parameters + ---------- + + get_solution_set: + If True, returns both the solution and the solution set for each iteration, default is False. + + Returns + ------- + If ``get_solution_set`` is ``False``, returns the solution as a list. + Otherwise, returns a tuple containing the solution and a list of solution sets for each iteration. + + Examples + -------- + >>> equations = [[2, 1, 1], [3, 10, 2], [2, 1, 4]] + >>> solutions = [5, 10, 9] + >>> GS_Solver = gauss_seidel(equations, solutions) + """ + + gs = GaussSeidel(system_of_equations, solution, n_iter, tolerance, initial_guess, n_decimal) + + return gs.solve() if not get_solution_set else (gs.solve(), gs.solution_set) diff --git a/src/num_solvers/iterative_solvers/system_of_equations/predictor_corrector_methods.py b/src/num_solvers/iterative_solvers/system_of_equations/predictor_corrector_methods.py new file mode 100644 index 0000000..66b7f9a --- /dev/null +++ b/src/num_solvers/iterative_solvers/system_of_equations/predictor_corrector_methods.py @@ -0,0 +1,130 @@ +"""Predictor-Corrector Methods. + +This module provides functionality to get iterative solutions using predictor-corrector method. Currently, this module +provides, + +- euler_trapezoidal: To solve first-order ODE via predictor-corrector method. + +In addition, this module also provides, separately, + +- euler_method: To solve a first-order ODE without predictor-corrector method. + +Created on Dec 20 10:55:46 2023 +""" + +from custom_inherit import doc_inherit + +from ... import DOC_STYLE, FList, Func, IFloat, LList, N_DECIMAL, OptIFloat +from ...__backend.errors_ import XFindNotDefined +from ...__backend.extra_ import linear_list, round_value_ + + +def trapezoidal_method(ode: Func, x_current: IFloat, y_current: IFloat, x_next: IFloat, y_next: IFloat, + step_size: IFloat, n_decimal: int = N_DECIMAL) -> IFloat: + """ + Calculate the next value using the Trapezoidal Rule in numerical integration. + + Parameters + ---------- + ode: + The ordinary differential equation function. + x_current: + Current x value. + y_current: + Current y value. + x_next: + Next x value. + y_next: + Next y value. + step_size: + Step size. + n_decimal: + Number of decimals to round the result to (default is N_DECIMAL). + + Returns + ------- + The next y value calculated using the Trapezoidal Rule. + """ + + def _trapezoidal(): + f_n = ode(x_current, y_current) + f_p = ode(x_next, y_next) + return y_current + (step_size / 2) * (f_n + f_p) + + return round_value_(_trapezoidal(), n_decimal) + + +def euler_method(ode: Func, x_initial: IFloat, y_initial: IFloat, step_size: IFloat = 0.1, x_find: OptIFloat = None, + n_decimal: int = 6, full_result: bool = False) -> LList or FList: + """ + Solve a first-order ordinary differential equation using the Euler method. + + Parameters + ---------- + ode : Callable + The ordinary differential equation in the form dy/dx = ode(x, y). + x_initial : float + Initial x-value. + y_initial : float + Initial y-value. + step_size : float, optional + Step size for numerical integration. If not provided, it is calculated based on num_steps and x_max. + x_find : float, optional + The value of x at which the value of function is to be determined. + n_decimal : int, optional + Number of decimal places to round the result, default is 6. + full_result : bool, optional + If True, return a tuple of arrays (x_values, y_values). If False, return the final y-value only. + + Returns + ------- + Result of the Euler method. If full_result is True, returns a tuple of arrays (x_values, y_values). + If full_result is False, returns the final y-value only. + """ + + if x_find is None: + raise XFindNotDefined('The value for x to find is not defined.') + + num_steps = int((x_find - x_initial) / step_size) + 1 + + x_values, y_values = [x_initial] * num_steps, [y_initial] * num_steps + + for i in range(1, num_steps): + y_values_ = y_values[i - 1] + step_size * ode(x_values[i - 1], y_values[i - 1]) + x_values_ = x_values[i - 1] + step_size + + y_values[i] = round_value_(y_values_, n_decimal) + x_values[i] = round_value_(x_values_, n_decimal) + + x_initial += step_size + + return [x_values, y_values] if full_result else y_values[-1] + + +@doc_inherit(euler_method, style=DOC_STYLE) +def euler_trapezoidal(ode: Func, x_initial: IFloat, y_initial: IFloat, step_size: IFloat = 0.1, + x_find: OptIFloat = None, n_decimal: int = 6, full_result: bool = False): + """Solve a first-order ordinary differential equation using the Euler method with trapezoidal correction.""" + + if x_find is None: + raise XFindNotDefined('The value for x to find is not defined.') + + num_steps = int((x_find - x_initial) / step_size) + 1 + + x_values = linear_list(x_initial, x_initial + (num_steps * step_size), num_steps + 1, n_decimal=n_decimal) + + y_values = [0] * len(x_values) + + predictor_values, corrector_values = [0] * len(x_values), [0] * len(x_values) + y_values[0], predictor_values[0], corrector_values[0] = [y_initial] * 3 + + for i in range(1, num_steps + 1): + predictor = euler_method(ode, x_values[i - 1], corrector_values[i - 1], step_size, n_decimal) + predictor_values[i] = round_value_(predictor, n_decimal) + + corrector = trapezoidal_method(ode, x_values[i - 1], y_values[i - 1], x_values[i], predictor_values[i], + step_size, n_decimal) + corrector_values[i] = corrector + y_values[i] = corrector + + return (x_values, (predictor_values, corrector_values)) if full_result else (x_values, corrector_values) diff --git a/src/num_solvers/matrix_decomposition/__init__.py b/src/num_solvers/matrix_decomposition/__init__.py new file mode 100644 index 0000000..7e67b65 --- /dev/null +++ b/src/num_solvers/matrix_decomposition/__init__.py @@ -0,0 +1 @@ +"""Created on Jan 09 05:07:28 2024""" diff --git a/src/num_solvers/matrix_decomposition/cholesky_decomposition.py b/src/num_solvers/matrix_decomposition/cholesky_decomposition.py new file mode 100644 index 0000000..a154cfa --- /dev/null +++ b/src/num_solvers/matrix_decomposition/cholesky_decomposition.py @@ -0,0 +1,121 @@ +"""Cholesky decomposition + +This module provides the functionality to decompose a symmetric and positive definite matrix A, using +Cholesky decomposition. + +Created on Jan 09 02:03:16 2024 +""" + +from math import sqrt + +from .. import LLList, LList +from ..__backend.errors_ import NonSymmetricMatrix, NotPositiveDefinite + +__all__ = ['cholesky_decomposition'] + + +def cholesky_decomposition(matrix_a: LList, n_decimal: int = 4) -> LLList: + """ + Performs the Cholesky decomposition on the given matrix. + + Parameters + ---------- + matrix_a: + The matrix to decompose. + n_decimal: + The number of digits to round off for the result. + + Returns + ------- + The Cholesky decomposition of the matrix. Matrix L, and L_star. + + Raises + ------ + NonSymmetricMatrix: + If the matrix is not symmetric. + NotPositiveDefinite: + If the matrix is not positive definite. + """ + + if not __is_symmetric(matrix_a): + raise NonSymmetricMatrix('The matrix is not symmetric. Can not perform Cholesky decomposition.') + + if not __is_positive_definite(matrix_a): + raise NotPositiveDefinite('The matrix is not positive definite. Can not perform Cholesky decomposition.') + + matrix_l = [[0] * len(matrix_a) for _ in range(len(matrix_a))] + + n_dimensions = len(matrix_a) + + matrix_l[0][0] = sqrt(matrix_a[0][0]) + for i in range(1, n_dimensions): + for k in range(n_dimensions): + if k == 0: + matrix_l[i][k] = matrix_a[k][i] / matrix_l[0][0] + elif k < i: + f1 = matrix_a[i][k] + f1_ = [matrix_l[i][j] * matrix_l[k][j] for j in range(i - 1)] + matrix_l[i][k] = (f1 - sum(f1_)) / matrix_l[k][k] + elif k == i: + f4 = sum((i**2 for i in matrix_l[i])) + matrix_l[i][k] = sqrt(matrix_a[i][i] - f4) + else: + matrix_l[i][k] = 0 + + matrix_l = [[round(element, n_decimal) for element in row] for row in matrix_l] + matrix_l_star = [[matrix_l[row][col] for row in range(n_dimensions)] for col in range(n_dimensions)] + + return [matrix_l, matrix_l_star] + + +def __is_symmetric(matrix: LList) -> bool: + """Checks if the given matrix is symmetric or not. + + Parameters + ---------- + matrix: + The matrix to check. + + Returns + ------- + True if the matrix is symmetric, else False. + """ + + n_dim = len(matrix) + matrix_transpose = [[matrix[row][col] for row in range(n_dim)] for col in range(n_dim)] + + return matrix == matrix_transpose + + +def __is_positive_definite(matrix: LList) -> bool: + """Checks if the matrix is positive definite or not. + + Parameters + ---------- + matrix: + The matrix to check. + + Returns + ------- + True if the matrix is positive definite, else False. + + """ + + n_dimensions = len(matrix) + + # Check if the matrix is square + if len(matrix) != len(matrix[0]): + return False + + # Check if the matrix is symmetric + if not __is_symmetric(matrix): + return False + + # Check if all leading principal minors have positive determinants + for i in range(1, n_dimensions + 1): + minor = [row[:i] for row in matrix[:i]] + determinant = sum((minor[j][j] for j in range(i))) + if determinant <= 0: + return False + + return True