Интегрирование уравнений движения системы с двумя степенями свободы в среде Python

Интегрирование уравнений движения системы с двумя степенями свободы в среде Python

Пример интегрирования в среде Python уравнений движения механической системы с двумя степенями свободы на примере курсовой работы “Исследование движения точки и механической системы” курса теоретической механики Самарского университета

Круглая пластина может вращаться вокруг оси, проходящей через точку \(О_1\) и перпендикулярной плоскости рисунка. В канале пластины движется шарик М, закреплённый на пружине со свободной длиной \(l_0\). Канал, в котором движется шарик, находится на расстоянии 2h = R от оси вращения.

Рассмотрим два режима движения системы:

  • движение системы при известном законе движения пластины;
  • движение системы как системы с двумя степенями свободы.

Система с одной степенью свободы

В первом случае круглая пластина вращается по известному закону \(\varphi(t)\) и механическая система имеет только одну степень свободы – движение шарика относительно пластины.

model.png

Уравнения движения

При движении пластины с известной угловой скоростью \(\dot{\varphi}(t) = \omega\) и ускорением \(\varepsilon\) уравнения относительного движения шарика относительно подвижной системы координат \(Oxy\), связанной с пластиной, имеют следующий вид: \[\left\{ \begin{aligned} & m \ddot x = G \sin \varphi - P - F_{f} + \Phi_e^n \sin \psi + \Phi_e^\tau \cos \psi \\ & m \ddot y = G \cos \varphi + N_y + \Phi_e^n \cos \psi - \Phi_e^\tau \sin \psi - \Phi_c \\ & m \ddot z = N_z \end{aligned} \right.\]

где \(m\) - масса шарика, \(P = c \dot x\) - сила упругости пружины, \(F_f = \text{sign} (V_r) f \sqrt{N_y^2 + N_z^2}\) - сила трения скольжения, \(G = m g\) - сила веса, \(\Phi_e^n\) - переносная центробежная сила инерции: \[\Phi_e^n = m \omega^2 (O_1M) = m \omega^2 \frac{O_1K}{\cos \psi}\]

переносная вращательная сила инерции: \[\Phi_e^\tau = m \varepsilon (O_1M) = m \varepsilon \frac{O_1K}{\cos \psi}\] \[AK = \sqrt{R^2 - h^2}, \quad KM = (x + l_0) - AK, \quad \psi = \arctan \frac{KM}{O_1 K}\]

Движение шарика происходит вдоль оси \(x\), поэтому \(y=0\). Из второго уравнения системы выразим реакцию \(N_y\): \[N_y = \Phi_e^\tau \sin \psi - G \cos \varphi - \Phi_e^n \cos \psi + \Phi_c\]

Дифференциальное уравнение второго порядка, описывающее движения шарика вдоль оси \(Ox\), будет иметь вид: \[\ddot x = g \sin \varphi - \frac{c}{m} x - \text{sign}(V_r) f (\Phi_e^\tau \sin \psi - G \cos \varphi - \Phi_e^n \cos \psi)\frac{1}{m} + \frac{\Phi_e^n}{m} \sin \psi + \frac{\Phi_e^\tau}{m} \cos \psi = f(x,\dot x, t)\]

Приведём это уравнение к системе дифференциальных уравнений первого порядка: \[\left\{ \begin{aligned} & \dot x = V_x \\ & \dot V_x = f(x, V_x, t) \end{aligned} \right.\]

Вектор состояния системы в момент времени t определяется координатным столбцом \[q(t) = \begin{bmatrix} x \\ V_x \end{bmatrix}\]

Дифференциальные уравнения движения в матричной форме: \[\frac{dq}{dt} = f(q, t)\]

Численно проинтегрируем полученные уравнения.

Функция правых частей

Для численного интегрирования уравнений движения системы в среде Python необходимо создать функцию правых частей f(q,t). В начале разрабатываемой программы подключаем модуль для работы с массивами (array):

import numpy as np

и функцию интегрирования дифференциальных уравнений solve_ivp из модуля scipy.integrate:

from scipy.integrate import solve_ivp

Для \(\varphi = \omega t\), \(\omega = \text{const}\) функция правых частей дифференциальных уравнений движения имеет вид:

def dqdt(t, q):    
    x   = q[0]
    vx  = q[1]
    
    # Параметры системы
    g   = 9.807   # ускорение свободного падения
    m   = 0.5     # масса шарика
    w   = 2*np.pi # угловая скорость пластины
    c   = 150     # жёсткость пружины
    L0  = 0.025   # свободная длина пружины
    R   = 0.1     # радиус пластины
    phi = w*t     # Угол поворота пластины    
    f   = 0.0     # Коэффициент трения
        
    # Силы реакции 
    Ny = m*(2*w*vx-g*np.cos(phi)-w*w*R)
    Nz = 0
    # Модуль силы реакции (для вычисления силы трения)
    N  = np.sqrt(Ny**2+Nz**2) 
    
    # Ускорение шарика
    ax  = g*np.sin(w*t) - x*c/m - f*np.sign(vx)*N/m + w*w*(L0+x-np.sqrt(3)*0.5*R)
        
    # аргумент функции - вектор q содержит координату и скорость,
    # возвращаем производную от q, т.е. скорость и ускорение    
    return (vx, ax)

Аргументы функции: t - время, q - вектор состояния системы для момента времени t. Для системы с одной степенью свободы это координата шарика x и его скорость vx.

Интегрирование

Численное интегрирование выполняется при помощи функции solve_ivp

Порядковые аргументы функции:

  1. имя функции правой части системы ДУ
  2. интервал интегрирования от t0 до tN
  3. массив начальных условий

Именованные аргументы:

  • rtol - относительная погрешность
  • atol - абсолютная погрешность

Ошибка на шаге интегрирования не превышает: atol + rtol * abs(y)

По умолчанию используется метод Рунге-Кутты с контролем точности и автоматическим выбором шага.

sol = solve_ivp(dqdt, [0, 5], [0, 0], rtol = 1e-6, atol = 1e-6)

Результат работы функции solve_ivp - структура sol с результатами интегрирования:

    message: The solver successfully reached the end of the integration interval.
    success: True
    status: 0
        t: [ 0.000e+00  1.000e-04 ...  4.992e+00  5.000e+00]
        y: [[ 0.000e+00 -1.215e-08 ...  8.183e-03  1.020e-02]
            [ 0.000e+00 -2.429e-04 ...  2.862e-01  2.480e-01]]
      sol: None
 t_events: None
 y_events: None
     nfev: 1586
     njev: 0
      nlu: 0

sol.t – Массив времени 1xN: первый элемент массива - начальное значение времени t0 последний элемент массива - конечное значение времени интегрирования tN.

sol.y – таблица значений перемещения sol.y[0] и скорости шарика sol.y[1].

Построение графиков

Для построения графиков необходимо подключить библиотеку matplotlib:

import matplotlib.pyplot as plt

Зададим параметры графиков по умолчанию: ширину и высоту графиков в дюймах, размер и тип шрифта

plt.rcParams["figure.figsize"] = (17/2.5,10/2.5)
plt.rcParams["font.size"] = 14
plt.rcParams["font.family"] = 'Arial'

Построим график изменения координаты шарика:

plt.plot(sol.t,sol.y[0])
plt.xlabel('t, c');
plt.ylabel('x, м');
plt.grid(True,linestyle='dotted');

Графики изменения координаты и скорости шарика можно показать на одном рисунке.

plt.plot(sol.t,sol.y[0],'k-', sol.t,sol.y[1],'k--');
plt.xlabel('t, c');
plt.ylabel('м, м/с');
plt.grid(True,linestyle='dotted');

В этом случае необходимо добавить “легенду”

plt.legend(('Координата','Скорость'));


По умолчанию “легенда” располагается на поле графика автоматически. Положение можно указать при помощи строкового параметра loc. Возможные значения:

  • ‘best’
  • ‘upper right’
  • ‘upper left’
  • ‘lower left’
  • ‘lower right’
  • ‘right’
  • ‘center left’
  • ‘center right’
  • ‘lower center’
  • ‘upper center’
  • ‘center’
plt.legend(('Координата','Скорость'), loc = 'upper left');

В подписях осей и легенде можно использовать LaTeX-разметку для математических символов

plt.ylabel('x, м, $V_x$, м/c');
plt.legend(('x','$V_x$'));

Для сохранения рисунка в графический файл используется функция savefig

plt.savefig("x_Vx.png",dpi=300)

Еще один вариант построения графиков двух функций на одном рисунке, используя две оси y с разным масштабом

fig, ax1 = plt.subplots()
line1 = ax1.plot(sol.t,sol.y[0],'k-', label='x');
ax1.set_xlabel('t, c');
ax1.set_ylabel('x, м');

ax2 = ax1.twinx()
line2 = ax2.plot(sol.t,sol.y[1],'k--', label='$V_x$')
ax2.set_ylabel('$V_x$, м/с')

model.png

Передача параметров в функцию правых частей

В приведенном выше примере параметры системы были определены в функции правых частей, что не является хорошей практикой программирования. При таком подходе при изменении параметров системы будет необходимо вносить изменения в эту функцию. Для исследования влияния параметров на движение системы удобнее параметры системы задавать вне тела функции правых частей. Для передачи параметров в функцию правых частей можно использовать объект типа namedtuple из модуля collections.

from collections import namedtuple

# Создаем структуру данных 
Parameters = namedtuple('Parameters', 'g, m, w, c, L0, R, f')

# Создаем набор параметров
p1 = Parameters(9.81,0.5,2*np.pi,150,0.025,0.1,0.0)
# или так
p1 = Parameters(g = 9.81, m = 0.5, w = 2*np.pi, c = 100, L0 = 0.025, R = 0.1, f = 0.0)

# Ещё один набор параметров (с другой жёсткостью пружины)
p2 = Parameters(g = 9.81, m = 0.5, w = 2*np.pi, c = 200, L0 = 0.025, R = 0.1, f = 0.0)

Следует отметить, что как и классический tuple (кортеж), именованный кортеж не изменяем, т.е. после создания нельзя изменить значение атрибута, это приведет к ошибке.

Перепишем функцию правых частей дифференциальных уравнений, добавив третий аргумент функции – структуру с параметрами системы

def dqdt(t, q, p):
    # Функция правых частей
    # t - текущее время
    # q - вектор состояния системы
    # p - параметры системы
    # для системы с одной степенью свободы это координата шарика и его скорость
    x  = q[0]
    vx = q[1]
    
    # Для удобства все параметры извлекаем из именованного кортежа в отдельные переменные
    g, m, w, c, L0, R, f = list(p)

    phi = w*t    # Угол поворота     

    Ny = m*(2*w*vx-g*np.cos(phi)-w*w*R)
    Nz = 0
    N  = np.sqrt(Ny**2+Nz**2)     
    # Вычисляем ускорение шарика
    ax  = g*np.sin(w*t) - x*c/m - f*np.sign(vx)*N/m + w*w*(L0+x-np.sqrt(3)*0.5*R)
   
    # q содержит координату и скорость,
    # возвращаем производную от q - скорость и ускорение    
    return (vx, ax)

Поскольку функция численного интегрирования solve_ivp первым аргументом ожидает функцию с двумя параметрами, используем лямбда-функцию двух аргументов, которая будет вызывать новую функцию правых частей с тремя аргументами, передавая в неё структуру с параметрами системы:

sol1 = solve_ivp(lambda t, q: dqdt(t, q, p1), [0, 5], [0, 0], rtol = 1e-6)

Решение для второго набора параметров

sol2 = solve_ivp(lambda t, q: dqdt(t, q, p2), [0, 5], [0, 0], rtol = 1e-6)

Сравним два решения:

plt.plot(sol1.t, sol1.y[0], 'k-', sol2.t, sol2.y[0],'k--');
plt.legend(['c=100 Н/м','c=200 Н/м']);
plt.ylabel('x, м');
plt.xlabel('t, c');
plt.title('Влияние жёсткости пружины');

Определение силы реакции при движении шарика в канале

Вариант 1.

Создадим функцию, которая будет вычислять реакцию заданного вектора состояния q и параметров p:

def Ny(t, q, p):
    # t - время
    # q - массив из двух элементов (x, Vx)
    g, m, w, c, L0, R, f = list(p)
    return 2*m*w*q[1] - m*g*np.cos(w*t) - m*w*w*R

Склеиваем таблицу времени (строку 1xN) с таблицей решения (матрица 2xN), транспонируем результат. В результате этих действий получим матрицу Nx3, каждая строка которой содержит решение для заданного момента времени:

# t0 x(t0) vx(t0)
# t1 x(t1) vx(t1)
# t2 x(t2) vx(t2)
# ...
# tN x(tN) vx(tN)
data1 = np.vstack((sol1.t,sol1.y)).T

Вызываем функцию Ny(t, q, p) для каждой строки полученной матрицы:

sol1_Ny = list(map( lambda row: Ny(row[0],row[1:],p1), data1 ))

Для второго набора параметров

data2 = np.vstack((sol2.t,sol2.y)).T;
sol2_Ny = list(map( lambda row: Ny(row[0],row[1:],p2), data2 ))

Построим графики сил реакции \(N_y\) для двух наборов параметров

plt.plot(sol1.t,sol1_Ny,'r-',sol2.t,sol2_Ny,'b-.');
plt.xlabel('t, c');
plt.ylabel('$N_y$, Н');
plt.legend(['c=100 Н/м','c=200 Н/м']);

Вариант 2.

Функция, которую можно вызвать сразу для всех значений таблицы решения.

def Ny(t, q, p):
    # t - массив t0, t1, t2, ... tN
    # q - массив [ [ x(t0), x(t1), x(t2), ..., x(tN) ], [ vx(t0), vx(t1), ..., vx(tN) ] ]    
    g, m, w, c, L0, R, f = list(p)
    return 2*m*w*q[1,:] - m*g*np.cos(w*t) - m*R*w**2
plt.plot(sol1.t,Ny(sol1.t,sol1.y,p1), sol2.t,Ny(sol2.t,sol2.y,p2));
plt.xlabel('t, c');
plt.ylabel('$N_y$, Н');
plt.legend(['c=100 Н/м','c=200 Н/м'],loc='lower right');

model.png

Аппроксимация функции sign гладкой функцией

Численные методы интегрирования дифференциальных уравнений обычно предполагают гладкость функции правых частей \(f(t,q)\) по своим аргументам. Аппроксимируем функцию sign, которая использовалась для определения направления действия силы трения, гладкой функцией на основе сигмоиды \[S(x) = \frac{1}{1+e^{-x}}\]

x = np.linspace(-10,10,100)

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))
plt.plot(x, sigmoid(x),'r-', x, np.sign(x), 'k--');
plt.grid()
plt.xlabel('x')
plt.legend(['S(x)','sign(x)']);

