Hai

Di sini kami menggambarkan karya bidang tertentu dan kemudian membuat beberapa fitur yang indah (semuanya SANGAT sederhana di sini).

Apa yang akan terjadi pada artikel ini?

Kasus umum:

Kami menggambarkan dasar, yaitu, bekerja dengan vektor (sepeda untuk mereka yang tidak memiliki numpy di tangan)
Kami menggambarkan titik materi dan bidang interaksi

Kasus khusus (berdasarkan umum):

Mari kita membuat visualisasi bidang vektor kekuatan medan elektromagnetik (gambar pertama dan ketiga)
Kami akan memvisualisasikan gerakan partikel dalam medan elektromagnetik

Temui aku di bawah luka!

Pemrograman Latar Belakang Teoritis

Vektor

Dasar dari semua fondasi adalah vektor (terutama dalam kasus kami). Oleh karena itu, dengan deskripsi vektor kita akan mulai. Apa yang kita butuhkan Operasi aritmatika pada vektor, jarak, modul, dan beberapa hal teknis. Vektor yang akan kita warisi dari daftar. Beginilah tampilannya diinisialisasi:

class Vector(list): def __init__(self, *el): for e in el: self.append(e)

Yaitu, sekarang kita dapat membuat vektor

 v = Vector(1, 2, 3)

Mari kita atur penambahan operasi aritmatika:

 class Vector(list): ... def __add__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] + other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self + other

Hebat:

 v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 + v2 >>> [3, 59, 26.2]

Kami juga mendefinisikan semua operasi aritmatika (kode lengkap vektor akan lebih rendah). Sekarang Anda membutuhkan fungsi jarak. Saya dapat membuat dist rustic (v1, v2) - tetapi ini tidak indah, jadi definisikan ulang% operator:

 class Vector(list): ... def __mod__(self, other): return sum((self - other) ** 2) ** 0.5

Hebat:

 v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 % v2 >>> 58.60068258988115

Kami juga membutuhkan beberapa metode untuk menghasilkan vektor yang lebih cepat dan hasil yang indah. Tidak ada yang rumit di sini, jadi inilah seluruh kode dari kelas Vector:

Semua Kode Kelas Vektor

 class Vector(list): def __init__(self, *el): for e in el: self.append(e) def __add__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] + other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self + other def __sub__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] - other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self - other def __mul__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] * other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self * other def __truediv__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] / other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self / other def __pow__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] ** other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self ** other def __mod__(self, other): return sum((self - other) ** 2) ** 0.5 def mod(self): return self % Vector.emptyvec(len(self)) def dim(self): return len(self) def __str__(self): if len(self) == 0: return "Empty" r = [str(i) for i in self] return "< " + " ".join(r) + " >" def _ipython_display_(self): print(str(self)) @staticmethod def emptyvec(lens=2, n=0): return Vector(*[n for i in range(lens)]) @staticmethod def randvec(dim): return Vector(*[random.random() for i in range(dim)])

Titik material

Di sini, secara teori, semuanya sederhana - titik memiliki koordinat, kecepatan, dan akselerasi (untuk kesederhanaan). Selain itu, ia memiliki massa dan serangkaian parameter khusus (misalnya, untuk medan elektromagnetik, muatan tidak melukai kami, tetapi tidak ada yang mengganggu Anda untuk menyetel setidaknya putaran).

Inisialisasi akan sebagai berikut:

 class Point: def __init__(self, coords, mass=1.0, q=1.0 speed=None, **properties): self.coords = coords if speed is None: self.speed = Vector(*[0 for i in range(len(coords))]) else: self.speed = speed self.acc = Vector(*[0 for i in range(len(coords))]) self.mass = mass self.__params__ = ["coords", "speed", "acc", "q"] + list(properties.keys()) self.q = q for prop in properties: setattr(self, prop, properties[prop])

