diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 9012915..b583b36 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -3,9 +3,16 @@ Change Log .. currentmodule:: idesolver +v1.2.0 +------ + +* :math:`d(x)` and :math:`k(x, s)` can now be matrices or vectors (PR :pr:`60`, thanks `nbrucy `_ for reviewing). +* `UnexpectedlyComplexValuedIDE` is no longer raised if any `TypeError` occurs during integration. + This may lead to less-explicit errors, but this is preferable to incorrect error reporting. + v1.1.0 ------ -* Add support for multidimensional IDEs (PR :pr:`35` resolves :issue:`28`, thanks `nbrucy `_!) +* Add support for multidimensional IDEs (PR :pr:`35` resolves :issue:`28`, thanks `nbrucy `_!). v1.0.5 ------ diff --git a/docs/source/index.rst b/docs/source/index.rst index 2cfe03f..28f2266 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -16,12 +16,15 @@ The algorithm is based on a scheme devised by `Gelmi and Jorquera `_. +IDESolver is open-source software developed on `GitHub `_. +If you encounter problems or would like to extend it for new use cases, please consider opening an issue or a pull request. + :doc:`quickstart` A brief tutorial in using IDESolver. :doc:`manual` - Details about the implementation of IDESolver. - Includes information about running the test suite. + Details about the implementation of IDESolver, extensions like solving multidimensional IDEs, and subtleties like stopping conditions. + Also includes information about running the test suite. :doc:`api` Detailed documentation for IDESolver's API. diff --git a/docs/source/manual.rst b/docs/source/manual.rst index 847bdf7..4491e0d 100644 --- a/docs/source/manual.rst +++ b/docs/source/manual.rst @@ -50,6 +50,85 @@ To be conservative and to make sure we don't over-correct, we'll combine :math:` The process then repeats: solve the IDE-turned-ODE with :math:`y^{(1)}` on the right-hand-side, see how different it is, maybe make a new guess, etc. +Multidimensional IDEs +--------------------- + +IDESolver can handle multidimensional IDEs, though I do not know how well the convergence properties hold. + +There are several ways to specify multidimensional IDEs, +with IDESolver flexibly interpreting the arguments to allow for IDEs to be expressed in as simple a notation as possible +at each level of complexity. +For example, the three definitions below for a two-dimensional IDE +(where :math:`\gamma_0(x)` and :math:`\gamma_1(x)` are the functions to be solved for) +are all treated as equivalent by IDESolver. + +Scalar :math:`d` and :math:`k`: + +.. code-block:: + + IDESolver( + x=np.linspace(0, 7, 100), + y_0=[0, 1], + c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]], + d=lambda x: -0.5, + f=lambda y: y, + lower_bound=lambda x: 0, + upper_bound=lambda x: x, + ) + +which could be written out as + +.. math:: + + \frac{d}{dx}\begin{bmatrix} \gamma_0 \\ \gamma_1 \end{bmatrix} & = \begin{bmatrix} \frac{1}{2} (\gamma_1 + 1) \\ -\frac{1}{2} \gamma_0 \end{bmatrix} - \frac{1}{2} \int_{0}^{x} \, \begin{bmatrix} \gamma_0(s) \\ \gamma_1(s) \end{bmatrix} \, ds + +Vector :math:`d` and :math:`k` (treated as a scalar multiplier for each dimension): + + +.. code-block:: + + IDESolver( + x=np.linspace(0, 7, 100), + y_0=[0, 1], + c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]], + d=lambda x: [-0.5, -0.5], + k=lambda x, s: [1, 1], + f=lambda y: y, + lower_bound=lambda x: 0, + upper_bound=lambda x: x, + ) + +which could be written out as + +.. math:: + + \frac{d}{dx}\begin{bmatrix} \gamma_0 \\ \gamma_1 \end{bmatrix} & = \begin{bmatrix} \frac{1}{2} (\gamma_1 + 1) \\ -\frac{1}{2} \gamma_0 \end{bmatrix} - \begin{bmatrix} \frac{1}{2} \\ \frac{1}{2} \end{bmatrix} \int_{0}^{x} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \, \begin{bmatrix} \gamma_0(s) \\ \gamma_1(s) \end{bmatrix} \, ds + + +Matrix :math:`d` and :math:`k` (treated using matrix multiplication): + + +.. code-block:: + + IDESolver( + x=np.linspace(0, 7, 100), + y_0=[0, 1], + c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]], + d=lambda x: [[-0.5, 0], [0, -0.5]], + k=lambda x, s: [[1, 0], [0, 1]], + f=lambda y: y, + lower_bound=lambda x: 0, + upper_bound=lambda x: x, + ) + +which could be written out as + +.. math:: + + \frac{d}{dx}\begin{bmatrix} \gamma_0 \\ \gamma_1 \end{bmatrix} & = \begin{bmatrix} \frac{1}{2} (\gamma_1 + 1) \\ -\frac{1}{2} \gamma_0 \end{bmatrix} - \begin{bmatrix} \frac{1}{2} & 0 \\ 0 & \frac{1}{2} \end{bmatrix} \int_{0}^{x} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \, \begin{bmatrix} \gamma_0(s) \\ \gamma_1(s) \end{bmatrix} \, ds + + + Stopping Conditions ------------------- diff --git a/idesolver/idesolver.py b/idesolver/idesolver.py index ad4d4d2..ca1f1a1 100644 --- a/idesolver/idesolver.py +++ b/idesolver/idesolver.py @@ -141,7 +141,7 @@ def __init__( Defaults to :math:`k(x, s) = 1`. f : The function :math:`F(y)`. - Defaults to :math:`f(y) = 0`. + Defaults to :math:`F(y) = 0`. lower_bound : The lower bound function :math:`\\alpha(x)`. Defaults to the first element of ``x``. @@ -320,7 +320,7 @@ def solve(self, callback: Optional[Callable] = None) -> np.ndarray: ) ) break - except (np.ComplexWarning, TypeError) as e: + except np.ComplexWarning as e: raise exceptions.UnexpectedlyComplexValuedIDE( "Detected complex-valued IDE. Make sure to pass y_0 as a complex number." ) from e @@ -368,7 +368,13 @@ def _solve_rhs_with_known_y(self, y: np.ndarray) -> np.ndarray: def integral(x): def integrand(s): - return self.k(x, s) * self.f(interpolated_y(s)) + k = self.k(x, s) + f = self.f(interpolated_y(s)) + + if k.ndim == 1: + return k * f + else: + return k @ f result = [] for i in range(self.y_0.size): @@ -383,7 +389,12 @@ def integrand(s): return coerce_to_array(result) def rhs(x, y): - return self.c(x, interpolated_y(x)) + (self.d(x) * integral(x)) + c = self.c(x, interpolated_y(x)) + d = self.d(x) + if d.ndim == 1: + return c + (d * integral(x)) + else: + return c + (d @ integral(x)) return self._solve_ode(rhs) diff --git a/setup.cfg b/setup.cfg index df79539..7443d52 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = idesolver -version = 1.1.0 +version = 1.2.0 description = A general purpose iterative numeric integro-differential equation (IDE) solver long_description = file: README.rst long_description_content_type = text/x-rst diff --git a/tests/test_solver_against_analytic_solutions.py b/tests/test_solver_against_analytic_solutions.py index 5e3d906..4a26852 100644 --- a/tests/test_solver_against_analytic_solutions.py +++ b/tests/test_solver_against_analytic_solutions.py @@ -116,7 +116,33 @@ upper_bound=lambda x: x, ), lambda x: [np.sin(x), np.cos(x)], - ) + ), + ( # equivalent to previous case, but with d(x) and k(x, s) as vectors + IDESolver( + x=np.linspace(0, 7, 100), + y_0=[0, 1], + c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]], + d=lambda x: [-0.5, -0.5], + k=lambda x, s: [1, 1], + f=lambda y: y, + lower_bound=lambda x: 0, + upper_bound=lambda x: x, + ), + lambda x: [np.sin(x), np.cos(x)], + ), + ( # equivalent to previous case, but with d(x) and k(x, s) as matrices + IDESolver( + x=np.linspace(0, 7, 100), + y_0=[0, 1], + c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]], + d=lambda x: [[-0.5, 0], [0, -0.5]], + k=lambda x, s: [[1, 0], [0, 1]], + f=lambda y: y, + lower_bound=lambda x: 0, + upper_bound=lambda x: x, + ), + lambda x: [np.sin(x), np.cos(x)], + ), ]