“Крутизной” сигмоиды можно управлять при помощи параметра-множителя k перед аргументом x \[S(x) = 2 \left( \frac{1}{1+e^{-k x}} - \frac{1}{2} \right)\]

def sigmoid(x , k = 1):
    return (1.0 / (1.0 + np.exp(-k*x)) - 0.5)*2
x = np.linspace(-0.005,0.005,100)

plt.plot(x,sigmoid(x,10000),x,sigmoid(x,2000),x,sigmoid(x,1000), x, np.sign(x),'k--');
plt.legend(['k=10 000','k=2 000','k=1 000','sign(x)']);

Перепишем функцию правых частей:

def dqdt(t, q, p):
    # Функция правых частей
    # t - текущее время
    # q - вектор состояния системы
    # p - параметры системы
    # для системы с одной степенью свободы это координата шарика и его скорость
    x  = q[0]
    vx = q[1]
    
    # Для удобства все параметры извлекаем из именованного кортежа в отдельные переменные
    g, m, w, c, L0, R, f = list(p)

    phi = w*t    # Угол поворота     

    Ny = m*(2*w*vx-g*np.cos(phi)-w*w*R)
    Nz = 0
    N  = np.sqrt(Ny**2+Nz**2)     
    # Вычисляем ускорение шарика
    ax  = g*np.sin(w*t) - x*c/m - f*sigmoid(vx,10000)*N/m + w*w*(L0+x-np.sqrt(3)*0.5*R)
   
    # q содержит координату и скорость,
    # возвращаем производную от q - скорость и ускорение    
    return (vx, ax)

