Wykres ruchu dwóch ciał

0

O ile podstawy czy nawet odczyt plików w miarę ogarnąłem, tak niestety wykresy na razie sprawiają mi trudność.
Zadanie wygląda tak:
**Wyznacz ruch dwóch ciał połączonych sprężynami o sztywności k=10 znajdujących się polu grawitacyjnym o przyśpieszeniu g⃗ =[0,−9.81][ms2]. Oba ciała są wstępnie ustawione w odległości l=10 od siebie i leżą na linii poziomej, a sprężyna w tym momencie jest nienaciągnięta. Prędkość początkowa ciała pierwszego wynosi u1=[0,5][ms] a drugiego $ u_2=[-3,3] []$. Masy obu ciał są równe i wynoszą m=1[kg] lub m=10[kg] (sprawdź jak zmienią się wyniki w obu przypadkach). Wykonaj wykres torów ruchu obu ciał, a następnie utwórz animację, w której przedstawisz ruch ciał w postaci punktów poruszających się po obliczonych torach. **

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation

k=10
l=10
m_1=1
m_10=10
kroki=500
czas_poczatek=0
czas_koniec=10

lokalizacja_a = [0,0]
lokalizacja_b = [-l,0]

predkosc_a=[0,0.5]
predkosc_b=[-3,3]

grawitacja=9.81

tablica=np.linspace(czas_poczatek,czas_koniec,kroki+1)

dt=tablica[1]-tablica[0]

przyspieszenie=np.zeros((kroki+1,2))

I co dalej?
Dane mam, ale jakim wzorem je ogarnąć? Jakąś funkcję trzeba stworzyć?
Może chodzi o coś takiego: https://pl.wikipedia.org/wiki/Algorytm_Verleta

0

Naprawdę nikt nie jest w stanie nic podpowiedzieć?

0

@aPsikus, musisz wykonać symulacje. Dobierz sobie jakieś małe dt, przypisz poczatkowe wartosci i policz w kazdej ramce dla kazdego obiektu (pseudo kod):

x += v * dt
v += a * dt
a = F / m

zauwaz ze x,v,a i F są wektorami

wzory na F sam musisz policzyc, masa powinna byc gdzies podana

druga opcja to policzyc to samo analitycznie i wtedy po prostu podstawic wartosci dla danego t

0

Masa jest podana.
m_1 = 1
m_10 = 10

a co mi da policzenie x i v?

0
import numpy as np
from numpy.linalg import norm
from scipy.integrate import odeint

import matplotlib.pyplot as plt
from matplotlib import animation


m, l, k = [10, 10, 10]
x1 = [0, 0]
x2 = [-l, 0]
u1 = [0, 50]
u2 = [-30, 30]
endTime=10
nSteps=500
nFrames=100

def equation(Y, t):
    '''
    Equations of motion
    dr = x2 - x1
    m a1 = dr k (1 - l/|dr|)
    m a2 = - dr k (1 - l/|dr|)
    '''
    x1 = Y[:2]
    u1 = Y[2:4]
    x2 = Y[4:6]
    u2 = Y[6:]

    dr = x2-x1

    dY = np.zeros(8)

    if x1[1] < 0 or x2[1] <0:
        return dY

    dY[2:4] = dr*k/m*(1-l/norm(dr))
    dY[3] -= 9.81
    dY[6:] = -dr*k/m*(1-l/norm(dr))
    dY[7] -= 9.81
    dY[:2] = u1
    dY[4:6] = u2

    return dY


Y0 = np.array([x1,u1,x2,u2]).reshape(8)


Y = odeint(equation, Y0, np.linspace(0, endTime, nSteps) )

fig = plt.figure()

plt.plot(Y[:,0], Y[:,1])
plt.plot(Y[:,4], Y[:,5])
p1, = plt.plot(Y0[0],Y0[1], '.', markersize=10, color="blue")
p2, = plt.plot(Y0[4],Y0[5], '.', markersize=10, color="green")
spring, = plt.plot(Y0[4],Y0[5], '-', markersize=10, color="red")


plt.axes().set_aspect('equal', 'datalim')
plt.axes().grid(True)



def anim(frame):
    t = int(len(Y)*float(frame)/nFrames)
    x1 = Y[t,:2]
    x2 = Y[t,4:6]

    p1.set_data(x1)
    p2.set_data(x2)
    spring.set_data([x1[0],x2[0]], [x1[1], x2[1]])



a = animation.FuncAnimation(fig, anim, frames=nFrames, interval=10, repeat=False)


# fig.set_size_inches(3, 3, forward=True)
# a.save('/home/user/Documents/dydaktyka/courses/figures/python_inst03/mass10.gif', writer='imagemagick', fps=100)

