ValueError: operands could not be broadcast together with shapes (100,) (99,) error in Python

Question

I am using the bvp solver in Python to solve a 4th order boundary value problem. The actual equations are not being shown to avoid any further complexity. The code that I have written for the same has been attached below.

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

def xmesh(k1):      # xmesh definition
 return np.linspace(0,L,k1)

def solveit(constant):
   
   def fun(x,y):             # this function returns all the derivatives of y(x)
    array=np.empty(100,)     # array is an 1-D array
    array.fill(1)
    def array_function():
        return array     
    var= array_function()    # var is an 1-D array
    rhs= var+y[0]            # rhs is an 1-D array
    return [y[1],y[2],y[3],rhs]

   def bc(ya,yb):            # boundary conditions 
    return np.array([ya[0],ya[1],yb[0],yb[1]])
  
   init= np.zeros((4,len(xmesh(100))))       # initial value for the bvp solver
   sol= solve_bvp(fun,bc,xmesh(100),init,tol=1e-6,max_nodes=5000)
   arr= sol.sol(xmesh(100))[0] 
   return arr
arr= solveit(0.1)

I bump into the following error: operands could not be broadcast together with shapes (100,) (99,), every time I try to run the above code. The stack trace of the above error has also been attached below.

ValueError                                Traceback (most recent call last)
<ipython-input-52-62db777e3281> in <module>()
     24    arr= sol.sol(xmesh(100))[0]
     25    return arr
---> 26 arr= solveit(0.1)

6 frames
<ipython-input-52-62db777e3281> in solveit(constant)
     21 
     22    init= np.zeros((4,len(xmesh(100))))       # initial value for the bvp solver
---> 23    sol= solve_bvp(fun,bc,xmesh(100),init,tol=1e-6,max_nodes=5000)
     24    arr= sol.sol(xmesh(100))[0]
     25    return arr

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in solve_bvp(fun, bc, x, y, p, S, fun_jac, bc_jac, tol, max_nodes, verbose, bc_tol)
   1084                                        fun_jac_wrapped, bc_jac_wrapped, x, h)
   1085         y, p, singular = solve_newton(n, m, h, col_fun, bc_wrapped, jac_sys,
-> 1086                                       y, p, B, tol, bc_tol)
   1087         iteration += 1
   1088 

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in solve_newton(n, m, h, col_fun, bc, jac, y, p, B, bvp_tol, bc_tol)
    439     n_trial = 4
    440 
--> 441     col_res, y_middle, f, f_middle = col_fun(y, p)
    442     bc_res = bc(y[:, 0], y[:, -1], p)
    443     res = np.hstack((col_res.ravel(order='F'), bc_res))

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in col_fun(y, p)
    324 
    325     def col_fun(y, p):
--> 326         return collocation_fun(fun, y, p, x, h)
    327 
    328     def sys_jac(y, p, y_middle, f, f_middle, bc0):

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in collocation_fun(fun, y, p, x, h)
    311     y_middle = (0.5 * (y[:, 1:] + y[:, :-1]) -
    312                 0.125 * h * (f[:, 1:] - f[:, :-1]))
--> 313     f_middle = fun(x[:-1] + 0.5 * h, y_middle, p)
    314     col_res = y[:, 1:] - y[:, :-1] - h / 6 * (f[:, :-1] + f[:, 1:] +
    315                                               4 * f_middle)

/usr/local/lib/python3.6/dist-packages/scipy/integrate/_bvp.py in fun_p(x, y, _)
    648     if k == 0:
    649         def fun_p(x, y, _):
--> 650             return np.asarray(fun(x, y), dtype)
    651 
    652         def bc_wrapped(ya, yb, _):

<ipython-input-52-62db777e3281> in fun(x, y)
     14         return array
     15     var= array_function()    # var is an 1-D array
---> 16     rhs= var+y[0]            # rhs is an 1-D array
     17     return [y[1],y[2],y[3],rhs]
     18 

ValueError: operands could not be broadcast together with shapes (100,) (99,) 

As the error suggests, it is being raised because I am trying to perform some mathematical operations on two arrays of different sizes. But I don't understand why this error is even being raised over here, considering that all the arrays defined above have the same shape of (100,). Any fix to the above problem would be highly appreciated. I am stuck with this error for quite some time now.

P.S.- The functions that I have defined in the code above have no physical meaning and are completely random. The above code is just a simplistic version of a fairly complicated code that I have written. So, I do not need a correct numerical solution to the above code. All I need is the code to work fine without any error. Also, I have previously used the bvp solver in Python with success.

Answer 1

When I use print(x.shape, y[0].shape) then at some moment both change size to 99 and I thin you should use x.shape to create your array

 array = np.empty(x.shape,) # array is an 1-D array

And this works for me. But I don't know if it gives expected results.

BTW: Later I saw x.shape, y[0].shape can change size even to 3564

---

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

L = 100

def xmesh(k1):
    """xmesh definition"""
    return np.linspace(0, L, k1)

def solveit(constant):
   
    def fun(x, y):             
        """returns all the derivatives of y(x)"""
        array = np.empty(x.shape,) # array is an 1-D array
        array.fill(1)
        rhs = array + y[0]         # rhs is an 1-D array
        return [y[1], y[2], y[3], rhs]
    
    def bc(ya, yb):
        """boundary conditions"""
        return np.array([ya[0], ya[1], yb[0], yb[1]])
      
    init = np.zeros((4, len(xmesh(100))))       # initial value for the bvp solver
    sol = solve_bvp(fun, bc, xmesh(100), init, tol=1e-6, max_nodes=5000)
    arr = sol.sol(xmesh(100))[0]
    
    return arr
    
arr = solveit(0.1)