Numpy操作数不能与形状（3，）（2，）（3，）一起广播

Question

我正在学习微分方程课程，并且代码是用于SIR模型的Runge Kutta方法。我不是一个出色的程序员，并且在Explicit_Runge_Kutta类的调用部分遇到了dY_j += dt*a[j,l]*ks[l]的问题。它返回：

ValueError                                Traceback (most recent call last)
<ipython-input-30-467a168a5b4e> in <module>
     69 T = 24*30 #Nmax*dt
     70 
---> 71 ts, ys = rk4(y0, t0, T, sir, Nmax)
     72 
     73 plt.plot(ts, ys)

<ipython-input-30-467a168a5b4e> in __call__(self, y0, t0, T, f, Nmax)
     26                 dY_j = np.zeros_like(y, dtype=np.double)
     27                 for l in range(j):
---> 28                     dY_j += dt*a[j,l]*ks[l]
     29 
     30                 ks[j] = f(t_j, y + dY_j)

ValueError: operands could not be broadcast together with shapes (3,) (2,) (3,)

这是我的完整代码：

class Explicit_Runge_Kutta:
    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c

    def __call__(self, y0, t0, T, f, Nmax):
        # Extract Butcher table
        a, b, c = self.a, self.b, self.c

        # Stages
        s = len(b)
        ks = [np.zeros_like(y0, dtype=np.double) for s in range(s)]

        # Start time-stepping
        ys = [y0]
        ts = [t0]
        dt = (T - t0)/Nmax

        while(ts[-1] < T):
            t, y = ts[-1], ys[-1]

            # Compute stages derivatives k_j
            for j in range(s):
                t_j = t + c[j]*dt
                dY_j = np.zeros_like(y, dtype=np.double)
                for l in range(j):
                    dY_j += dt*a[j,l]*ks[l]

                ks[j] = f(t_j, y + dY_j)

            # Compute next time-step
            dy = np.zeros_like(y, dtype=np.double)
            for j in range(s):
                dy += dt*b[j]*ks[j]

            ys.append(y + dy)
            ts.append(t + dt)

        return (np.array(ts), np.array(ys))

class SIR:
    def __init__(self, beta, gamma):
        self.beta = beta
        self.gamma = gamma

    def __call__(self, t, y):
        return np.array([-beta*y[0]*y[1], beta*y[0]*y[1]-gamma*y[0]*y[1]])

beta = 10/(40*8*24)
gamma = 3/(15*24)
sir = SIR(beta, gamma)

a = np.array([[0, 0, 0, 0],
              [0.5, 0, 0, 0],
              [0, 0.5, 0, 0],
              [0, 0, 1, 0]])
b = np.array([1/6, 1/3, 1/3, 1/6])
c = np.array([0, 0.5, 0.5, 1])

rk4 = Explicit_Runge_Kutta(a,b,c)

t0 = 0
y0 = [50, 1, 0]

dt = 0.1 #6 min
D = 30 #days
Nmax = 10*30*24 #steps
T = 24*30 #Nmax*dt

ts, ys = rk4(y0, t0, T, sir, Nmax)

plt.plot(ts, ys)
plt.plot(ts, ys[:,0]+ys[:,1]+ys[:,2])
plt.legend(["S","I","R", "S+I+R"])

抱歉，如果这还不足以解决此问题，我本人将无法理解所有这些信息。 Explicit_Runge_Kutta是在我们的任务中分配给我们的，其余部分受我们的jupyterlab讲义中的内容严重影响。

Answer 1

您的SIR.__call__仅返回大小为2的向量，因此ks[l]的大小仅为2。我想应该是

def __call__(self, t, y):
        return np.array([-beta*y[0]*y[1], beta*y[0]*y[1]-gamma*y[0]*y[1], gamma*y[0]*y[1]])

代替