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

Question

我正在尝试模拟具有单一自由度的质量 spring 系统。对于时间积分，我使用 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 returns 一个列表而不是一个 numpy 数组时它起作用：

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

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

    a = -K/M*u

    return [v, a]

Answer 2

您的问题是您尝试重新使用状态向量。现在，如果积分器在这个简单的欧拉步骤中做类似的事情，

y += f(t,y)*dt

然后更改 y 作为调用 f 的副作用将显着改变结果。

如果您假设您可能有多个系统实例，如第一行所示

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

那么你应该关心扁平化结果导数向量

    return np.concatenate([v,a])

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

Python3 ode solver with non-zero initial conditions failing

ode

python-3.x