Решение и графики

p1 = Parameters(g = 9.81, m = 0.5, w = 2*np.pi, c = 100, L0 = 0.025, R = 0.1, f = 0.1)
sol = solve_ivp(lambda t, q: dqdt(t, q, p1), [0, 5], [0, 0], rtol = 1e-7, method = 'BDF')

plt.plot(sol.t,sol.y[0],'k-',sol.t,sol.y[1],'k--');
plt.xlabel('t, c');
plt.ylabel('x, м, $V_x$, м/c');
plt.legend(('x','$V_x$'));

Система с двумя степенями свободы

Во втором рассматриваемом случае механическая системы имеет две степени свободы, при этом движение происходит под действием только потенциальных сил – силы упругости пружины и силы тяжести. Трение при движении шарика внутри канала и при вращении пластины вокруг оси \(O_1\) не учитывается.

model.png

Уравнения движения запишем в форме уравнений Лагранжа второго рода для угла поворота пластины и координаты шарика относительно пластины: \[\begin{aligned} & \frac{d}{dt} \frac{\partial T}{\partial \dot \varphi} - \frac{\partial T}{\partial \varphi} = - \frac{\partial \Pi}{\partial \varphi} \\ & \frac{d}{dt} \frac{\partial T}{\partial \dot x} - \frac{\partial T}{\partial x} = - \frac{\partial \Pi}{\partial x} \end{aligned}\]

