Python solve_bvp报错：operands could not be broadcast together (100,) (99,)

Let's break down why you're hitting this ValueError and fix it quickly:

The Root Cause #

You're hardcoding a fixed-length array (array=np.empty(100,)) in your fun function, but solve_bvp doesn't always pass 100 points to fun during its internal calculations.

When you start with an initial mesh of 100 points, the solver will internally compute values at midpoints of each mesh interval (that's 99 points total) as part of its collocation method. At that point, y[0] will have shape (99,), but your var array is stuck at (100,)—hence the broadcast error.

The Fix #

Instead of hardcoding the array length, generate it based on the number of points in the x argument passed to fun (which matches the number of columns in y). Here's the adjusted code with the key change highlighted:

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

# Added missing L definition (required for xmesh to work)
L = 1.0

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)
        # ✅ Changed from fixed 100 to len(x) to match current input points
        array=np.empty(len(x),) 
        array.fill(1)
        def array_function():
            return array
        var= array_function() # var is now same length as y[0]
        rhs= var+y[0] # Now shapes match!
        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)

Extra Note #

I also added a definition for L since it was missing in your original code (it's used in xmesh but wasn't declared—this would have caused another error once the shape issue was fixed).

This adjustment ensures var always matches the length of y[0], no matter how many points solve_bvp passes to fun during its internal iterations.

内容的提问来源于stack exchange，提问作者Anweshan