Интеграция Runge kutta для отклика подвеса

Question

Я должен использовать интеграцию Runge Kutta, чтобы решить и построить реакцию ускорения подвески транспортного средства для определенного интервала времени. Когда на графике изображен отклик, он должен иметь несколько колебаний, но из-за демпфирования амплитуда колебаний будет уменьшаться до тех пор, пока отклик ускорения не исчезнет.

Я не могу заставить этот ответ вымереть. Глядя на мой окончательный результат, компонент "z" в массиве продолжает оставаться равным его начальному условию, нулю. Я думаю, что это может быть проблемой, но не знаю, как ее решить. ODEint это то, что я никогда не использовал раньше.

Вот мой код на сегодняшний день:

Не беспокойтесь о том, как я настроил свои матрицы и поработаю с ними (если только здесь не будет ошибки), но они соответствуют моим ручным подсчетам.

    # -*- coding: utf-8 -*-
    """
    Created on Tue Oct 18 11:57:20 2016

    @author: FrancoH
    """
    import numpy as np
    import numpy.linalg as linalg
    from scipy.integrate import odeint
    import matplotlib.pyplot as plt

    m = 1760     #(Suspended mass)
    k_f = 66757  #(Front Spring Coefficient)
    k_r = 154816 #(Rear Spring Coefficient)
    c_f = 3040   #N.s/m
    c_r = 3040
    J = 3737.067 #moment of inertia
    l_f = 1.43   #front length
    l_r = 1.47
    w = 39.969   #omega


    #time array   
    t0 = 0.0                # Start time 
    tf = 3                  # end time 
    intervals = 1000             
    tspan = np.linspace(t0,tf,intervals)


    #matrixes
    c_1 = c_f+c_r
    c_2 = -1*(l_r)*(c_r)+(l_f)*(c_f)
    c_3 = -1*(l_r)*(c_r)+(l_f)*(c_f)
    c_4 = ((l_r)**2)*(c_r)+((l_f)**2)*(c_f)

    k_1 = k_f + k_r
    k_2 = -1*(l_r)*(k_r) + (l_f)*(k_f)
    k_3 = -1*(l_r)*(k_r) + (l_f)*(k_f)
    k_4 = ((l_r)**2)*(k_r)+((l_f)**2)*(k_f)

    m_matrix = np.zeros((2,2))
    c_matrix = np.zeros((2,2))
    k_matrix = np.zeros((2,2))

    for i in range(0,1):

        np.put(m_matrix[i],[i],m)
        np.put(m_matrix[i+1],[i+1],J)

        np.put(c_matrix[i],[i],c_1)
        np.put(c_matrix[i],[i+1],c_2)
        np.put(c_matrix[i+1],[i],c_3)
        np.put(c_matrix[i+1],[i+1],c_4)

        np.put(k_matrix[i],[i],k_1)
        np.put(k_matrix[i],[i+1],k_2)
        np.put(k_matrix[i+1],[i],k_3)
        np.put(k_matrix[i+1],[i+1],k_4)

    #print c_matrix
    #print k_matrix

    m_inv = linalg.inv(m_matrix)  

    E_1_1 = [[0,0],[0,0]]
    E_1_2 = [[J,J],[J,J]]
    E_2_1 = -1*m_inv*k_matrix
    E_2_2 = -1*m_inv*c_matrix

    # to calculate the entries of matrix F

    F_1_1 = -1*k_r
    F_1_2 = -1*c_r
    F_2_1 = k_r*l_r 
    F_2_2 = c_f*l_f

    F_matrix = [[F_1_1,F_1_2],[F_2_1,F_2_2]]
    MF_matrix = m_inv*F_matrix 

    # to calculate entries for Q matrix 

    #y_matrix = [[np.sin(w*tspan)],[w*np.cos(w*tspan)]]

    """ Rearx"""

    def fr(y,t,params):

        z , x = y 
        A,B,C,D,w = params
        derivs = [z,A*x+B*z+C*np.sin(w*t+np.pi)+D*w*np.cos(w*t+np.pi)]
        return derivs

Выше: z - первая производная от x по времени. Функция справа - вторая производная от x по времени, но теперь она переписывается как первая производная от x по времени.

    # initial values
    z0 = 0.0
    x0 = 0.0

    # params
    params = [E_2_1[0,0],E_2_2[0,0],MF_matrix[0][0],MF_matrix[0][1],w]

    # Bundle initial values
    y0 = [z0,x0]

Ниже:
Здесь я вызываю функцию "odeint". "fr" относится к функции выше, где мое уравнение движения второго порядка было преобразовано в первый порядок. Это я сделал вручную.

    # call ode solver
    SolutionRearx = odeint(fr,y0,tspan,args=(params,))

    plt.plot(tspan,SolutionRearx[:,1],'r')

Массив SolutionRearx возвращает все первые значения в виде нулей, и это не является желаемым результатом. Это приводит к тому, что вторые значения также являются ошибочными.

Я не знаю, является ли мой подход к методу Рунге Кутты даже правильным, и, кроме того, я совершенно не умею заниматься кодированием.

Answer 1

Это может быть так же просто, как исправить первую строку в fr Функция ODE для

x, z = y

вместо z,x = y,

На выходе производной функции первый компонент затем является производной x, который zи второй компонент является производной от z которая является второй производной от x который определяется исходным дифференциальным уравнением второго порядка, так что теперь это правильно, как есть.