Scipysolve_bvp() ändert die Arraygröße während der Iteration ⇐ Python

[Guest]

Ich versuche, eine ODE 4. Ordnung eines Balkens auf einem elastischen Fundament zu lösen. Ich wende dieses Problem auf einen seitlich belasteten Pfahl an, der in mehreren Bodenschichten eingebettet ist, die im Code angegeben werden können. Ein Array c enthält die Steifigkeitsparameter des Systems, die ich elementweise mit y[0] multiplizieren muss.

Es tritt jedoch irgendwann im Quellcode ein Problem auf y[0] wird um 1 Element kürzer und ich bekomme das Fehler:

Code: Select all

ValueError: operands could not be broadcast together with shapes (40,) (39,)

Weiß jemand, wie man das beheben kann? Ich habe versucht,
hinzuzufügen

Code: Select all

if c.size < y[0].size:
np.delete(c,-1)

direkt vor der return-Anweisung in fun(x,y) aber ohne Erfolg.

Code: Select all

# ODE-Solver ---------------------------------------------------

from scipy.integrate import solve_bvp
import numpy as np

#Input ---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*

layers = [2,2,6] #m
M_Es = (10**3)*[45,20,60] #MN/m2

E_c = 33.6*10**6

F = 1000 #kN
D = 1 #m
nodes = 40

#---------------------------------------------------

#total length of the pile:
L = sum(layers) #m

#x-coordinates of the pile
x = np.linspace(0,L,nodes)

#Empty arr. for soil's compression modulus:
M_E = np.zeros(nodes)

#soil layer boundaries' depths:
depth = np.zeros(len(layers)+1)

for i in range(1,len(depth)):
depth[i] = sum(layers[:i])

i=0

#assign M_E to the nodes of the current soil layer "j"
#while x lies between two boundaries assign the i-th value in "M_Es" array
for j in range(len(depth)-1):
while x[i] >= sum_layers[j] and x[i] < sum_layers[j+1]:
M_E[i] = M_Es[j]
i += 1

#calculate the spring stiffness of the soil "k_sh" and the flexural rigidity of the pile "EI"
k_sh = 1.4*M_E/D #kN/m3
EI = E_c*0.25*np.pi*(D*0.5)**4 #kNm2

c = k_sh*D/(EI)

def fun(x, y):
return np.vstack((y[1],y[2],y[3],np.multiply(-c,y[0])))

def bc(ya, yb):
return np.array([ya[2], yb[2], ya[3]+F/EI, yb[3]])

y_a = np.zeros((4, x.size))

res_a = solve_bvp(fun, bc, x, y_a)