Losowe (lub pseudolosowe) spacery to mój nieleczony zajob moje hobby, o którym uparcie od czasu do czasu tu piszę ku własnej rozrywce, a przerażeniu znudzeniu Czytelnika.
Dzisiaj jednak usuniemy z równania element losowy (i nawet pseudolosowy) - i rozwiążemy sobie geometryczne zadanie z elementami spaceru. Celem wpisu jest nie tyle samo zdanie, co zabawa w programowanie. Najpierw będziemy rysować, a potem trochę policzymy.
Zadanie brzmi następująco: startujemy z jakiegoś punktu na płaszczyźnie. Idziemy jeden metr do przodu. Zakręcamy w lewo o X stopni. Idziemy dwa metry do przodu. Zakręcamy w lewo o X stopni. Idziemy trzy metry do przodu. Zakręcamy w lewo o X stopni. Idziemy cztery metry do przodu. Zakręcamy...
Krótko mówiąc: w N-tym kroku przemieszczamy się o N metrów po linii prostej, a potem obracamy się o X stopni w lewo.
W pewnym momencie zauważamy, że po pewnej liczbie kroków udało nam się przeciąć pierwszy (ten jednometrowy) odcinek trasy naszego spaceru.
Ponieważ jesteśmy wybitnie uzdolnieni matematycznie, natychmiast zauważamy, że niechcący trafiliśmy na najmniejszą możliwą wartość kąta zakrętu X, przy której jest to możliwe.
Pytanie: ile wynosi X? Odpowiedzi należy udzielić z dokładnością do - dajmy na to - trzech miejsc dziesiętnych po przecinku. Albo ośmiu. Wsioryba.
Zadanie można rozwiązać co najmniej na dwa sposoby: analitycznie (i wtedy dostaniemy wynik pełen średnio skomplikowanej trygonometrii), albo przez symulację. Jako że moje gruczoły matematyczne są już dość mocno wydudkane, a także, ponieważ obiecałem, że będziemy kodować - w naturalny sposób wybieram bramkę numer dwa.
Mój pomysł jest następujący: zasymulować dużo spacerów o długości, dajmy na to, stu kroków, dla zakresu kąta X od, dajmy na to, 5° do 175°. Każdy spacer zapisać jako osobną klatkę animacji (warto też gdzieś wyświetlić aktualną wartość X na każdej klatce), a następnie poskładać klatki do, za przeproszeniem, kupy - i sobie tę animację obejrzeć. To powinno nam dać odpowiedź z dokładnością do (być może?) jednego stopnia. A potem napisać sobie osobny skrypt, który symuluje to samo, ale w węższym (za to "gęstszym") zakresie kątów.
Zaczniemy zatem od tej prostszej części.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, FFMpegWriter
# Parametry
ile_kroków = 100
kąty = np.arange(5, 176, 1) #Zakręty, od 5 do 175 stopni, co 1 stopień
# Przygotowanie wykresu
wykres, oś = plt.subplots(figsize=(6, 6))
oś.set_aspect('equal') # Równe jednostki na obu osiach
oś.set_xlim(-200, 200) # Zakładamy, że 200 wystarczy
oś.set_ylim(-200, 200) # Co najwyżej część trasy ucieknie poza kadr
trasa, = oś.plot([], [], lw=1)
# Pojedyncza trasa spaceru (chwilowo pusta).
# Przecinek, żeby potem uniknąć indeksowania (trasa[0]) - bo plot() zwraca listę a my potrzebujemy tylko jednego elementu.
# LW: grubość linii.
def generuj_ścieżkę(kąt):
x, y = [0], [0] # x, y: aktualna pozycja
kąt_bieżący = 0
długość_kroku = 1
for _ in range(ile_kroków):
kąt_radiany = np.deg2rad(kąt_bieżący) # Konwersja stopni na radiany
dx = długość_kroku * np.cos(kąt_radiany) # Przesunięcie x i y
dy = długość_kroku * np.sin(kąt_radiany)
x.append(x[-1] + dx) # Zapamiętujemy nową pozycję
y.append(y[-1] + dy)
kąt_bieżący += kąt # Zakręcamy
długość_kroku += 1 # Wydłużamy krok
return x, y # Zwracamy dwie listy
def klatka_animacji(numer_klatki):
# funkcja wywoływana dla każdej klatki animacji
kąt = kąty[numer_klatki]
x, y = generuj_ścieżkę(kąt)
trasa.set_data(x, y)
oś.set_title(f"Kąt: {kąt}°")
# Auto-zoom - nie polecam
# margines = 1
# xmin, xmax = min(x), max(x)
# ymin, ymax = min(y), max(y)
# oś.set_xlim(xmin - margines, xmax + margines)
# oś.set_ylim(ymin - margines, ymax + margines)
return trasa,
# Tworzymy animację
animacja = FuncAnimation(wykres, klatka_animacji, frames=len(kąty), blit=True)
# blit=True przyspiesza rendering, rysując tylko elementy, które się zmieniły od ostatniej klatki
# Idealne do animacji ze stałym tłem
# Zapisz animację
writer = FFMpegWriter(fps=20)
animacja.save("spacery.mp4", writer=writer)
plt.close()
print("Gotowe - sprawdź plik spacery.mp4")
Efekt:
Hmmm. Całkiem ładne, ale szczerze mówiąc gówno widać. Zmniejszamy liczbę kroków do 10 i zawężamy zakres kątów do 90-170. Robimy też zbliżenie:
Ha. Pauzując w okolicy 137 stopni widzimy, że czwarty odcinek dotyka pierwszego.
A więc możemy skrócić trasę do czterech kroków i zrobić jeszcze lepsze zbliżenie. A także nieco zagęścić kąty:
Wyraźnie widać, że styk następuje w okolicach 138 stopni - ale ciężko ustalić dokładniejszy kąt.
Tu przechodzimy do części drugiej, czyli metody numeryczne.
Pomysł jest następujący: napisać funkcję, która na wejściu dostaje kąt X, a na wyjściu zwróci jedynkę jeżeli czwarty odcinek trasy przecina pierwszy, lub minus jedynkę - jeżeli nie przecina. A następnie zapuścić jakąś interpolację kwadratową czy inną bisekcję, żeby znaleźć konkretny kąt z zadaną dokładnością.
from shapely.geometry import LineString
from scipy.optimize import brentq
import numpy as np
def generuj_trasę(kąt_zakrętu, liczba_kroków=4):
punkty = [(0, 0)]
kąt = 0
długość_kroku = 1
for _ in range(liczba_kroków):
kąt_rad = np.deg2rad(kąt)
dx = długość_kroku * np.cos(kąt_rad)
dy = długość_kroku * np.sin(kąt_rad)
nowy_punkt = (punkty[-1][0] + dx, punkty[-1][1] + dy)
punkty.append(nowy_punkt)
kąt += kąt_zakrętu
długość_kroku += 1
return [LineString([punkty[i], punkty[i+1]]) for i in range(len(punkty)-1)]
def przecina_czy_nie(phi):
segments = generuj_trasę(phi)
return 1 if segments[0].intersects(segments[3]) else -1
# szukamy kąta, w którym 2 linie przecinają się
szukany_kąt = brentq(przecina_czy_nie, 130, 140, xtol=1e-5)
# tu się dzieje magia właściwa: brentq szuka miejsca zerowego zadanej funkcji
# w odróżnieniu od tradycyjnej bisekcji, brentq używa interpolacji kwadratowej
# więc jest szybsza a także bezpieczeniejsza
# Parametry:
# 130, 140 - przedział startowy
# xtol: tolerancja (dokładność) - tu 0.00001
# przecina_czy_nie: funkcja, której miejsca zerowego szukamy
print(f"Przecięcie nastąpiło przy kącie ≈ {szukany_kąt:.5f}°")
Wynik:
Przecięcie nastąpiło przy kącie ≈ 138.59038°
Jak dla mnie jest to dokładność całkiem wystarczająca.
Bonusik: można wykazać, że w momencie styku odcinek numer jeden dotyka odcinka numer cztery dokładnie w jego połowie, a więc odcinki numer 2, 3, oraz połowa odcinka nr 4 tworzą trójkąt równoramienny o ramionach długości 2 i podstawie 3. Ale już mi zabrakło rozpędu, żeby to tutaj rozpisać.
Nudne, prawda?
Jeżeli chcesz do komentarza wstawić kod, użyj składni:
[code]
tutaj wstaw swój kod
[/code]
Jeżeli zrobisz literówkę lub zmienisz zdanie, możesz edytować komentarz po jego zatwierdzeniu.