Оцените скорость на пружине с помощью итеративного подхода

Question

Проблема:

Рассмотрим систему с массой и пружиной, как показано на рисунке ниже . Жесткость пружины и масса объекта известны. Поэтому, если пружина растянута, сила, которую оказывает пружина, может быть рассчитана по закону Гука , а мгновенное ускорение может быть оценено по Ньютону aws движения. Интегрирование ускорения в два раза дает расстояние, на которое будет двигаться пружина, и вычитание этого значения из начальной длины приводит к новому положению для расчета ускорения и повторного запуска l oop. Поэтому, когда ускорение линейно уменьшается, скорость выравнивается при определенном значении (вверху справа) Все, что после этой точки, сжатие и замедление пружины в этом случае не учитывается.

Мой вопрос, как бы на go о кодировании в python. Пока что я написал какой-то псевдокод.

instantaneous_acceleration = lambda x: 5*x/10    # a = kx/m
delta_time = 0.01     #10 milliseconds
a[0] = instantaneous_acceleration(12)     #initial acceleration when stretched to 12 m
v[0] = 0        #initial velocity 0 m/s
s[0] = 12       #initial length 12 m
i = 1
while a[i] > 12:
    v[i] = a[i-1]*delta_time + v[i-1]      #calculate the next velocity
    s[i] = v[i]*delta_time + s[i-1]        #calculate the next position
    a[i] = instantaneous_acceleration (s[i])          #use the position to derive the new accleration
    i = i + 1

Любая помощь или советы приветствуются.

Answer 1

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

x'' = -kx/m
x'' + (k/m)x = 0
r^2 + k/m = 0
r^2 = -(k/m)
r = i*sqrt(k/m)
x(t) = A*e^(i*sqrt(k/m)t)
     = A*cos(sqrt(k/m)t + B) + i*A*sin(sqrt(k/m)t + B)
     = A*cos(sqrt(k/m)t + B)

Из начальных условий мы знаем, что

x(0) = 12 = A*cos(B)
v(0) = 0 = -sqrt(k/m)*A*sin(B)

Второе из этих уравнений верно, только если мы выберем A = 0 или B = 0 или B = Pi.

• если A = 0, то первое уравнение не имеет решения.
• если B = 0, первое уравнение имеет решение A = 12.
• если B = Pi, первое уравнение имеет решение A = -12.

Мы, вероятно, предпочитаем B = 0 и A = 12. Это дает

x(t) =  12*cos(sqrt(k/m)t)
v(t) = -12*sqrt(k/m)*sin(sqrt(k/m)t)
a(t) = -12*(k/m)cos(sqrt(k/m)t)

Таким образом, в любое дополнительное время t [n + 1] = t [n] + dt, мы можем просто вычислить точное положение, скорость и ускорение для t [n] без какого-либо накопления дрейфа или неточности.

Все это говорит, если вы заинтересованы в том, как численно найти x (t) и v (t) и a (t) при произвольном обыкновенном дифференциальном уравнении, ответ гораздо сложнее. Есть много хороших способов сделать то, что можно назвать численной интеграцией. Метод Эйлера самый простой:

// initial conditions

        t[0] = 0
        x[0] = …
       x'[0] = …
         …
  x^(n-1)[0] = …
    x^(n)[0] = 0