plt.show()```

Tylko czy ktoś może mi to wytłumaczyć tak jak chłop krowie na rowie?
1. Co to jest Y? Skąd lub na podstawie czego się bierze?
2. Co to jest t w equation?
0

Przede wszystkim warunki zadania nie określają w jakiej płaszczyźnie jest zdefiniowany pierwotny ruch ciał (mamy tylko dwie zmienne).

  • Jeśli jest ona pozioma, to zadanie jest dużo prostsze, bo położenie 'y' można wyliczyć dokładnie ze wzoru na ruch w polu grawitacyjnym dla dowolnego t, a 'x' i 'z' liczyć przybliżeniami dla każdego kroku, jakby ciała przebywały w nieważkości.
  • Jeśli jest pionowa, to wystarczy liczyć 'x' i 'y' przybliżane w kolejnych krokach czasowych, z uwzględnieniem siły grawitacji i siły sprężystości.

Kolejne pytanie brzmi: jaki jest wzór na siłę z jaką działa nasza sprężyna? Sztywność k=10 niewiele (mi przynajmniej) mówi, więc musisz mieć tam jakieś wzory/tabele podające f zależnie od jej sztywności i rozciągnięcia.

Jak będziesz mieć powyższe, to można pomyśleć nad dalszymi wzorami, bo na razie nie wiadomo, co tutaj i jak liczyć.

0

@Freja Draco: A opublikowane przeze mnie rozwiązane (swoją drogą znalezione na GitHub) nie odpowiada na te pytania?
Próbuję wszystko rozkminić samodzielnie i po prostu utknąłem :(

0

@aPsikus: a powyższy kod z GitHub został napisany do twojego zadania?

0

@Freja Draco: tak, ten kod to rozwiązanie zadania.
Choć nie chcę go kopiować a zrozumieć co i jak robi.

Nie do końca go rozumiem.

Z resztą pytania już zadawałem.

0
aPsikus napisał(a):

Tylko czy ktoś może mi to wytłumaczyć tak jak chłop krowie na rowie?

Da się, o ile rozumiesz podstawy wyżej matematyki. Krowa raczej nie.

aPsikus napisał(a):

Co to jest Y?

Zmienna, która przechowuje dla danego przedziału czasu, w tym przypadku np.linspace(0, endTime, nSteps), parametry ruchu.

aPsikus napisał(a):

Co to jest t w equation?

Y i t (czas) są ze sobą powiązane ogólnym równaniem różniczkowym w postaci:

eq.png

W docstringu masz napisane, jaka jest chwilowa siła działająca na każdą z mas:

m a1 = dr k (1 - l/|dr|)
m a2 = - dr k (1 - l/|dr|)

Następnie te równania są rozwiązywane w tej linii:

Y = odeint(equation, Y0, np.linspace(0, endTime, nSteps))

przy zadanych parametrach początkowych Y0 (czyli odległość między punktami, wektory prędkości i położenia). Oczywiście w tym przykładzie t to czas, ale równie dobrze możesz mieć równanie różniczkowe po dx zamiast po dt.

0

Dzięki @Pyxis

Powiedz proszę, czy z treści zadania wynika jak je rozwiązać, czy też to równanie wyszło tylko na podstawie podanego rozwiązania?

Tylko nadal niewiele rozumiem.
Tworzymy tablicę wypełnioną ośmioma zerami.

dY = np.zeros(8)

Po czym część elementów tablicy zastępujemy bądź obliczeniami, grawitacją lub prędkością?
Dlaczego akurat 8? Skąd wiadomo które elementy jak zamienić?

I podstawowe pytanie: dlaczego prędkość poruszania jest razy 10 względem treści zadania?
**Prędkość początkowa ciała pierwszego wynosi u1=[0,5][ms] a drugiego $ u_2=[-3,3]**

u1 = [0, 50]
u2 = [-30, 30]

I jeszce to: odeint z biblioteki scipy. Skąd wiadomo, że trzeba tego użyć?

Jak widzisz rozpisałem sobie dane z zadania, ale już to jakim równaniem to liczyć i jak nie byłem wstanie zrozumieć.
Mogę skopiować gotowca i oddać za zaliczenie. Ale to mnie nie nauczy niczego...

0

Straszny wysyp pytań ostatnio czyżby zaliczenia na studiach :D

0

Oczywiście :)
Albo projekty własne.

A nie wszystko wyjaśniają, albo zamiast na ćwiczenia czas spędzało się gdzie indziej.

0
aPsikus napisał(a):

Powiedz proszę, czy z treści zadania wynika jak je rozwiązać (...)

Tak, wynika, fizyka to nie zgadywanka.

aPsikus napisał(a):

(...) odeint z biblioteki scipy. Skąd wiadomo, że trzeba tego użyć?

Zdaje się, że chyba nie masz podstaw fizyki. Zacznij od tego i kop dalej. Moduł odeint rozwiązuje równania różniczkowe.

aPsikus napisał(a):

Mogę skopiować gotowca i oddać za zaliczenie. Ale to mnie nie nauczy niczego...

Nauczyć musisz się sam. Na tym polegają studia, by drążyć temat samodzielnie. To nie szkoła średnia. Jakby to była jedna linijka, to bym wyjaśnił, ale nie zamierzam robić wykładu z fizyki + matematyki, bo ewidentnie brakuje Ci podstaw.

1 użytkowników online, w tym zalogowanych: 0, gości: 1