Dan untuk memindahkan, melumpuhkan dan mempercepat poin kami, kami akan menulis metode berikut:

 class Point: ... def move(self, dt): self.coords = self.coords + self.speed * dt def accelerate(self, dt): self.speed = self.speed + self.acc * dt def accinc(self, force): #        self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0

Bagus sekali, intinya sendiri sudah selesai.

Kode Poin (dengan output yang bagus)

 class Point: def __init__(self, coords, mass=1.0, q=1.0 speed=None, **properties): self.coords = coords if speed is None: self.speed = Vector(*[0 for i in range(len(coords))]) else: self.speed = speed self.acc = Vector(*[0 for i in range(len(coords))]) self.mass = mass self.__params__ = ["coords", "speed", "acc", "q"] + list(properties.keys()) self.q = q for prop in properties: setattr(self, prop, properties[prop]) def move(self, dt): self.coords = self.coords + self.speed * dt def accelerate(self, dt): self.speed = self.speed + self.acc * dt def accinc(self, force): self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0 def __str__(self): r = ["Point {"] for p in self.__params__: r.append(" " + p + " = " + str(getattr(self, p))) r += ["}"] return "\n".join(r) def _ipython_display_(self): print(str(self))

Hasil:

Bidang interaksi

Kami menyebut bidang interaksi sebagai objek yang mencakup sekumpulan semua titik material dan mengerahkannya. Kami akan mempertimbangkan kasus khusus dari alam semesta kita yang indah, jadi kita akan memiliki satu interaksi khusus (tentu saja, ini mudah diperluas). Deklarasikan konstruktor dan tambahkan titik:

 class InteractionField: def __init__(self, F): # F -   , F(p1, p2, r), p1, p2 - , r -    self.points = [] self.F = F def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs))

Sekarang bagian yang menyenangkan adalah mendeklarasikan fungsi yang mengembalikan "ketegangan" pada saat itu. Meskipun konsep ini mengacu pada interaksi elektromagnetik, dalam kasus kami ini adalah beberapa vektor abstrak di mana kita akan memindahkan titik. Dalam hal ini, kita akan memiliki properti dari titik q, dalam kasus tertentu - muatan titik (secara umum - apa pun yang kita inginkan, bahkan vektor). Jadi, apa ketegangan di titik C? Sesuatu seperti ini:

v e c E (v e c C) = j u m l a h v e c F_{i} (v e c C)

$\ vec E (\ vec C) = \ jumlah {\ vec F_i (\ vec C)}$

Yaitu, ketegangan di titik

v e c C

$\ vec C$ sama dengan jumlah gaya semua titik material yang bekerja pada beberapa titik satuan.

 class InteractionField: ... def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) #       (   ,       F,      ) for p in self.points: if coord % p.coords < 10 ** (-10): #      ,        ,   ,     (  ) continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) #    -  proj = proj + (coord - p.coords) / d * fmod #  return proj

Pada titik ini, Anda sudah dapat menggambar bidang vektor, tetapi kami akan melakukannya di akhir. Sekarang mari kita mengambil langkah dalam interaksi kita

 class InteractionField: ... def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * pq) p.accelerate(dt) p.move(dt)

Semuanya sederhana di sini. Untuk setiap titik, kami menentukan tegangan pada koordinat ini dan kemudian menentukan kekuatan akhir pada titik material ETU:

v e c F = v e c E * q

$\ vec F = \ vec E * q$

Tentukan fungsi yang hilang.

Semua Kode InteractionField

 class InteractionField: def __init__(self, F): self.points = [] self.F = F def move_all(self, dt): for p in self.points: p.move(dt) def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) for p in self.points: if coord % p.coords < 10 ** (-10): continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj + (coord - p.coords) / d * fmod return proj def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * pq) p.accelerate(dt) p.move(dt) def clean_acc(self): for p in self.points: p.clean_acc() def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs)) def gather_coords(self): return [p.coords for p in self.points]

Kasus khusus. Gerakan partikel dan visualisasi bidang vektor

