初始条件非零的Python3 ode求解器

Question

我正在尝试模拟具有单个自由度的质量弹簧系统。对于时间积分，我使用scipy中的ode函数。另外，我正在将数值与解析解进行比较。

运行我的脚本所得到的频率响应与分析解决方案相同，但是在第一步之后，幅度出现了跳跃。我无法在代码中找到错误源。有人有主意吗？

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode


def main():

    # spring stiffness, mass, initial displacement
    K, M, u0 = 1.

    v0 = 0.

    # max time and time step size
    T = 10.
    dt = 1E-5

    u, v, t = ODE(K,M,T,dt,u0,v0)

    # analytical solution
    u_analytical = u0 * np.cos(np.sqrt(K/M)*t)

    plt.plot(t,u_analytical)
    plt.plot(t,u)
    plt.legend(['analytic', 'numeric'])
    plt.show()



def MSD(t,states,K,M):

    u, v = np.reshape(states, (2, -1))

    a = -K/M*u

    states[0] = v
    states[1] = a

    return states


def ODE(K,M,T,dt,u0,v0):

    states0 = np.array([u0, v0])

    r = ode(MSD)

    r.set_integrator('vode', method='bdf', order=5, nsteps=1000, rtol=1E-8)

    r.set_initial_value(states0, 0)
    r.set_f_params(K,M)

    time_array = np.linspace(0, T, T/dt + 1)
    u = [u0]
    v = [v0]

    for t in range(1, len(time_array)):

        states = r.integrate(time_array[t])

        u.append(states[0])
        v.append(states[1])

        assert r.successful()

    u = np.asarray(u)
    v = np.asarray(v)

    return u, v, time_array

main()

Answer 1

我找到了解决问题的方法。

当函数MBD返回列表而不是numpy数组时，它将起作用：

def MSD_fixed(t,states,K,M):

    u, v = np.reshape(states, (2, -1))

    a = -K/M*u

    return [v, a]