在root或fsolve中使用solve_ivp会给出ValueError：使用序列设置数组元素

Question

我需要求解具有特定参数 p 的方程组，然后需要找到能够给出所需结果的 p 值。我的代码看起来像（简化版本）

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import root

def system(t, y, alpha):
    phi, psi, N = y
    dphi_dt = psi
    dpsi_dt = -3 * H(phi, psi, alpha) * psi - dV_dphi(phi, alpha)
    dN_dt = H(phi, psi, alpha)
    return [dphi_dt, dpsi_dt, dN_dt]

# initial conditions
psi0 = -dV_dphi(phi0, alpha)/(3*H(phi0,0, alpha))  # Initial value of psi
N0 = 0.0    # Initial value of N

# time span
t0 = 0
dt = 0.01  # Time step
t = t0

# Initialize lists to store the results
t_values = [t0]
phi_values = [phi0]
psi_values = [psi0]
N_values = [N0]


# Initial conditions
y = [phi0, psi0, N0]

# Solve the system of differential equations for one step
sol = solve_ivp(system, [t, t+dt], y, args=(alpha,))
y = sol.y[:, -1]
t = sol.t[-1]

# Extract solutions
phi, psi, N = y

# calculate extra info
# save data into a dictionary (containing lists) e.g.
data = {'H':[1,2,3,4,5]}
return data

这段代码运行完美，给了我想要的结果。

稍后，我需要修复参数 alpha 并找到 phi0 给我的 data['H'][0] == h，例如。我正在用

def H(phi0):
    data = inflate(phi0, alpha)
    p = data['H'][0] - h
    return p

sol = root(H, 10.0)

其中 data['H'] 是我之前计算的列表之一。正如预期的那样，函数 H 返回一个数字。当调用root（或fsolve，结果是相同的）时，代码调用inflate，然后inflate调用solve_ivp。此时，solve_ivp 抱怨道：

def check_arguments(fun, y0, support_complex):
      5     """Helper function for checking arguments common to all solvers."""
----> 6     y0 = np.asarray(y0)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

对 inflate 函数的调用甚至从未到达返回点。 solve_ivp 一次也没有结束。据我了解，在solve_ivp的迭代中，序列以某种方式作为初始条件y0内的参数传递。仅当在根内部调用solve_ivp时才会发生这种情况，而如果我单独调用solve_ivp则不会发生这种情况。

我不知道为什么会发生这种情况或如何解决这个问题。任何帮助表示赞赏

我尝试更改 fsolve 的 root 和 odeint 的solve_ivp。使用序列设置数组元素时返回相同的错误。但错误的来源不同。我还尝试验证函数 H 是否正常工作并返回一个数字，确实如此。通过阅读文档，我确信当在 inflate 内部调用solve_ivp 时会出现错误。函数 inflate 永远不会到达其返回点，因此仅在 root 内部进行一次迭代。我确信solve_ivp在inflate中可以正常工作，因为它确实给出了所需的结果。仅当在 root 内部调用时才会出现问题。

我尝试更改solve_ivp方法，但有些方法给了我同样的错误，而其他方法花了很长时间才完成，我放弃了。

Answer 1

我想，我可以在这里看到问题，但我无法提供完整的代码，因为您在复制的脚本中省略了“H”和“dV_dphi”函数的定义。

因此，在您的脚本中，“phi0”和“psi0”被定义为 numpy 数组，这就是当“solve_ivp”尝试将 y0 转换为 numpy 数组时导致问题的原因。

要解决此问题，应确保“N0”、“phi0”和“psi0”是标量，而不是 numpy 数组。为此，可以使用“浮动”。

定义“系统”后，我将继续如下代码，并合并这些更改：

    ###your code including def system(t, y, alpha)...

    def inflate(phi0, alpha):
    phi0 = float(phi0)  # Ensure phi0 is a scalar
    psi0 = float(-dV_dphi(phi0, alpha) / (3 * H(phi0, 0, alpha)))  # Ensure psi0 is a scalar
    N0 = 0.0  # Initial value of N
    
    # initial conditions
    y = [phi0, psi0, N0]

    # time span
    t0 = 0
    dt = 0.01  # Time step
    t = t0


    print(f"Initial conditions: y = {y}")


    sol = solve_ivp(system, [t, t + dt], y, args=(alpha,))
    y = sol.y[:, -1]
    t = sol.t[-1]


    phi, psi, N = y


    data = {'H': [H(phi, psi, alpha)]} 


    print(f"Solved values: phi = {phi}, psi = {psi}, N = {N}")
    print(f"Data: {data}")

    return data

def H_function(phi0, alpha, h):
    data = inflate(phi0, alpha)
    p = data['H'][0] - h
    return p

# Example values 
alpha = 1.0
h = 1.0

# finding the optimal phi0
sol = root(H_function, 10.0, args=(alpha, h))

print(f"Root finding result: {sol}")


# extracting the optimal phi0
optimal_phi0 = sol.x[0]
print(f"Optimal phi0: {optimal_phi0}")

我对上述函数使用了任意定义，并且代码有效。