Jadi kami sampai pada yang paling menarik. Mari kita mulai dengan ...

Memodelkan gerakan partikel dalam medan elektromagnetik

 u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 + 0.1)) for i in range(3): u.append(Vector.randvec(2) * 10, q=random.random() - 0.5)

Sebenarnya, koefisien k harus sama dengan semacam miliaran (9 * 10 ^ (- 9)), tetapi karena akan padam pada saat t -> 0, saya segera memutuskan untuk membuat keduanya angka yang memadai. Oleh karena itu, dalam fisika kita k = 300'000. Dan dengan yang lainnya, saya pikir sudah jelas.

r ** 2 + 0,1

- ini adalah penghindaran pembelahan oleh 0. Kita, tentu saja, dapat dibingungkan, memecahkan sistem perbedaan yang besar, tetapi pertama-tama tidak ada persamaan gerak untuk lebih dari 2 tubuh, dan kedua ini jelas tidak termasuk dalam konsep "artikel untuk pemula"

Selanjutnya, kami menambahkan sepuluh titik (ruang 2 dimensi) dengan koordinat dari 0 hingga 10 di sepanjang setiap sumbu. Juga, kami memberi setiap poin biaya dari -0,25 hingga 0,25. Sekarang buat lingkaran dan gambar poin sesuai dengan koordinatnya (dan jejaknya):

 X, Y = [], [] for i in range(130): u.step(0.0006) xd, yd = zip(*u.gather_coords()) X.extend(xd) Y.extend(yd) plt.figure(figsize=[8, 8]) plt.scatter(X, Y) plt.scatter(*zip(*u.gather_coords()), color="orange") plt.show()

Apa yang seharusnya terjadi:

Bahkan, gambar di sana akan sepenuhnya acak, karena lintasan setiap titik tidak dapat diprediksi pada saat perkembangan mekanika.

Visualisasi bidang vektor

Semuanya sederhana di sini. Kita perlu melalui koordinat dengan beberapa langkah dan menggambar vektor di masing-masing ke arah yang benar.

 fig = plt.figure(figsize=[5, 5]) res = [] STEP = 0.3 for x in np.arange(0, 10, STEP): for y in np.arange(0, 10, STEP): inten = u.intensity(Vector(x, y)) F = inten.mod() inten /= inten.mod() * 4 #     res.append(([x - inten[0] / 2, x + inten[0] / 2], [y - inten[1] / 2, y + inten[1] / 2], F)) for r in res: plt.plot(r[0], r[1], color=(sigm(r[2]), 0.1, 0.8 * (1 - sigm(r[2])))) #        plt.show()

Kira-kira kesimpulan ini seharusnya diperoleh.

Anda dapat memperpanjang vektor sendiri, ganti * 4 dengan * 1.5:

Kami bermain dengan dimensi dan pemodelan

Buat ruang lima dimensi dengan 200 poin dan interaksi yang tidak tergantung pada kuadrat jarak, tetapi pada derajat ke-4.

 u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 4 + 0.1)) for i in range(200): u.append(Vector.randvec(5) * 10, q=random.random() - 0.5)

Sekarang semua koordinat, kecepatan, dll. Didefinisikan dalam lima dimensi. Sekarang mari kita memodelkan sesuatu:

 velmod = 0 velocities = [] for i in range(100): u.step(0.0005) velmod = sum([p.speed.mod() for p in u.points]) #       velocities.append(velmod) plt.plot(velocities) plt.show()

Ini adalah grafik dari jumlah semua kecepatan pada waktu tertentu. Seperti yang Anda lihat, seiring waktu mereka perlahan-lahan berakselerasi.

Ya, itu adalah instruksi singkat tentang cara membuat hal yang begitu sederhana. Tetapi apa yang terjadi jika Anda bermain dengan bunga:

