Allow d(x) and k(x, s) to be matrices

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 <https://github.com/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 <https://github.com/nbrucy>`_!)
+* Add support for multidimensional IDEs (PR :pr:`35` resolves :issue:`28`, thanks `nbrucy <https://github.com/nbrucy>`_!).
 
 v1.0.5
 ------

docs/source/index.rst

@@ -16,12 +16,15 @@ The algorithm is based on a scheme devised by `Gelmi and Jorquera <https://doi.o
 
 If you use IDESolver in your work, please consider `citing it <https://doi.org/10.21105/joss.00542>`_.
 
+IDESolver is open-source software developed on `GitHub <https://github.com/JoshKarpel/idesolver>`_.
+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.

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
 -------------------
 

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)
 

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

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)],
+    ),
 ]