我需要求解具有特定参数 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 给我的数据['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 """用于检查所有求解器共有的参数的辅助函数。""" ----> 6 y0 = np.asarray(y0)
ValueError: 用序列设置数组元素。请求的数组在 1 维之后具有不均匀的形状。检测到的形状是 (3,) + 非均匀部分。
对 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方法,但有些方法给了我同样的错误,而其他方法花了很长时间才完成,我放弃了。
出现的错误“ValueError:使用序列设置数组元素。请求的数组在维度 1 之后形状不一致。检测到的形状是 (3,) + 形状不一致的部分”表明正在尝试使用形状不同的序列来创建 NumPy 数组,特别是当
solve_ivp
在
root
内被调用来更新
y0
时。这通常发生在
solve_ivp
的连续调用中,
y0
没有以一致的形状传递时。
没有在
H(phi0)
函数内返回
phi0
的更新值,这可能会导致在
root
求解器进行迭代时将错误的值传递给
solve_ivp
。
以下是如何解决此问题的方法:
-
从
H(phi0)
函数返回phi0
: 修改的H(phi0)
函数,以返回计算出的phi0
值,以及data['H'][0] - h
的差值。 -
在
root
中更新phi0
: 在的root
调用中,使用从H(phi0)
返回的更新后的phi0
值更新phi0
。
以下是更新后的代码:
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import root
# ...(其他函数和定义)
def H(phi0):
data = inflate(phi0, alpha)
p = data['H'][0] - h
return p, phi0 # 返回 p 和更新后的 phi0
sol = root(lambda phi0: H(phi0)[0], 10.0) # 只使用 p 进行求根
phi0 = sol.x[0] # 从结果中获取更新后的 phi0
通过返回更新后的
phi0
并将其传递回
root
,确保了
solve_ivp
始终收到形状一致的
y0
数组,从而防止了
ValueError
。
此外,仔细检查的
inflate
函数中
y
的更新,尤其是在循环中或连续调用
solve_ivp
时。确保
y
的形状在每次传递给
solve_ivp
时都与初始条件
y0
的形状一致(在这种情况下为 (3,))。