Как обработать «ValueError: установка элемента массива с последовательностью». ошибка

Question

Когда k!=0 resolve_ivp дает мне эту ошибку ValueError: setting an array element with a sequence.

Я думаю, проблема в том, что k!=0 fkm становится массивом, а я пытаюсь вставить sol array

Но я не мог найти способ решить проблему. Есть идеи?

Изменить: я пробовал 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.

Answer 1

Когда я помещаю:

U = np.array([[0],[1],[0],[1]])
print(k, dU_dy(0, U))
print(np.array(dU_dy(0, U), float))

перед вызовом solve_ivp, я получаю:

1323:~/mypy$ python3 stack61661864.py 
0 [array([1]), array([0]), array([1]), array([-0.69518879])]
[[ 1.        ]
 [ 0.        ]
 [ 1.        ]
 [-0.69518879]]
0 [array([1]), array([0]), array([1]), array([-0.3475944])]
[[ 1.       ]
 [ 0.       ]
 [ 1.       ]
 [-0.3475944]]
0 [array([1]), array([0]), array([1]), array([-0.2317296])]
[[ 1.       ]
 [ 0.       ]
 [ 1.       ]
 [-0.2317296]]
1 [array([1]), array([0]), array([1]), array([-0.0819576 ,  0.0122602 ,  0.00429568,  0.00517298,  0.00733516,
        0.00895909,  0.00857527,  0.00368501,  0.060713  ,  1.47233495])]
TypeError: only size-1 arrays can be converted to Python scalars

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "stack61661864.py", line 72, in <module>
    wk += np.outer(Y(k,m),sin(m*pi*x/a))
  File "stack61661864.py", line 57, in Y
    print(np.array(dU_dy(0, U), float))
ValueError: setting an array element with a sequence.

Он отлично работает для k=0, sol - это список из 4 (1,) массивов фигур, что составляет массив (4,1).

Но для k=1 4-й член представляет собой массив фигур (10,). Этот sol нельзя превратить в массив.

Для k=0, fkm = 4*Q/(m*pi*D0). Это должен быть скаляр. Но для других k это fkm = np.zeros((N_nodes))

Как я уже сказал, вам нужно убедиться, что функция ivp/bvp возвращает массив правильной формы. Я тестировал sol с множеством операторов print. Вызов dU_dy сам по себе, как я здесь демонстрирую, является хорошим способом его проверки (при условии, что входные данные совпадают с тем, что отправляет решатель).

Проверьте формы и проверьте их снова!

Answer 2

Могу я предложить, чтобы оба ваших случая k==0 и k!=0 возвращали один и тот же тип данных и форму, это поможет решателю. Я не знаю специфики вашей проблемы, но если fkm должен быть массивом в одном случае, тогда это должен быть массив в другом случае, и вы можете уточнить, что это означает.

Я пытался отредактировать вашу функцию Y, но не понимая специфики проблемы, я не уверен, что смогу помочь больше:

def Y(k,m):
    if (k==0):
        # for k==0 make fkm a list of N_nodes with same value
        fkm = [4*Q/(m*pi*D0)]*N_nodes
    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))
        # then add all N_node values of fkm to solver (so all 10 are there in this case)
        sol = [U[1],U[2],U[3]] + (fkm + 2*(m*pi/a)**2*U[2]-(m*pi/a)**4*U[0]).tolist()
        return sol
    def BCs(y0, yb):
        Y0, Y0y, Y0yy, Y0yyy = y0
        Yb, Yby, Ybyy, Ybyyy = yb
        return [Y0, Y0yy, Yb, Ybyy]

# then N_node guesses are need for solver
Y_guess = solve_ivp(dU_dy,(0,b),[0,1,0] +[1]*N_nodes,t_eval=y,vectorized=False).y
Y, Yy, Yyy, Yyyy = solve_bvp(dU_dy, BCs, y, Y_guess, max_nodes=N_nodes).y

Для k==0 и k==1 solve_ivp линия проходит без заминок, однако теперь нужно посмотреть solve_bvp.