// iterative step

  x^(n)[k+1] = f(x^(n-1)[k], …, x'[k], x[k], t[k])
x^(n-1)[k+1] = x^(n-1)[k] + dt * x^(n)[k]
             …
     x'[k+1] = x'[k] + dt * x''[k]
      x[k+1] = x[k] + dt * x'[k]
      t[k+1] = t[k] + dt

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

Более точная версия этого, как правило Симпсона, делает то же самое, но принимает среднее значение над последний квант времени (а не значение любой конечной точки; в приведенном выше примере используется начало интервала). Гарантируется, что среднее значение за интервал будет ближе к истинному значению за интервал, чем любая конечная точка (если только функция не была постоянной в течение этого интервала, и в этом случае Симпсон по крайней мере так же хорош).

Вероятно, Лучшими стандартными методами численного интегрирования для ODE (при условии, что для большей стабильности вам не нужны какие-то методы типа leapfrog) являются методы Runge Kutta. Адаптивный метод временного шага Рунге-Кутты с достаточным порядком обычно делает свое дело и дает вам точные ответы. К сожалению, математика для объяснения методов Рунге-Кутты, вероятно, слишком сложна и трудоемка, чтобы освещать ее здесь, но вы можете найти информацию об этих и других передовых методах в Интернете или, например, в «Числовых рецептах», серии книг по численным методам, которая содержит множество очень полезные примеры кода.

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

Answer 2

У вас есть ошибка знака силы, для пружины или любого другого колебания она всегда должна быть противоположна направлению возбуждения. Исправление этого дает мгновенное колебание. Тем не менее, ваше условие l oop теперь никогда не будет выполнено, поэтому вы также должны адаптировать его.

Вы можете немедленно увеличить порядок вашего метода, повысив его с текущего симплектического c метода Эйлера до Чехарда-Верле. Вам нужно только изменить интерпретацию v[i], чтобы она была скоростью t[i]-dt/2. Затем первое обновление использует ускорение в середине на t[i-1], чтобы вычислить скорость на t[i-1]+dt/2=t[i]-dt/2 от скорости на t[i-1]-dt/2, используя формулу средней точки. Затем в следующей строке обновление позиции представляет собой аналогичную формулу средней точки, использующую скорость в среднее время между временами позиции. Все, что вам нужно изменить в коде, чтобы получить это преимущество, это установить начальную скорость на единицу времени t[0]-dt/2, используя расширение Тейлора на t[0].

instantaneous_acceleration = lambda x: -5*x/10    # a = kx/m
delta_time = 0.01     #10 milliseconds
s0, v0 = 12, 0 #initial length 12 m, initial velocity 0 m/s
N=1000
s = np.zeros(N+1); v = s.copy(); a = s.copy()
a[0] = instantaneous_acceleration(s0)     #initial acceleration when stretched to 12 m
v[0] = v0-a[0]*delta_time/2        
s[0] = s0       
for i in range(N):
    v[i+1] = a[i]*delta_time + v[i]               #calculate the next velocity
    s[i+1] = v[i+1]*delta_time + s[i]             #calculate the next position
    a[i+1] = instantaneous_acceleration (s[i+1])  #use the position to derive the new acceleration
#produce plots of all these functions
t=np.arange(0,N+1)*delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s,v,a)):
    g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();

Это, очевидно, правильно колебание. Точное решение - 12*cos(sqrt(0.5)*t), используя его и его производные для вычисления ошибок в численном решении (помните, что скачкообразное изменение скоростей) дает через

w=0.5**0.5; dt=delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s-12*np.cos(w*t),v+12*w*np.sin(w*(t-dt/2)),a+12*w**2*np.cos(w*t))): 
    g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();

график ниже, показывая ошибки в ожидаемый размер delta_time**2.

Answer 3

Аналитический подход - это самый простой способ получить скорость простой системы, которая подчиняется закону Гука.

Однако, если вы хотите физически точный численный / итерационный подход, я настоятельно рекомендую использовать методы, такие как стандартный Эйлер или методы Рунге-Кутты (предложено Патриком87). [Исправление: метод OPs является симплектическим c методом 1-го порядка, если знак члена ускорения исправлен.]

Вы, вероятно, хотите использовать гамильтонов подход и подходящий интегратор симплект c, такой как чехарда второго порядка (также предложенная Патриком87).

Для закона Хукса вы можете express гамильтониан H = T (p) + V (q), где p - импульс (связанный со скоростью), а q - позиция (связанная с тем, как далеко находится строка из равновесия).

У вас есть кинетическая c энергия T и потенциальная энергия V

T (p) = 0,5 * p ^ 2 / m

V (q) = 0,5 * k * q ^ 2

Вам просто нужны производные от этих двух выражений для моделирования системы

dT / dp = p / m

dV / dq = k * q

Я привел подробный пример (хотя для другой двумерной системы), включая реализацию метода 1-го и 4-го порядка здесь: https://zymplectic.com/case3.html по методу 0 и методу 1

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