对非齐次ODE使用scipy solve_bvp

Question

我正在尝试解决以下4阶BVP

y''=K-C*y

我的x变量是一个有100个节点的linspace。正如你所看到的，K是一个长度为100的向量，它使方程是非齐次的。然而，当我按下solve时，会出现以下错误：

Cell In [11], line 18, in fun(x, y)
17 def fun(x, y):
---> 18     ans = vector-np.multiply(C,y[0])
19     return np.vstack((y[1],y[2],y[3],ans))
ValueError: operands could not be broadcast together with shapes (100,) (99,)

为什么解算器会突然将y的长度更改1，以及如何修复此错误？

编辑：我必须补充一点，当K不存在时，解算器工作良好，即方程是齐次的。

from scipy.integrate import solve_bvp
import numpy as np

L = 10 
nodes = 100
A = 1000
B = 1500
C = 0.05
x = np.linspace(0,L,nodes)
vector = np.ones(nodes)
def fun(x, y):
ans = vector-np.multiply(C,y[0])
return np.vstack((y[1],y[2],y[3],ans))
def bc(ya, yb):
return np.array([ya[2], yb[2], ya[3]+A/B, yb[3]])
y_a = np.zeros((4, x.size))
res_a = solve_bvp(fun, bc, x, y_a)

res1 = res_a.sol(x)[0]
res2 = res_a.sol(x)[1]
res3 = B*res_a.sol(x)[2]
res4 = B*res_a.sol(x)[3]

Answer 1

解算器在第一轮中为第一细分的nodes-1=99段建立多项式近似系统。

不能保证细分在以后的解算器循环中保持不变。因此，ODE右侧函数必须使用任意的x数组。这意味着作为函数表给出的参数需要针对通用x阵列进行插值。numpy.interp中有用于瞬时插值的过程，scipy.interpolate.interp1d中有用于生成插值函数的过程。