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

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 <https://github.com/nbrucy>`_!)

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

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