我在使用scipy的solve_bvp时遇到问题

v1uwarro 于 2023-08-05 发布在其他

关注(0)|答案(1)|浏览(199)

几周来，我一直试图将ODE的这一部分翻译成python，但在运行它时，我经常遇到不同的错误。我认为我对solve_bvp的理解是一般的，但我不能弄清楚我做错了什么。我一直在阅读库的文档，了解各个函数是如何工作的，但几周后我觉得自己并没有真正理解它。
问题：

的数据
这就是我试图重新创建的数学，下面是我想出的代码。

import numpy as np
import scipy

#Constants#
P1 = 0.05
P2 = 0.1
Tau = 0.5
W = 1
Uh = 1
Uv = 1.562
Bv1 = 1.59
Bv2 = 1.59
A1 = 1.2
A2 = 0.8
Roe1 = 0.177
Roe2 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
Qv = 100
H = 1

#g is lamda this variable will change with time#
#This is variable test area#

g = 1

BH1 = 1
BH2 = 2

x1 = np.array([0, 1])
x2 = (0, 1)
t = np.linspace(0, 1, num=365)
#############################################

i = input('What is the value of i?')
i = ([int(i)])
#Functions#

def Bhinv(i):
    if i == 1:
        l = BH1*min(Nv/(Qv*int(Nh)), 1)+int(Roe1)*H*min(Nv/(Qh*int(Nh)), 1)
        return l
    if i == 2:
        l = BH2*min(Nv/(Qv*Nh), 1)+Roe2*H*min(Nv/(Qh*Nh), 1)
        return l
    else:
        return print("Null")

def sys1(x1, x2, i):
    if i.all() == 1:
        dx1dt = ((int(P1)/g)-1)*Uh*x1[0] + \
            (((int(Bhinv(1)))/g)*(int(Nh)/Nv)*x2[0])
        dx2dt = ((int(Bv1)/g)*(Nv/Nh)*x1[0])-int(Uv)*x2[0]
        return dx1dt, dx2dt

    if i.all() == 2:
        dx1dt = ((P2/g)-1)*Uh*x1(t)+(((BH2)/g)*(Nh/Nv)*x2(t))
        dx2dt = ((Bv2/g)*(Nv/Nh)*x1(t))-Uv*x2(t)
        return dx1dt, dx1dt

# BRN(Bhi) (i = 1,2)

def step1(x1, x2):
    x1 = 1
    x2 = 0
    return np.array([1, 0, 0])

y = ([0, 0],
     [0, 0])
v1s1 = scipy.integrate.solve_bvp(sys1, step1, x1, y, x2, y)

字符串
我正在处理的当前错误是IndexError: too many indices for array。

scipy

来源：https://stackoverflow.com/questions/76615402/i-am-having-trouble-with-scipys-solve-bvp

1条答案

按热度按时间

vuktfyat1#

这不是边值问题，而是初值问题。可以使用scipy.integrate.solve_ivp解决这些问题。由于这是一个很大的问题（您需要针对i=1和i=2求解它，并优化值beta_h1、beta_h2和lambda以实现某些条件），因此我对所有这些未知数使用值1，并仅针对i=1的情况实现求解器。

import numpy as np
from scipy.integrate import solve_ivp

p1 = 0.05
p2 = 0.1
tau = 0.5
omega = 1
mu_h = 1
mu_v = 1.562
beta_v1 = 1.59
beta_v2 = 1.59
a1 = 1.2
a2 = 0.8
rho1 = 0.177
rho2 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
Qv = 100
H = 1

beta_h1 = 1
beta_h2 = 2

lamb = 1

def beta_h1_func(beta_h):
    return beta_h*min(Nv/Qv/Nh, 1) + rho1*H*min(Nv/Qh/Nh, 1)

def beta_v1_func(t):
    p = a1 if t < tau else (omega - a1*tau)/(omega-tau)
    return beta_v1*min(Qv*Nh/Nv, 1)*p

def ode1(t, x):
    x1, x2 = x
    dx1dt = (rho1/lamb - 1)*mu_h*x1 + (beta_h1_func(beta_h1)*Nh/lamb/Nv)*x2
    dx2dt = (beta_v1_func(t)*Nv/lamb/Nh)*x1 - mu_v*x2
    return [dx1dt, dx2dt]

x01 = np.array([1, 0])
x02 = np.array([0, 1])
t = np.linspace(0, omega, 100)
sol1 = solve_ivp(ode1, [t.min(), t.max()], x01, t_eval=t)
v1 = sol1.y[:, -2:-1]
sol2 = solve_ivp(ode1, [t.min(), t.max()], x02, t_eval=t)
v2 = sol2.y[:, -2:-1]

rho = np.max(np.linalg.eigvals(np.hstack((v1, v2))))
print(rho)  # 0.6812599731032912

字符串
这个代码并不完美。你有大量的变量，通常会传递给函数，但我采取了懒惰的路线，把它们留在全局范围内。这不是一个好的实践，但这是我为了简单而做的。

赞(0）回复(0）举报 2023-08-05

我来回答

我在使用scipy的solve_bvp时遇到问题

1条答案

相关问题

热门标签

最新问答