Интеграция уравнения свободы множественной степени свободы

Question

Я хотел бы найти отклик простой динамической системы с 2 степенями свободы, определяемой ODE: m ddx(t) + k ddx(t) = f(t). Система возбуждается с силой f1 = f01 * sin(omega*t) в точке 1:

Матрицы массы и жесткости определяются следующим образом:

import numpy as np 

m1 = 10 
m2 = 8 

k1 = 1e6 
k2 = 3e6 
k3 = 7e6 

m = np.diag([m1,m2]) 
k = np.array([[k1 + k2, -k2], 
       [-k2, k2 + k3]]) 

omega = 500 # excitation angular frequency 

массовая матрица обратного:

minv = np.linalg.inv(m) 

функции интеграции:

def cantilever_beam_ode_with_external_excitation(t, initial_state): 
    """ 
    ODE for a cantilever beam with external excitation (force()) 
    EoM: m ddx + k x = f 
    """ 
    # unpack the state vector 
    x = initial_state[0] # vector 
    xd = initial_state[1] # vector, unimportant, since no damping 

    f = force(t) 

    xdd = minv.dot(-k.dot(x) + f) 

    # return the two state derivatives 
    return [xd, xdd] 

def force(t): 
    force = np.zeros(m.shape[0]) 
    f0 = 100 
    force[0] = f0 * np.sin(omega*t) 
    return force 

Параметры интеграции и интеграции:

# initial conditions 
x0 = np.zeros(m.shape[0]) 
dx0 = np.zeros(m.shape[0]) 
initial_state = (x0, dx0) 

integrator = 'vode' 
state_ode_f = integrate.ode(cantilever_beam_ode_with_external_excitation) 
state_ode_f.set_integrator(integrator, nsteps=500, rtol=1e-3, atol=1e-5, first_step=1e-2) 

t_start = 0 
t_final = 1 

state_ode_f.set_initial_value(initial_state, t_start) 

time = np.array([t_start]) 
velocity = np.expand_dims(initial_state[0], axis=0) 
acceleration = np.expand_dims(initial_state[1], axis=0) 

print("#\t Time\t\t Timestep \t state_ode_f.successful()") 
i = 1 

while state_ode_f.t < (t_final): 
    state_ode_f.integrate(t_final, step=True) 
    time = np.append(time, state_ode_f.t) 
    velocity = np.append(velocity, np.expand_dims(state_ode_f.y[0], axis=0), axis=0) 
    acceleration = np.append(acceleration, np.expand_dims(state_ode_f.y[1], axis = 0), axis=0) 
    print("{0} \t {1:0.5f}\t {2:0.4e} \t {3}".format(i, state_ode_f.t, time[-1] - time[-2], state_ode_f.successful())) 
    i += 1 

я получаю следующее сообщение в результате:

C:\Anaconda3\python.exe D:/faks/LADISK/python/CMS_sparse/integrators/FEM_simple.py 
C:\Anaconda3\lib\site-packages\scipy\integrate\_ode.py:853: UserWarning: vode: Illegal input detected. (See printed message.) 
# Time  Timestep state_ode_f.successful() 
    'Unexpected istate=%s' % istate)) 
4500  0.00000  0.0000e+00  False 
DVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2) 
     In above message, I1 =  116 I2 =  52 
4501  0.00000  0.0000e+00  False 
DVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2) 
     In above message, I1 =  116 I2 =  52 
4502  0.00000  0.0000e+00  False 
DVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2) 
     In above message, I1 =  116 I2 =  52 
4503  0.00000  0.0000e+00  False 
DVODE-- RWORK length needed, LENRW (=I1), exceeds LRW (=I2) 
     In above message, I1 =  116 I2 =  52 
... 

Любые идеи о том, как решить эту проблему приветствуются.

Спасибо!

Answer 1

Я думаю, проблема в том, что вы определяете состояние как список (или кортеж) массивов (каждый 3 элемента). С другой стороны, ode ожидает скаляр, массив или список (скаляров).

я могу заставить его работать (но не могу гарантировать, что правильные значения), повернув состояние в 1 больше массива:

state_array = np.hstack(initial_state) 

И обернуть функцию с тем, которая перестраивает и hstacks результатов

def fn(t, state): 
    # try raveled state 
    state = state.reshape((2,-1)) 
    dd = cantilever_beam_ode_with_external_excitation(t, state) 
    return np.hstack(dd) 

state_ode_f = integrate.ode(fn) 
state_ode_f.set_integrator(integrator, nsteps=500, rtol=1e-3, atol=1e-5, first_step=1e-2) 

t_start = 0 
t_final = 1 
delta_t = 0.1 

state_ode_f.set_initial_value(state_array, t_start) 

for i in range(20): 
    state_ode_f.integrate(t_final, step=True) 
    print(i,state_ode_f.t) 
    print(state_ode_f.y.reshape((2,-1))) 

производит значения, как:

(0, 0.001) 
[[ 0.00000000e+00 0.00000000e+00 -1.42857119e-08] 
[ 0.00000000e+00 0.00000000e+00 -1.42857119e-05]] 
(1, 0.00125) 
[[ 0.00000000e+00 9.76562337e-10 -2.00892823e-08] 
[ 0.00000000e+00 4.10156182e-06 -2.34374957e-05]] 
(2, 0.0015) 
[[ 0.00000000e+00 3.02734325e-09 -2.82366021e-08] 
[ -1.46484351e-07 1.03759748e-05 -3.57561308e-05]] 
(3, 0.00175) 
[[ -1.31835916e-10 8.00781113e-09 -4.13295119e-08] 
[ -6.93054082e-07 2.14355430e-05 -5.38674361e-05]] 
... 
(19, 0.006110558142273646) 
[[ -2.41150395e-06 1.60094544e-05 -1.55746906e-05] 
[ -3.78656683e-03 2.36296597e-02 -2.20818279e-02]]