積分多自由度運動方程

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]]