【发布时间】：2020-05-07 15:37:29
【问题描述】：
当k!=0solve_ivp 给我这个错误ValueError: setting an array element with a sequence.
我认为问题是当k!=0 fkm 变成一个数组并且我试图放入sol 数组时
但我找不到解决问题的方法。有什么想法吗？
编辑：我尝试了vectorized=True，但没有成功。
import numpy as np
from scipy.integrate import solve_bvp, solve_ivp
from numpy import sin, cos, pi
from scipy.special import binom
E = 200e9
nu = 0.3
Q = -10e3
a = 1
b = 2
t0 = 0.01
D0 = E*t0**3/(12*(1-nu**2))
e = 0
def t(y): return t0*(1+e*y)
def D(y): return D0*(1+e*y)**3
def D(i,y): return D0*binom(3,i)*y**i
def Dy(i,y): return D0*binom(3,i)*(i)*y**(i-1)
def Dyy(i,y): return D0*binom(3,i)*(i)*(i-1)*y**i-2
N = 3
N_nodes = 10
x = np.linspace(0, a, N_nodes)
y = np.linspace(0, b, N_nodes)
Y_arr = np.zeros((N_nodes, 5, N_nodes))
def Y(k,m):
if (k==0):
fkm = 4*Q/(m*pi*D0)
else:
fkm = np.zeros((N_nodes))
for i in range(1,k+1):
Y, Yy, Yyy, Yyyy, Yyyyy = Y_arr[k-i]
A = D(i,y)*Yyyyy + 2*Dy(i,y)*Yyyy-2*(m*pi/a)**2*Dy(i,y)*Yy
B = (Dyy(i,y) - 2*D(i,y)*(m*pi/a)**2)*Yyy
C = (D(i,y)*(m*pi/a)**4-nu*(m*pi/a)**2*Dyy(i,y))*Y
fkm += -1/D0*(A+B+C)
def dU_dy(y,U):
sol = np.zeros((4,N_nodes))
sol = [U[1],U[2],U[3],fkm + 2*(m*pi/a)**2*U[2]-(m*pi/a)**4*U[0]]
return sol
def BCs(y0, yb):
Y0, Y0y, Y0yy, Y0yyy = y0
Yb, Yby, Ybyy, Ybyyy = yb
return [Y0, Y0yy, Yb, Ybyy]
Y_guess = solve_ivp(dU_dy,(0,b),[0,1,0,1],t_eval=y,vectorized=False).y
Y, Yy, Yyy, Yyyy = solve_bvp(dU_dy, BCs, y, Y_guess, max_nodes=N_nodes).y
Yyyyy = fkm + 2*(m*pi/a)**2*Yyy-(m*pi/a)**4*Y
Y_arr[k] = np.array([Y, Yy, Yyy,Yyyy, Yyyyy])
return Y
w = np.zeros((N_nodes, N_nodes))
for k in range(N):
wk = np.zeros((N_nodes,N_nodes))
for m in range(1,N+1):
wk += np.outer(Y(k,m),sin(m*pi*x/a))
w += wk*(e**k)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-58-01b7605c068b> in <module>()
38 wk = np.zeros((N_nodes,N_nodes))
39 for m in range(1,N+1):
---> 40 wk += np.outer(Y(k,m),sin(m*pi*x/a))
41 w += wk*(e**k)
<ipython-input-58-01b7605c068b> in Y(k, m)
26 return [Y0, Y0yy, Yb, Ybyy]
27
---> 28 Y_guess = solve_ivp(dU_dy,(0,b),[0,1,0,1],t_eval=y,vectorized=True).y
29 Y, Yy, Yyy, Yyyy = solve_bvp(dU_dy, BCs, y, Y_guess, max_nodes=N_nodes).y
30 Yyyyy = fkm + 2*(m*pi/a)**2*Yyy-(m*pi/a)**4*Y
d:\python\python36\lib\site-packages\scipy\integrate\_ivp\ivp.py in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
541 method = METHODS[method]
542
--> 543 solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
544
545 if t_eval is None:
d:\python\python36\lib\site-packages\scipy\integrate\_ivp\rk.py in __init__(self, fun, t0, y0, t_bound, max_step, rtol, atol, vectorized, first_step, **extraneous)
93 self.max_step = validate_max_step(max_step)
94 self.rtol, self.atol = validate_tol(rtol, atol, self.n)
---> 95 self.f = self.fun(self.t, self.y)
96 if first_step is None:
97 self.h_abs = select_initial_step(
d:\python\python36\lib\site-packages\scipy\integrate\_ivp\base.py in fun(t, y)
137 def fun(t, y):
138 self.nfev += 1
--> 139 return self.fun_single(t, y)
140
141 self.fun = fun
d:\python\python36\lib\site-packages\scipy\integrate\_ivp\base.py in fun_single(t, y)
124 if vectorized:
125 def fun_single(t, y):
--> 126 return self._fun(t, y[:, None]).ravel()
127 fun_vectorized = self._fun
128 else:
d:\python\python36\lib\site-packages\scipy\integrate\_ivp\base.py in fun_wrapped(t, y)
19
20 def fun_wrapped(t, y):
---> 21 return np.asarray(fun(t, y), dtype=dtype)
22
23 return fun_wrapped, y0
d:\python\python36\lib\site-packages\numpy\core\_asarray.py in asarray(a, dtype, order)
83
84 """
---> 85 return array(a, dtype, copy=False, order=order)
86
87
ValueError: setting an array element with a sequence.
【问题讨论】：
-
问题出在
dU_dy的返回值上。根据文档，它应该是一个形状为 (n,) 或 (n,k) 的数值数组。您返回一个列表，强制solve_ivp为您进行转换，但它失败了（不一定在第一次迭代时）。 -
@hpaulj 如何返回一个形状为 (n,k) 的数组而不失败？
-
找出失败的原因，然后解决这个问题！
-
@hpaulj 我试图解决两天这就是我在这里的原因。我在文档中找不到如何返回 (n,k) 数组。我试过
vectorized=True