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

Я пытаюсь для каждого из 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()

python scipy dictionary ode runge-kutta
person Karina Sakurai    schedule 24.03.2020    source источник
comment
Вы не очищаете список параметров между случаями, как вы думаете, как решатель извлечет текущую силу в момент времени t из списка, переданного в качестве параметра? Тем более, что это задача не решателя, а функции производных. spring ничего не делает в этом отношении.   -  person Lutz Lehmann    schedule 25.03.2020

Ответы (3)


arrow_upward
0
arrow_downward

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

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

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
person Alan    schedule 25.03.2020

arrow_upward
0
arrow_downward

Проблема с возвращаемым значением из 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

to

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

и проверь что получилось..

person DrBwts    schedule 25.03.2020

arrow_upward
0
arrow_downward

Вы совершенно ошибаетесь, передавая массив значений 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

Более питонично использовать присваивания кортежей для распаковки кортежей параметров, выполняя в одной строке то, что ранее было распределено по нескольким, а также делая порядок параметров более наглядным.

Удалите некоторые ненужные сложности, такие как наличие дополнительного механизма подсчета, используйте механизм enumerate(list). Кроме того, не используйте явно созданный список параметров par, когда гораздо проще создать его специально в необязательном аргументе 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-й пример дает хаотично колеблющийся результат.

введите здесь описание изображения

person Lutz Lehmann    schedule 25.03.2020