Блог

Главная - Вопросы по программированию - Runge Kutta4 на python для сложной системы

Runge Kutta4 на python для сложной системы

#python

#python

Вопрос:

Я пытаюсь решить сложную систему, в которой вы можете визуализировать ее как несколько точек / узлов, где между этими узлами подключена пружинно-демпферная система, каждая точка переносит силы от всех других подключенных пружин и амортизаторов в дополнение к гравитационным силам на них, поскольку каждая пружина и амортизатор имеют определенную массу. Я использую классы для создания сетки и начальных условий, но я не уверен, как точно рассчитать новые положения и ускорения, используя runge kutta 4 в этой части определяется runge kutta

 rows = 5
columns = 6

class Runga_kutta4():
def __init__(self, node, u0, v0, t):
    
    self.u0 = u0
    self.v0 = v0
    self.t = t
    self.u = u = 0, 0
    
    self.ux = u[0]
    self.uy = u[1]

    self.v = v = 0, 0

    self.vx = v[0]
    self.vy = v[1]

    f = Forces(u0, u, v0, v)
    self.Node_Forces = f.nodeforces(node)

    self.dt = t[1] - t[0]
    results = self.calculation()
    return results

# Returns the acceleration a
def acceleration(self, Node_Forces):
    """
    F = m *a
    a = F/m
    F_sys = F_externe - (F_damping   Springs) - F_g
    """
    a_list = []
    for (f, m) in zip(Node_Forces, Masses.Lattice_Mass()):
        ax = f[0]/m[0]
        ay = f[1]/m[1]
        a_list.append((ax, ay))
    return a_list.reshape(5, 6)

def calculation(self):

    for i in range(self.t.size - 1):
        # F at time step t / 2
        f_t_05_x = (self.Node_Forces[0][i   1] - self.Node_Forces[0][i]) / 2   self.Node_Forces[0][i]
        f_t_05_y = (self.Node_Forces[1][i   1] - self.Node_Forces[1][i]) / 2   self.Node_Forces[1][i]

        u1x = self.ux[i]
        v1x = self.vx[i]
                    
        u1y = self.uy[i]
        v1y = self.vy[i]
        
        a1x = self.acceleration(self.Node_Forces[0][i])
        a1y = self.acceleration(self.Node_Forces[1][i])
        
        u2x = self.ux[i]   v1x * self.dt / 2
        v2x = self.vx[i]   a1x * self.dt / 2
        
        u2y = self.uy[i]   v1y * self.dt / 2
        v2y = self.vy[i]   a1y * self.dt / 2

        a2x = self.acceleration(f_t_05_x)
        a2y = self.acceleration(f_t_05_y)
        
        u3x = self.ux[i]   v2x * self.dt / 2
        v3x = self.vx[i]   a2x * self.dt / 2
        
        u3y = self.uy[i]   v2y * self.dt / 2
        v3y = self.vy[i]   a2y * self.dt / 2

        a3x = self.acceleration(f_t_05_x)
        a3y = self.acceleration(f_t_05_y)

        u4x = self.ux[i]   v3x * self.dt
        v4x = self.vx[i]   a3x * self.dt

        u4y = self.uy[i]   v3y * self.dt
        v4y = self.vy[i]   a3y * self.dt

        a4x = self.acceleration(self.Node_Forces[0][i   1])
        a4y = self.acceleration(self.Node_Forces[1][i   1])
        
        self.ux[i   1] = self.ux[i]   self.dt / 6 * (v1x   2 * v2x   2 * v3x   v4x)
        self.vx[i   1] = self.vx[i]   self.dt / 6 * (a1x   2 * a2x   2 * a3x   a4x)
        self.uy[i   1] = self.uy[i]   self.dt / 6 * (v1y   2 * v2y   2 * v3y   v4y)
        self.vy[i   1] = self.vy[i]   self.dt / 6 * (a1y   2 * a2y   2 * a3y   a4y)

        self.u = (self.ux, self.uy)
        self.v = (self.vx, self.vy)

    return self.u, self.v


l = Lattice(3)
t0, te, dt = 0, 3, 0.001                            # s
t = np.linspace(t0, te, round((te-t0)/dt   1))

for node in l.latticeNodes():
    position0 = 0, 0
    velocity0 = 0, 0
    state0 = np.append(position0, velocity0)
    new_state = Runga_kutta4(node, position0, velocity0, t)

visualise(l)

фото системы

Комментарии:

1. Существует библиотека с именем scipy (которая довольно стандартная), которая решает ODE и имеет некоторые методы рунге-кутты. В частности, я бы полагался на их реализацию runge kutta 4 в качестве отправной точки для сравнения. Вот рабочий пример решателя ODE от scipy.

2. я не хочу использовать эту библиотеку, поскольку я работаю с разными измерениями и разными уравнениями, поэтому я написал свою собственную функцию, но она не работает

3. Хорошо, круто — можете ли вы записать систему дифференциальных уравнений в latex и добавить изображение уравнений к вашему вопросу? Это значительно облегчило бы понимание того, что вы пытаетесь сделать. Если система уравнений написана четко, гораздо более вероятно, что кто-то с разным опытом набросится на вас с нужной библиотекой / советом / codefix, чтобы помочь вам.

4. итак, теперь я пытаюсь использовать другой подход, создавая метод runge kutta только для одного узла вместо сетки узлов, поэтому мой код теперь выглядит так

5. я отредактировал основной код в вопросе

Метки: