Repository: scipy/scipy
Issue category: Bug Report

BUG: integrate.trapezoid: broadcasting failure after #21524
### Describe your issue.

`scipy.integrate.trapezoid` fails when 2D arrays of identical shape are passed as `x` and `y` arguments.

This failure started in recent nightlies, probably from https://github.com/scipy/scipy/pull/21524.

The same example works with scipy 1.14.1 and earlier.

@Armavica @lucascolley 

### Reproducing Code Example

```python
In [1]: import numpy as np, scipy

In [2]: x = np.random.randn(3, 10)

In [3]: scipy.integrate.trapezoid(x, x)
```


### Error message

```shell
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 scipy.integrate.trapezoid(x, x)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/integrate/_quadrature.py:142, in trapezoid(y, x, dx, axis)
    140     shape = [1] * y.ndim
    141     shape[axis] = y.shape[axis]
--> 142     x = xp.reshape(x, shape)
    143     d = x[tuple(slice1)] - x[tuple(slice2)]
    144 try:

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/_lib/array_api_compat/_internal.py:28, in get_xp.<locals>.inner.<locals>.wrapped_f(*args, **kwargs)
     26 @wraps(f)
     27 def wrapped_f(*args, **kwargs):
---> 28     return f(*args, xp=xp, **kwargs)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/scipy/_lib/array_api_compat/common/_aliases.py:382, in reshape(x, shape, xp, copy, **kwargs)
    380     y.shape = shape
    381     return y
--> 382 return xp.reshape(x, shape, **kwargs)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/numpy/core/fromnumeric.py:285, in reshape(a, newshape, order)
    200 @array_function_dispatch(_reshape_dispatcher)
    201 def reshape(a, newshape, order='C'):
    202     """
    203     Gives a new shape to an array without changing its data.
    204
   (...)
    283            [5, 6]])
    284     """
--> 285     return _wrapfunc(a, 'reshape', newshape, order=order)

File ~/.pyenv/versions/3.12.3/envs/py3.12.3/lib/python3.12/site-packages/numpy/core/fromnumeric.py:59, in _wrapfunc(obj, method, *args, **kwds)
     56     return _wrapit(obj, method, *args, **kwds)
     58 try:
---> 59     return bound(*args, **kwds)
     60 except TypeError:
     61     # A TypeError occurs if the object does have such a method in its
     62     # class, but its signature is not identical to that of NumPy's. This
   (...)
     66     # Call _wrapit from within the except clause to ensure a potential
     67     # exception has a traceback chain.
     68     return _wrapit(obj, method, *args, **kwds)

ValueError: cannot reshape array of size 30 into shape (1,10)
```


### SciPy/NumPy/Python version and system information

```shell
1.15.0.dev0+git20241118.0ea90cf 1.26.4 sys.version_info(major=3, minor=12, micro=3, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  lapack:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.13.6
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 15.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 15.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 13.3.0
  pythran:
    include directory: ../../../../../../private/var/folders/g6/rgtlsw6n123b0gt5483s5_cm0000gn/T/pip-build-env-zkcj2f79/overlay/lib/python3.12/site-packages/pythran
    version: 0.17.0
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /private/var/folders/g6/rgtlsw6n123b0gt5483s5_cm0000gn/T/cibw-run-7dkhp5o6/cp312-macosx_arm64/build/venv/bin/python
  version: '3.12'
```



Functions to edit:
scipy/integrate/_quadrature.py:trapezoid

Patch:
diff --git a/scipy/integrate/_quadrature.py b/scipy/integrate/_quadrature.py
index 4f423da63895..65faaaa79b85 100644
--- a/scipy/integrate/_quadrature.py
+++ b/scipy/integrate/_quadrature.py
@@ -136,11 +136,16 @@ def trapezoid(y, x=None, dx=1.0, axis=-1):
         d = dx
     else:
         x = _asarray(x, xp=xp, subok=True)
-        # reshape to correct shape
-        shape = [1] * y.ndim
-        shape[axis] = y.shape[axis]
-        x = xp.reshape(x, shape)
-        d = x[tuple(slice1)] - x[tuple(slice2)]
+        if x.ndim == 1:
+            d = x[1:] - x[:-1]
+            # make d broadcastable to y
+            slice3 = [None] * nd
+            slice3[axis] = slice(None)
+            d = d[tuple(slice3)]
+        else:
+            # if x is n-D it should be broadcastable to y
+            x = xp.broadcast_to(x, y.shape)
+            d = x[tuple(slice1)] - x[tuple(slice2)]
     try:
         ret = xp.sum(
             d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0,

Test patch: diff --git a/scipy/integrate/tests/test_quadrature.py b/scipy/integrate/tests/test_quadrature.py
index e49ded996a2e..ec554d3b0bc8 100644
--- a/scipy/integrate/tests/test_quadrature.py
+++ b/scipy/integrate/tests/test_quadrature.py
@@ -309,14 +309,6 @@ def test_ndim(self, xp):
         r = trapezoid(q, x=z[None, None,:], axis=2)
         xp_assert_close(r, qz)
 
-        # n-d `x` but not the same as `y`
-        r = trapezoid(q, x=xp.reshape(x[:, None, None], (3, 1)), axis=0)
-        xp_assert_close(r, qx)
-        r = trapezoid(q, x=xp.reshape(y[None,:, None], (8, 1)), axis=1)
-        xp_assert_close(r, qy)
-        r = trapezoid(q, x=xp.reshape(z[None, None,:], (13, 1)), axis=2)
-        xp_assert_close(r, qz)
-
         # 1-d `x`
         r = trapezoid(q, x=x, axis=0)
         xp_assert_close(r, qx)
@@ -325,6 +317,33 @@ def test_ndim(self, xp):
         r = trapezoid(q, x=z, axis=2)
         xp_assert_close(r, qz)
 
+    @skip_xp_backends('jax.numpy',
+                      reasons=["JAX arrays do not support item assignment"])
+    @pytest.mark.usefixtures("skip_xp_backends")
+    def test_gh21908(self, xp):
+        # extended testing for n-dim arrays
+        x = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
+        y = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
+
+        out0 = xp.linspace(200, 380, 10)
+        xp_assert_close(trapezoid(y, x=x, axis=0), out0)
+        xp_assert_close(trapezoid(y, x=xp.asarray([0, 10., 20.]), axis=0), out0)
+        # x needs to be broadcastable against y
+        xp_assert_close(
+            trapezoid(y, x=xp.asarray([0, 10., 20.])[:, None], axis=0),
+            out0
+        )
+        with pytest.raises(Exception):
+            # x is not broadcastable against y
+            trapezoid(y, x=xp.asarray([0, 10., 20.])[None, :], axis=0)
+
+        out1 = xp.asarray([ 40.5, 130.5, 220.5])
+        xp_assert_close(trapezoid(y, x=x, axis=1), out1)
+        xp_assert_close(
+            trapezoid(y, x=xp.linspace(0, 9, 10), axis=1),
+            out1
+        )
+
     @skip_xp_invalid_arg
     def test_masked(self, xp):
         # Testing that masked arrays behave as if the function is 0 where