使用scipy.integrate.solve_bvp解决ODE时，我的函数的返回不能传播到solve_bvp的适当数组形状中

Question

我正在尝试用solve_bvp解决一个二阶常微分方程。我已经将二阶常微分方程分成两个一阶常微分方程系统。我有一组根据x(网格)值变化的常量。所以我将它们作为一个形状数组(N，)传递给我的函数numdens。在尝试运行solve_bvp时，我得到一个错误，即返回具有不同的形状，即(N，)和(N-1，)，因此不能广播到一个数组中。但是当我在函数外部手动检查每个返回时，它的形状是(N，)。如果我在不改变常量的情况下运行求解器，我会得到一个类似于正确的解。

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

E_0 = 1 * 0.0000016021773 #erg: gcm^2/s^2
m_H = 1.6*10**(-24) #g
c = 3e11 #cm
sigma_c = 2*10**(-23)
n_0 = 1*10**(20)   #1/cm^3
v_0 = (2*E_0/m_H)**(0.5)  #cm/s
T = 10**7
b = 20.3
n_eq = b*T**3
n_s = 2.03*10**(19)
Q = 1 

def velocity(v,x):
    dvdx = -sigma_c*n_0*v_0*((8*v_0*v-7*v**2-v_0**2)/(2*v*c))
    return dvdx


n_num = 100
x_num = np.linspace(-1*10**(6),3*10**(6), n_num)

sol_velo = odeint(velocity,0.999999999999*v_0,x_num)

sol_new = np.reshape(sol_velo,n_num)

def constants(v):
    D1 = (c*v/(3*n_0*v_0*sigma_c))
    D2 = ((v**2-8*v_0*v+v_0**2)/(6*v))
    D3 = sigma_c*n_0*v_0*((8*v_0*v-7*v**2-v_0**2)/(2*v*c))
    return D1,D2,D3

def numdens(x,y):
    v = sol_new
    D1,D2,D3 = constants(v)
    return np.vstack((y[1],(-D2*y[1]-D3*y[0]+Q*((1-y[0])/n_eq))/(D1)))

def bc_num(ya, yb):
    return np.array([ya[0]-n_s,yb[0]-n_eq])


y_num = np.array([np.linspace(n_s, n_eq, n_num),np.linspace(n_s, n_eq, n_num)])


sol_num = solve_bvp(numdens, bc_num, x_num, y_num)




plt.plot(sol_num.x, sol_num.y[0], label='$n(x)$')
plt.plot(x_num, sol_velo-v_0/7, label='$v(x)$')
plt.yscale('log')
plt.grid(alpha=0.5)
plt.legend(framealpha=1)
plt.show()

Answer 1

需要考虑到BVP解算器使用自适应网格。也就是说，在优化初始栅格上的初始猜测后，求解器会识别具有过大误差的区域，并在其中创建新的网格节点。据我所知，相反的情况没有实现，即使在某些应用程序中，减少特别“好”的网段上的网格节点的数量可能是明智的。

因此，您正在做的numdens函数是不可理解的，它的功能必须与传递给ODE求解器的任何其他函数完全相同。如果我必须提出一些快速解决方案，并且不知道您想要解决的底层问题是什么，我会将v的赋值更改为

v = np.interp(x,x_num,sol_velo)

因为这至少应该产生一个正确格式的数组。