Кинетическая энергия системы имеет вид: \[T = \frac{1}{2} m_2\left(\dot{\varphi}^2 \left(l_0-\frac{\sqrt{3} R}{2}+x\right)^2+\left(\dot x - R \dot{\varphi}\right)^2\right)+\frac{3}{8} m_1 R^2 \dot{\varphi}^2\]

Потенциальная энергия системы: \[\Pi = \frac{1}{2} \left(c x^2+g m_2 \sin \varphi \left(-2 l_0 + \sqrt{3} R - 2 x\right) - g R (m_1+2 m_2) \cos \varphi\right)\]

Уравнения движения: \[\left\{ \begin{align} & \left[ m_2 \left(l_0-\frac{\sqrt{3} R}{2}+x\right)^2+\frac{3 m_1 R^2}{4}+m_2 R^2 \right] \cdot \ddot{\varphi} - m_2 R \cdot \ddot x = \frac{1}{2} \left(g m_2 \cos \varphi \left(2 l_0-\sqrt{3} R+2 x\right)-g R (m_1+2 m_2) \sin \varphi \right) - 2 m_2 \dot{\varphi} \dot{x} \left(l_0-\frac{\sqrt{3} R}{2}+x\right) \\ & -m_2 R \cdot \ddot \varphi + m_2 \cdot \ddot x = - c x + m_2 g \sin \varphi + m_2 \left( l_0 - \frac{\sqrt{3}}{2}R + x \right) \dot{\varphi}^2 \end{align} \right.\]

Полученную систему дифференциальных уравнений можно записать в матричной форме \[\left\{ \begin{aligned} & a_{11} \cdot \ddot \varphi + a_{12} \cdot \ddot x = b_1 \\ & a_{21} \cdot \ddot \varphi + a_{22} \cdot \ddot x = b_2 \end{aligned} \right.\]

или \[\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \cdot \begin{bmatrix} \ddot \varphi \\ \ddot x \end{bmatrix} = \begin{bmatrix} b_{1} \\ b_{2} \end{bmatrix}\]

где \[\begin{aligned} & a_{11} = \left[ m_2 \left(l_0-\frac{\sqrt{3} R}{2}+x\right)^2+\frac{3 m_1 R^2}{4}+m_2 R^2 \right], \quad a_{12} = - m_2 R \\ & a_{21} = -m_2 R, \quad a_{22} = m_2 \\ & b_1 = \frac{1}{2} \left(g m_2 \cos \varphi \left(2 l_0-\sqrt{3} R+2 x\right)-g R (m_1+2 m_2) \sin \varphi \right) - 2 m_2 \dot{\varphi} \dot{x} \left(l_0-\frac{\sqrt{3} R}{2}+x\right) \\ & b_2 = - c x + m_2 g \sin \varphi + m_2 \left( l_0 - \frac{\sqrt{3}}{2}R + x \right) \dot{\varphi}^2 \end{aligned}\]

Дифференциальные уравнения движения системы записаны в форме системы линейных уравнений относительно ускорений. Это систему можно решить аналитически и привести систему уравнений к форме Коши. Однако, для системы с большим числом степеней свободы, аналитическое решение и приведение системы к форме Коши будет громоздким, поэтому удобнее решение системы линейных уравнений проводить численно в функции правых частей.

Параметры системы

# Создаем структуру данных 
Parameters = namedtuple('Parameters', 'g, m1, m2, c, L0, R')

# Создаем набор параметров
p = Parameters(g=9.81, m1=2, m2=0.5, c=200, L0=0.03, R=0.1)

Функция правых частей

def dqdt(t, q, p):
    # Функция правых частей
    # t - текущее время
    # q - вектор состояния системы
    # p - параметры системы
    # для системы с одной степенью свободы это координата шарика и его скорость    
    phi = q[0]
    w   = q[1]
    x   = q[2]
    vx  = q[3]
    
    # Для удобства все параметры извлекаем из именованного кортежа в отдельные переменные
    g, m1, m2, c, L0, R = list(p)
    
    # Коэффициенты при ускорения в системе линейных уравнений
    # относительно производных
    a11 = 3*m1*R*R/4 + m2*R*R + m2*(L0 - np.sqrt(3)*R/2 + x)**2
    a12 = -m2*R
    
    a21 = -m2*R
    a22 =  m2
    
    # Элементы матрицы правых частей
    b1  = 0.5*g*( m2*(2*L0-np.sqrt(3)*R+2*x)*np.cos(phi) - (m1+2*m2)*R*np.sin(phi) ) - 2*m2*(L0 - 0.5*np.sqrt(3)*R + x)*vx*w
    b2  = -c*x + m2*(L0-0.5*np.sqrt(3)*R+x)*w*w + g*m2*np.sin(phi)
    
    # Составляем матрицы системы линейный уравнений относительно вторых производных
    A   = np.array( [[a11, a12], [a21, a22]] )
    B   = np.array( [b1,b2])

    # Решаем систему линейных уравнений
    d2q = np.linalg.solve(A, B)   
    
    # Получаем угловое ускорение 
    d2phi = d2q[0]
        # и линейное относительное ускорение шарика
    d2x   = d2q[1]
    
    # q содержит координату и скорость,
    # возвращаем производную от q - скорость и ускорение    
    return (w, d2phi, vx, d2x)

Начальные условия

\[\varphi(0) = 0.2, \quad x(0) = 0, \quad \dot{\varphi}(0) = 0.02, \quad \dot{x}(0) = 0.\]
q0  = [0.2, 0, 0.02, 0]

Интегрирование

sol = solve_ivp(lambda t, q: dqdt(t, q, p), [0, 5], q0, rtol = 1e-7)

График изменения угла поворота пластины

plt.plot(sol.t, sol.y[0], 'k-');
plt.xlabel('t, c');
plt.ylabel('$\\varphi$, рад');
plt.grid(True,linestyle='dotted');

График перемещения шарика

plt.plot(sol.t, sol.y[2], 'k-');
plt.xlabel('t, c');
plt.ylabel('$x$, м');
plt.grid(True,linestyle='dotted');

Проверим корректность построенной модели. Рассматриваемая механическая система с двумя степенями свободы является консервативной: при движении этой системы сумма её потенциальной П и кинетической Т энергий не изменяется.

Функция вычисления кинетической энергии системы

def kinetic_energy(t, q, p):
    phi = q[0]
    w   = q[1]
    x   = q[2]
    vx  = q[3]    
    
    g, m1, m2, c, L0, R = list(p)
    
    T = (3.0/8.0)*m1*R*R*w**2 + 0.5*m2*(w*w*(L0-0.5*np.sqrt(3)*R+x)**2+(vx-R*w)**2)
    
    return T

График изменения кинетической энергии

plt.plot(sol.t, kinetic_energy(sol.t,sol.y,p), 'k-');
plt.xlabel('t, c')
plt.ylabel('T, Дж');
plt.grid(True,linestyle='dotted');

Функция вычисления потенциальной энергии системы

def potential_energy(t, q, p):
    phi = q[0]
    w   = q[1]
    x   = q[2]
    vx  = q[3]    
    
    g, m1, m2, c, L0, R = list(p)

    P = 0.5*(x*x*c + m2*g*(-2*L0+np.sqrt(3)*R-2*x)*np.sin(phi) - g*(m1+2*m2)*R*np.cos(phi))
    
    return P

График изменения потенциальной энергии

plt.plot(sol.t, potential_energy(sol.t,sol.y,p), 'k-');
plt.xlabel('t, c')
plt.ylabel('П, Дж');
plt.grid(True,linestyle='dotted');

Полная энергия

E = kinetic_energy(sol.t,sol.y,p) + potential_energy(sol.t,sol.y,p)

plt.plot(sol.t, kinetic_energy(sol.t,sol.y,p), sol.t, potential_energy(sol.t,sol.y,p), sol.t, E);
plt.grid(True,linestyle='dotted');
plt.xlabel('t, c')
plt.ylabel('Энергия, Дж')
plt.legend(['T','П','T+П']);

График иллюстрирует сохранение полной энергии системы, что свидетельствует о корректности построенной модели механической системы с двумя степенями свободны.


© 2023. All rights reserved.

Powered by Hydejack v9.1.6