Semua kode dengan demo

 import random class Vector(list): def __init__(self, *el): for e in el: self.append(e) def __add__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] + other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self + other def __sub__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] - other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self - other def __mul__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] * other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self * other def __truediv__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] / other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self / other def __pow__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] ** other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self ** other def __mod__(self, other): return sum((self - other) ** 2) ** 0.5 def mod(self): return self % Vector.emptyvec(len(self)) def dim(self): return len(self) def __str__(self): if len(self) == 0: return "Empty" r = [str(i) for i in self] return "< " + " ".join(r) + " >" def _ipython_display_(self): print(str(self)) @staticmethod def emptyvec(lens=2, n=0): return Vector(*[n for i in range(lens)]) @staticmethod def randvec(dim): return Vector(*[random.random() for i in range(dim)]) class Point: def __init__(self, coords, mass=1.0, q=1.0, speed=None, **properties): self.coords = coords if speed is None: self.speed = Vector(*[0 for i in range(len(coords))]) else: self.speed = speed self.acc = Vector(*[0 for i in range(len(coords))]) self.mass = mass self.__params__ = ["coords", "speed", "acc", "q"] + list(properties.keys()) self.q = q for prop in properties: setattr(self, prop, properties[prop]) def move(self, dt): self.coords = self.coords + self.speed * dt def accelerate(self, dt): self.speed = self.speed + self.acc * dt def accinc(self, force): self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0 def __str__(self): r = ["Point {"] for p in self.__params__: r.append(" " + p + " = " + str(getattr(self, p))) r += ["}"] return "\n".join(r) def _ipython_display_(self): print(str(self)) class InteractionField: def __init__(self, F): self.points = [] self.F = F def move_all(self, dt): for p in self.points: p.move(dt) def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) for p in self.points: if coord % p.coords < 10 ** (-10): continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj + (coord - p.coords) / d * fmod return proj def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * pq) p.accelerate(dt) p.move(dt) def clean_acc(self): for p in self.points: p.clean_acc() def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs)) def gather_coords(self): return [p.coords for p in self.points] #  import matplotlib.pyplot as plt import numpy as np import time #     if False: u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 + 0.1)) for i in range(10): u.append(Vector.randvec(2) * 10, q=(random.random() - 0.5) / 2) X, Y = [], [] for i in range(130): u.step(0.0006) xd, yd = zip(*u.gather_coords()) X.extend(xd) Y.extend(yd) plt.figure(figsize=[8, 8]) plt.scatter(X, Y) plt.scatter(*zip(*u.gather_coords()), color="orange") plt.show() def sigm(x): return 1 / (1 + 1.10 ** (-x/1000)) #    if False: u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 + 0.1)) for i in range(3): u.append(Vector.randvec(2) * 10, q=random.random() - 0.5) fig = plt.figure(figsize=[5, 5]) res = [] STEP = 0.3 for x in np.arange(0, 10, STEP): for y in np.arange(0, 10, STEP): inten = u.intensity(Vector(x, y)) F = inten.mod() inten /= inten.mod() * 1.5 res.append(([x - inten[0] / 2, x + inten[0] / 2], [y - inten[1] / 2, y + inten[1] / 2], F)) for r in res: plt.plot(r[0], r[1], color=(sigm(r[2]), 0.1, 0.8 * (1 - sigm(r[2])))) plt.show() #    5-  if False: u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 4 + 0.1)) for i in range(200): u.append(Vector.randvec(5) * 10, q=random.random() - 0.5) velmod = 0 velocities = [] for i in range(100): u.step(0.0005) velmod = sum([p.speed.mod() for p in u.points]) velocities.append(velmod) plt.plot(velocities) plt.show()

Artikel selanjutnya mungkin tentang pemodelan yang lebih kompleks, dan mungkin persamaan fluida dan Navier-Stokes.
UPD: Artikel yang ditulis oleh rekan saya di sini

Terima kasih kepada MomoDev karena membantu membuat video.

Simulasi dasar interaksi fisik khusus dalam python + matplotlib