Python3 Оде решатель с ненулевыми начальными условиями - PullRequest
Новая
       22

Python3 Оде решатель с ненулевыми начальными условиями

0 голосов
/ 02 мая 2020

Я пытаюсь смоделировать систему массовых пружин с единой степенью свободы. Для временной интеграции я использую функцию ode от scipy. Также я сравниваю цифры с аналитическим решением.

Запуск моего скрипта приводит к той же частотной характеристике, что и аналитическое решение, но после первого шага по амплитуде наблюдается скачок. Я не могу найти источник ошибки в моем коде. У кого-нибудь есть идея?

enter image description here

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()

Ответы [ 2 ]

0 голосов
/ 03 мая 2020

Ваша проблема в том, что вы пытаетесь повторно использовать вектор состояния. Теперь, если интегратор делает что-то похожее на этом простом шаге Эйлера, 

y += f(t,y)*dt

, тогда изменение y в качестве побочного эффекта вызова f резко изменит результат.

Если вы предполагаете, что у вас возможно несколько экземпляров системы, как указано в первой строке 

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

, то вам следует позаботиться о выравнивании результирующего вектора производных 

    return np.concatenate([v,a])
0 голосов
/ 03 мая 2020

Я нашел решение моей проблемы.

Когда функция MBD возвращает список вместо массива numpy, она работает: 

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

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

    a = -K/M*u

    return [v, a]
...