Как правильно внедрить решатель оды системы принудительной массы с помощью переменной, зависящей от времени, используя scipy.integrate.odeint

Question

Я пытаюсь, для каждого из 5 случаев, численно интегрировать через функцию odeint, функцию массы пружины, причем параметр F изменяется во времени. Тем не менее, он представляет ошибку:

строка 244, в odeint int (bool (tfirst)))

ValueError: установка элемента массива с последовательностью.

from scipy.integrate import odeint as ode
import matplotlib.pyplot as graph
import numpy as np

def spring(y,t,par):

    m = par[0]
    k = par[1]
    c = par[2]
    F = par[3] 

    ydot=[0,0]
    ydot[0] = y[1]
    F = np.array(F)
    ydot[1] = F/m-c*y[1]/m-k*y[0]/m
    return ydot

#cases:[ti,tf,x0,xdot0,c]
A=[0.0,0.3,0.1,0.0,37.7]
B=[0.0,3.0,0.0,1.0,37.7]
C=[0.0,3.0,0.1,1.0,377]
D=[0.0,5.0,0.0,0.0,37.7]
E=[0.0,3.0,0.0,0.0,37.7]

cases = [A, B, C, D, E] #0,1,2,3,4

m=10.0
k=3553.0
par = [m,k]
h = 0.01
cont = 0

for case in cases: 
    x0 = case[2]
    xdot0 = case[3]
    y = [x0,xdot0]
    par.append(case[4])
    ti = case[0]
    tf = case[1]
    t = np.arange(ti, tf,h)

    F = []
    for time in t:
        if cont == 3:
            F.append(1000*np.sin(np.pi*time+np.pi/2))
        elif (cont == 4) and (time >= 0.5): 
            F.append(1000)
        else:
            F.append(0)
    cont = cont + 1
    par.append(F) #F is a vector

    Y = ode(spring,y,t,args=(par,))

    Xn = Y[:,0]
    Vn = Y[:,1]

    graph.plot(t,Xn)
    graph.show()

    graph.plot(t,Vn)
    graph.show()

    graph.plot(Xn,Vn)
    graph.show()

Answer 1

Проблема в том, что возвращаемое значение из spring(y,t,par)

odeint ожидает 2 единичных значения для ydot[0] & ydot[1], которые будут возвращены из spring(), но вместо этого получает одно значение в ydot[0], а затем массив длиной len(F) в ydot[1].

Измените

ydot[1] = F / m - c * y[1] / m - k * y[0] / m

на

ydot[1] = F[0] / m - c * y[1] / m - k * y[0] / m

и проверьте, что произойдет. .

Answer 2

Вы совершенно не правы, передавая массив значений F в производную функцию, не имея возможности связать ее записи со временем t, в котором оценивается функция. Помните, что аргументы y и t являются вектором состояния и (единичным, скалярным) временем, которое решатель должен вычислить на следующем этапе в численном методе. Самый простой способ - передать F как функцию, чтобы можно было напрямую вычислить ее значение.

def spring(y,t,par):

    m,k,c,F = par 
    return [ y[1], F(t)/m-c*y[1]/m-k*y[0]/m ]

#cases:[ti,tf,x0,xdot0,c]
A=[0.0,0.3,0.1,0.0,37.7]
B=[0.0,3.0,0.0,1.0,37.7]
C=[0.0,3.0,0.1,1.0,377]
D=[0.0,5.0,0.0,0.0,37.7]
E=[0.0,3.0,0.0,0.0,37.7]

cases = [A, B, C, D, E] #0,1,2,3,4

m=10.0
k=3553.0
h = 0.01

Для pythoni c удобнее использовать назначения кортежей для распаковки кортежей параметров, выполняя в одной строке то, что ранее было распределено по нескольким, также улучшило видимость порядка параметров.

Удалите некоторые ненужные сложности, такие как наличие дополнительного механизма подсчета, используйте механизм enumerate(list). Кроме того, не используйте явно составленный список параметров par, когда его гораздо проще создать ad-ho c в необязательном аргументе args (это было бы иначе, если бы были сотни систематически сконструированных записей).

В общем случае F может быть функцией интерполяции, сгенерированной методами из scipy.interpolate. Здесь достаточно преобразовать данные уравнения в лямбда-выражения, чтобы определить соответствующие скалярные функции.

for cont,case in enumerate(cases): 
    ti,tf,x0,xdot0,c = case 
    y = [x0,xdot0]
    t = np.arange(ti, tf, h)

    F = lambda t: 0
    if cont == 3:
            F = lambda t: 1000*np.sin(np.pi*t+np.pi/2)
    elif (cont == 4): 
            F = lambda t:  1000 if t >= 0.5 else 0

    Y = ode(spring,y,t,args=([m,k,c,F],))

    Xn,Vn = Y.T

    graph.figure(figsize=(9,3))
    graph.subplot(1,3,1); graph.plot(t,Xn)
    graph.subplot(1,3,2); graph.plot(t,Vn)

    graph.subplot(1,3,3); graph.plot(Xn,Vn)
    graph.tight_layout(); graph.show()

Все остальное происходит так же, как и раньше, с использованием вспомогательных участков для сборки различных графиков в одном изображении. Например, 4-й пример дает хаотически колеблющийся результат

Answer 3

В коде есть другие проблемы, которые возникнут при попытке запуска после устранения этой проблемы. Но в целях точности я отвечу только на ваш вопрос.

Проблема в этой строке кода:

ydot[1] = F/m-c*y[1]/m-k*y[0]/m

F здесь в настоящее время подается как список. В математике, если вы делите вектор на скаляр, предполагается, что вы выполняете поэлементное деление. Python списки не повторяют это поведение, но numpy массивы делают. Так что просто преобразуйте свой список в массив numpy следующим образом:

F = np.array(F)
ydot[1] = F/m-c*y[1]/m-k*y[0]/m