From af583dcab154366f518d38323dcd4dccdf055c6c Mon Sep 17 00:00:00 2001 From: Josh Karpel Date: Sun, 3 Oct 2021 14:32:02 -0500 Subject: [PATCH 1/4] let k be a matrix --- docs/source/changelog.rst | 7 +++++++ idesolver/idesolver.py | 12 +++++++++--- setup.cfg | 2 +- tests/test_solver_against_analytic_solutions.py | 15 ++++++++++++++- 4 files changed, 31 insertions(+), 5 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 9012915..44a1be8 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -3,6 +3,13 @@ Change Log .. currentmodule:: idesolver +v1.2.0 +------ + +* :math:`k(x, s)` can now be a matrix. +* `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 `_!) diff --git a/idesolver/idesolver.py b/idesolver/idesolver.py index ad4d4d2..a4869fa 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.size == 1: + return k * f + else: + return k @ f result = [] for i in range(self.y_0.size): 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..1ba3ee4 100644 --- a/tests/test_solver_against_analytic_solutions.py +++ b/tests/test_solver_against_analytic_solutions.py @@ -116,7 +116,20 @@ upper_bound=lambda x: x, ), lambda x: [np.sin(x), np.cos(x)], - ) + ), + ( # equivalent to previous case, but with k(x, s) as a matrix instead of a scalar + 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, + k=lambda x, s: np.array([[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)], + ), ] From 1169b3194585de9a805ffcaedc5fae40bb0bc453 Mon Sep 17 00:00:00 2001 From: Josh Karpel Date: Sun, 3 Oct 2021 14:43:50 -0500 Subject: [PATCH 2/4] let d(x) be a matrix too --- docs/source/changelog.rst | 2 +- idesolver/idesolver.py | 7 ++++++- tests/test_solver_against_analytic_solutions.py | 6 +++--- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 44a1be8..6b804dc 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -6,7 +6,7 @@ Change Log v1.2.0 ------ -* :math:`k(x, s)` can now be a matrix. +* :math:`d(x)` and :math:`k(x, s)` can now be matrices. * `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. diff --git a/idesolver/idesolver.py b/idesolver/idesolver.py index a4869fa..975e621 100644 --- a/idesolver/idesolver.py +++ b/idesolver/idesolver.py @@ -389,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.size == 1: + return c + (d * integral(x)) + else: + return c + (d @ integral(x)) return self._solve_ode(rhs) diff --git a/tests/test_solver_against_analytic_solutions.py b/tests/test_solver_against_analytic_solutions.py index 1ba3ee4..4798dbf 100644 --- a/tests/test_solver_against_analytic_solutions.py +++ b/tests/test_solver_against_analytic_solutions.py @@ -117,13 +117,13 @@ ), lambda x: [np.sin(x), np.cos(x)], ), - ( # equivalent to previous case, but with k(x, s) as a matrix instead of a scalar + ( # 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, - k=lambda x, s: np.array([[1, 0], [0, 1]]), + 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, From 0001bd6f7652f4c94b51a850c4e4eb584644449f Mon Sep 17 00:00:00 2001 From: Josh Karpel Date: Mon, 4 Oct 2021 18:13:56 -0500 Subject: [PATCH 3/4] allow vector d and k as well --- idesolver/idesolver.py | 4 ++-- tests/test_solver_against_analytic_solutions.py | 13 +++++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/idesolver/idesolver.py b/idesolver/idesolver.py index 975e621..ca1f1a1 100644 --- a/idesolver/idesolver.py +++ b/idesolver/idesolver.py @@ -371,7 +371,7 @@ def integrand(s): k = self.k(x, s) f = self.f(interpolated_y(s)) - if k.size == 1: + if k.ndim == 1: return k * f else: return k @ f @@ -391,7 +391,7 @@ def integrand(s): def rhs(x, y): c = self.c(x, interpolated_y(x)) d = self.d(x) - if d.size == 1: + if d.ndim == 1: return c + (d * integral(x)) else: return c + (d @ integral(x)) diff --git a/tests/test_solver_against_analytic_solutions.py b/tests/test_solver_against_analytic_solutions.py index 4798dbf..4a26852 100644 --- a/tests/test_solver_against_analytic_solutions.py +++ b/tests/test_solver_against_analytic_solutions.py @@ -117,6 +117,19 @@ ), 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), From 14c8fe7ad330811df692fc27a852a2d27ded51da Mon Sep 17 00:00:00 2001 From: Josh Karpel Date: Mon, 4 Oct 2021 18:39:48 -0500 Subject: [PATCH 4/4] give the docs some love --- docs/source/changelog.rst | 4 +- docs/source/index.rst | 7 +++- docs/source/manual.rst | 79 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 4 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 6b804dc..b583b36 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -6,13 +6,13 @@ Change Log v1.2.0 ------ -* :math:`d(x)` and :math:`k(x, s)` can now be matrices. +* :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 -------------------