Hướng dẫn animation scatter plot python - hoạt hình phân tán âm mưu python
Matplotlib Phiên bản 1.1 đã thêm một số công cụ để tạo hình ảnh động thực sự lắt léo. Bạn có thể tìm thấy một số hình ảnh động ví dụ tốt trên trang ví dụ Matplotlib. Tôi nghĩ rằng tôi sẽ chia sẻ ở đây một số điều tôi đã học được khi chơi xung quanh với các công cụ này. Hoạt hình cơ bảnCác công cụ hoạt hình trung tâm xung quanh lớp cơ sở Trước tiên, chúng tôi sẽ sử dụng """ Matplotlib Animation Example author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ import numpy as np from matplotlib import pyplot as plt from matplotlib import animation # First set up the figure, the axis, and the plot element we want to animate fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2) # initialization function: plot the background of each frame def init(): line.set_data([], []) return line, # animation function. This is called sequentially def animate(i): x = np.linspace(0, 2, 1000) y = np.sin(2 * np.pi * (x - 0.01 * i)) line.set_data(x, y) return line, # call the animator. blit=True means only re-draw the parts that have changed. anim = animation.FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html anim.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show()00000 download """ Matplotlib Animation Example author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ import numpy as np from matplotlib import pyplot as plt from matplotlib import animation # First set up the figure, the axis, and the plot element we want to animate fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2) # initialization function: plot the background of each frame def init(): line.set_data([], []) return line, # animation function. This is called sequentially def animate(i): x = np.linspace(0, 2, 1000) y = np.sin(2 * np.pi * (x - 0.01 * i)) line.set_data(x, y) return line, # call the animator. blit=True means only re-draw the parts that have changed. anim = animation.FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html anim.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show() Hãy bước qua điều này và xem những gì đang xảy ra. Sau khi nhập các phần yêu cầu của fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)2 và fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)3, tập lệnh đã thiết lập cốt truyện: fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2) Ở đây chúng tôi tạo một cửa sổ hình, tạo một trục duy nhất trong hình và sau đó tạo đối tượng dòng của chúng tôi sẽ được sửa đổi trong hình ảnh động. Lưu ý rằng ở đây chúng tôi chỉ cần vẽ một dòng trống: Chúng tôi sẽ thêm dữ liệu vào dòng sau. Tiếp theo, chúng tôi sẽ tạo các chức năng làm cho hoạt hình xảy ra. fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)4 là chức năng sẽ được gọi để tạo khung cơ sở mà hình ảnh động diễn ra. Ở đây chúng tôi chỉ sử dụng một hàm đơn giản đặt dữ liệu dòng thành không có gì. Điều quan trọng là chức năng này trả về đối tượng dòng, bởi vì điều này cho người hoạt hình đối tượng trên lô để cập nhật sau mỗi khung hình: def init(): line.set_data([], []) return line, Phần tiếp theo là chức năng hoạt hình. Nó có một tham số duy nhất, số khung fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)5 và vẽ sóng hình sin với sự thay đổi phụ thuộc vào fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)5: # animation function. This is called sequentially def animate(i): x = np.linspace(0, 2, 1000) y = np.sin(2 * np.pi * (x - 0.01 * i)) line.set_data(x, y) return line, Lưu ý rằng một lần nữa ở đây, chúng tôi trả lại một tuple của các đối tượng cốt truyện đã được sửa đổi. Điều này cho khung hoạt hình những phần của cốt truyện nên được hoạt hình. Cuối cùng, chúng tôi tạo đối tượng hoạt hình: anim = animation.FuncAnimation(fig, animate, init_func=init, frames=100, interval=20, blit=True) Đối tượng này cần phải tồn tại, vì vậy nó phải được gán cho một biến. Chúng tôi đã chọn hoạt hình 100 khung hình với độ trễ 20ms giữa các khung. Từ khóa fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)7 là một từ khóa quan trọng: điều này cho hoạt hình chỉ vẽ lại các phần của cốt truyện đã thay đổi. Thời gian lưu với fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)8 có nghĩa là hình ảnh động hiển thị nhanh hơn nhiều. Chúng tôi kết thúc với một lệnh lưu tùy chọn, và sau đó là lệnh hiển thị kết quả. Đây là những gì tập lệnh tạo ra: Khung này để tạo và lưu hoạt hình rất mạnh mẽ và linh hoạt: nếu chúng ta đưa một số vật lý vào chức năng fig = plt.figure() ax = plt.axes(xlim=(0, 2), ylim=(-2, 2)) line, = ax.plot([], [], lw=2)9, khả năng là vô tận. Dưới đây là một vài ví dụ về một số hình ảnh động vật lý mà tôi đã chơi xung quanh. Con lắc képMột trong những ví dụ được cung cấp trên trang ví dụ Matplotlib là một hình ảnh động của một con lắc kép. Ví dụ này hoạt động bằng cách tính toán vị trí con lắc trong hơn 10 giây, và sau đó hoạt hình kết quả. Tôi đã thấy điều này và tự hỏi liệu Python có đủ nhanh để tính toán động lực học không. Hóa ra nó là: Double Pendulum double_pendulum.py download""" General Numerical Solver for the 1D Time-Dependent Schrodinger's equation. adapted from code at http://matplotlib.sourceforge.net/examples/animation/double_pendulum_animated.py Double pendulum formula translated from the C code at http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ from numpy import sin, cos import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib.animation as animation class DoublePendulum: """Double Pendulum Class init_state is [theta1, omega1, theta2, omega2] in degrees, where theta1, omega1 is the angular position and velocity of the first pendulum arm, and theta2, omega2 is that of the second pendulum arm """ def __init__(self, init_state = [120, 0, -20, 0], L1=1.0, # length of pendulum 1 in m L2=1.0, # length of pendulum 2 in m M1=1.0, # mass of pendulum 1 in kg M2=1.0, # mass of pendulum 2 in kg G=9.8, # acceleration due to gravity, in m/s^2 origin=(0, 0)): self.init_state = np.asarray(init_state, dtype='float') self.params = (L1, L2, M1, M2, G) self.origin = origin self.time_elapsed = 0 self.state = self.init_state * np.pi / 180. def position(self): """compute the current x,y positions of the pendulum arms""" (L1, L2, M1, M2, G) = self.params x = np.cumsum([self.origin[0], L1 * sin(self.state[0]), L2 * sin(self.state[2])]) y = np.cumsum([self.origin[1], -L1 * cos(self.state[0]), -L2 * cos(self.state[2])]) return (x, y) def energy(self): """compute the energy of the current state""" (L1, L2, M1, M2, G) = self.params x = np.cumsum([L1 * sin(self.state[0]), L2 * sin(self.state[2])]) y = np.cumsum([-L1 * cos(self.state[0]), -L2 * cos(self.state[2])]) vx = np.cumsum([L1 * self.state[1] * cos(self.state[0]), L2 * self.state[3] * cos(self.state[2])]) vy = np.cumsum([L1 * self.state[1] * sin(self.state[0]), L2 * self.state[3] * sin(self.state[2])]) U = G * (M1 * y[0] + M2 * y[1]) K = 0.5 * (M1 * np.dot(vx, vx) + M2 * np.dot(vy, vy)) return U + K def dstate_dt(self, state, t): """compute the derivative of the given state""" (M1, M2, L1, L2, G) = self.params dydx = np.zeros_like(state) dydx[0] = state[1] dydx[2] = state[3] cos_delta = cos(state[2] - state[0]) sin_delta = sin(state[2] - state[0]) den1 = (M1 + M2) * L1 - M2 * L1 * cos_delta * cos_delta dydx[1] = (M2 * L1 * state[1] * state[1] * sin_delta * cos_delta + M2 * G * sin(state[2]) * cos_delta + M2 * L2 * state[3] * state[3] * sin_delta - (M1 + M2) * G * sin(state[0])) / den1 den2 = (L2 / L1) * den1 dydx[3] = (-M2 * L2 * state[3] * state[3] * sin_delta * cos_delta + (M1 + M2) * G * sin(state[0]) * cos_delta - (M1 + M2) * L1 * state[1] * state[1] * sin_delta - (M1 + M2) * G * sin(state[2])) / den2 return dydx def step(self, dt): """execute one time step of length dt and update state""" self.state = integrate.odeint(self.dstate_dt, self.state, [0, dt])[1] self.time_elapsed += dt #------------------------------------------------------------ # set up initial state and global variables pendulum = DoublePendulum([180., 0.0, -20., 0.0]) dt = 1./30 # 30 fps #------------------------------------------------------------ # set up figure and animation fig = plt.figure() ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2)) ax.grid() line, = ax.plot([], [], 'o-', lw=2) time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes) energy_text = ax.text(0.02, 0.90, '', transform=ax.transAxes) def init(): """initialize animation""" line.set_data([], []) time_text.set_text('') energy_text.set_text('') return line, time_text, energy_text def animate(i): """perform animation step""" global pendulum, dt pendulum.step(dt) line.set_data(*pendulum.position()) time_text.set_text('time = %.1f' % pendulum.time_elapsed) energy_text.set_text('energy = %.3f J' % pendulum.energy()) return line, time_text, energy_text # choose the interval based on dt and the time to animate one step from time import time t0 = time() animate(0) t1 = time() interval = 1000 * dt - (t1 - t0) ani = animation.FuncAnimation(fig, animate, frames=300, interval=interval, blit=True, init_func=init) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html #ani.save('double_pendulum.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show()download """ General Numerical Solver for the 1D Time-Dependent Schrodinger's equation. adapted from code at http://matplotlib.sourceforge.net/examples/animation/double_pendulum_animated.py Double pendulum formula translated from the C code at http://www.physics.usyd.edu.au/~wheat/dpend_html/solve_dpend.c author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ from numpy import sin, cos import numpy as np import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib.animation as animation class DoublePendulum: """Double Pendulum Class init_state is [theta1, omega1, theta2, omega2] in degrees, where theta1, omega1 is the angular position and velocity of the first pendulum arm, and theta2, omega2 is that of the second pendulum arm """ def __init__(self, init_state = [120, 0, -20, 0], L1=1.0, # length of pendulum 1 in m L2=1.0, # length of pendulum 2 in m M1=1.0, # mass of pendulum 1 in kg M2=1.0, # mass of pendulum 2 in kg G=9.8, # acceleration due to gravity, in m/s^2 origin=(0, 0)): self.init_state = np.asarray(init_state, dtype='float') self.params = (L1, L2, M1, M2, G) self.origin = origin self.time_elapsed = 0 self.state = self.init_state * np.pi / 180. def position(self): """compute the current x,y positions of the pendulum arms""" (L1, L2, M1, M2, G) = self.params x = np.cumsum([self.origin[0], L1 * sin(self.state[0]), L2 * sin(self.state[2])]) y = np.cumsum([self.origin[1], -L1 * cos(self.state[0]), -L2 * cos(self.state[2])]) return (x, y) def energy(self): """compute the energy of the current state""" (L1, L2, M1, M2, G) = self.params x = np.cumsum([L1 * sin(self.state[0]), L2 * sin(self.state[2])]) y = np.cumsum([-L1 * cos(self.state[0]), -L2 * cos(self.state[2])]) vx = np.cumsum([L1 * self.state[1] * cos(self.state[0]), L2 * self.state[3] * cos(self.state[2])]) vy = np.cumsum([L1 * self.state[1] * sin(self.state[0]), L2 * self.state[3] * sin(self.state[2])]) U = G * (M1 * y[0] + M2 * y[1]) K = 0.5 * (M1 * np.dot(vx, vx) + M2 * np.dot(vy, vy)) return U + K def dstate_dt(self, state, t): """compute the derivative of the given state""" (M1, M2, L1, L2, G) = self.params dydx = np.zeros_like(state) dydx[0] = state[1] dydx[2] = state[3] cos_delta = cos(state[2] - state[0]) sin_delta = sin(state[2] - state[0]) den1 = (M1 + M2) * L1 - M2 * L1 * cos_delta * cos_delta dydx[1] = (M2 * L1 * state[1] * state[1] * sin_delta * cos_delta + M2 * G * sin(state[2]) * cos_delta + M2 * L2 * state[3] * state[3] * sin_delta - (M1 + M2) * G * sin(state[0])) / den1 den2 = (L2 / L1) * den1 dydx[3] = (-M2 * L2 * state[3] * state[3] * sin_delta * cos_delta + (M1 + M2) * G * sin(state[0]) * cos_delta - (M1 + M2) * L1 * state[1] * state[1] * sin_delta - (M1 + M2) * G * sin(state[2])) / den2 return dydx def step(self, dt): """execute one time step of length dt and update state""" self.state = integrate.odeint(self.dstate_dt, self.state, [0, dt])[1] self.time_elapsed += dt #------------------------------------------------------------ # set up initial state and global variables pendulum = DoublePendulum([180., 0.0, -20., 0.0]) dt = 1./30 # 30 fps #------------------------------------------------------------ # set up figure and animation fig = plt.figure() ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2)) ax.grid() line, = ax.plot([], [], 'o-', lw=2) time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes) energy_text = ax.text(0.02, 0.90, '', transform=ax.transAxes) def init(): """initialize animation""" line.set_data([], []) time_text.set_text('') energy_text.set_text('') return line, time_text, energy_text def animate(i): """perform animation step""" global pendulum, dt pendulum.step(dt) line.set_data(*pendulum.position()) time_text.set_text('time = %.1f' % pendulum.time_elapsed) energy_text.set_text('energy = %.3f J' % pendulum.energy()) return line, time_text, energy_text # choose the interval based on dt and the time to animate one step from time import time t0 = time() animate(0) t1 = time() interval = 1000 * dt - (t1 - t0) ani = animation.FuncAnimation(fig, animate, frames=300, interval=interval, blit=True, init_func=init) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html #ani.save('double_pendulum.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show() Ở đây chúng tôi đã tạo ra một lớp lưu trữ trạng thái của con lắc kép (được mã hóa theo góc của mỗi cánh tay cộng với vận tốc góc của mỗi cánh tay) và cũng cung cấp một số chức năng để tính toán động lực học. Các hàm hoạt hình giống như trên, nhưng chúng ta chỉ có chức năng cập nhật phức tạp hơn một chút: nó không chỉ thay đổi vị trí của các điểm mà còn thay đổi văn bản để theo dõi thời gian và năng lượng (năng lượng phải không đổi nếu toán học của chúng ta là chính xác: đó là sự an ủi như vậy). Video dưới đây chỉ kéo dài mười giây, nhưng bằng cách chạy tập lệnh, bạn có thể xem con lắc dao động một cách hỗn loạn cho đến khi máy tính xách tay của bạn hết điện: Các hạt trong hộpMột hình ảnh động khác mà tôi tạo ra là sự va chạm đàn hồi của một nhóm các hạt trong một hộp dưới lực hấp dẫn. Các va chạm là đàn hồi: chúng bảo tồn năng lượng và động lượng 2D, và các hạt nảy ra một cách thực tế khỏi các bức tường của hộp và rơi vào ảnh hưởng của một lực hấp dẫn không đổi: """ Animation of Elastic collisions with Gravity author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ import numpy as np from scipy.spatial.distance import pdist, squareform import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib.animation as animation class ParticleBox: """Orbits class init_state is an [N x 4] array, where N is the number of particles: [[x1, y1, vx1, vy1], [x2, y2, vx2, vy2], ... ] bounds is the size of the box: [xmin, xmax, ymin, ymax] """ def __init__(self, init_state = [[1, 0, 0, -1], [-0.5, 0.5, 0.5, 0.5], [-0.5, -0.5, -0.5, 0.5]], bounds = [-2, 2, -2, 2], size = 0.04, M = 0.05, G = 9.8): self.init_state = np.asarray(init_state, dtype=float) self.M = M * np.ones(self.init_state.shape[0]) self.size = size self.state = self.init_state.copy() self.time_elapsed = 0 self.bounds = bounds self.G = G def step(self, dt): """step once by dt seconds""" self.time_elapsed += dt # update positions self.state[:, :2] += dt * self.state[:, 2:] # find pairs of particles undergoing a collision D = squareform(pdist(self.state[:, :2])) ind1, ind2 = np.where(D < 2 * self.size) unique = (ind1 < ind2) ind1 = ind1[unique] ind2 = ind2[unique] # update velocities of colliding pairs for i1, i2 in zip(ind1, ind2): # mass m1 = self.M[i1] m2 = self.M[i2] # location vector r1 = self.state[i1, :2] r2 = self.state[i2, :2] # velocity vector v1 = self.state[i1, 2:] v2 = self.state[i2, 2:] # relative location & velocity vectors r_rel = r1 - r2 v_rel = v1 - v2 # momentum vector of the center of mass v_cm = (m1 * v1 + m2 * v2) / (m1 + m2) # collisions of spheres reflect v_rel over r_rel rr_rel = np.dot(r_rel, r_rel) vr_rel = np.dot(v_rel, r_rel) v_rel = 2 * r_rel * vr_rel / rr_rel - v_rel # assign new velocities self.state[i1, 2:] = v_cm + v_rel * m2 / (m1 + m2) self.state[i2, 2:] = v_cm - v_rel * m1 / (m1 + m2) # check for crossing boundary crossed_x1 = (self.state[:, 0] < self.bounds[0] + self.size) crossed_x2 = (self.state[:, 0] > self.bounds[1] - self.size) crossed_y1 = (self.state[:, 1] < self.bounds[2] + self.size) crossed_y2 = (self.state[:, 1] > self.bounds[3] - self.size) self.state[crossed_x1, 0] = self.bounds[0] + self.size self.state[crossed_x2, 0] = self.bounds[1] - self.size self.state[crossed_y1, 1] = self.bounds[2] + self.size self.state[crossed_y2, 1] = self.bounds[3] - self.size self.state[crossed_x1 | crossed_x2, 2] *= -1 self.state[crossed_y1 | crossed_y2, 3] *= -1 # add gravity self.state[:, 3] -= self.M * self.G * dt #------------------------------------------------------------ # set up initial state np.random.seed(0) init_state = -0.5 + np.random.random((50, 4)) init_state[:, :2] *= 3.9 box = ParticleBox(init_state, size=0.04) dt = 1. / 30 # 30fps #------------------------------------------------------------ # set up figure and animation fig = plt.figure() fig.subplots_adjust(left=0, right=1, bottom=0, top=1) ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-3.2, 3.2), ylim=(-2.4, 2.4)) # particles holds the locations of the particles particles, = ax.plot([], [], 'bo', ms=6) # rect is the box edge rect = plt.Rectangle(box.bounds[::2], box.bounds[1] - box.bounds[0], box.bounds[3] - box.bounds[2], ec='none', lw=2, fc='none') ax.add_patch(rect) def init(): """initialize animation""" global box, rect particles.set_data([], []) rect.set_edgecolor('none') return particles, rect def animate(i): """perform animation step""" global box, rect, dt, ax, fig box.step(dt) ms = int(fig.dpi * 2 * box.size * fig.get_figwidth() / np.diff(ax.get_xbound())[0]) # update pieces of the animation rect.set_edgecolor('k') particles.set_data(box.state[:, 0], box.state[:, 1]) particles.set_markersize(ms) return particles, rect ani = animation.FuncAnimation(fig, animate, frames=600, interval=10, blit=True, init_func=init) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html #ani.save('particle_box.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show()6 download """ Animation of Elastic collisions with Gravity author: Jake Vanderplas email: website: http://jakevdp.github.com license: BSD Please feel free to use and modify this, but keep the above information. Thanks! """ import numpy as np from scipy.spatial.distance import pdist, squareform import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib.animation as animation class ParticleBox: """Orbits class init_state is an [N x 4] array, where N is the number of particles: [[x1, y1, vx1, vy1], [x2, y2, vx2, vy2], ... ] bounds is the size of the box: [xmin, xmax, ymin, ymax] """ def __init__(self, init_state = [[1, 0, 0, -1], [-0.5, 0.5, 0.5, 0.5], [-0.5, -0.5, -0.5, 0.5]], bounds = [-2, 2, -2, 2], size = 0.04, M = 0.05, G = 9.8): self.init_state = np.asarray(init_state, dtype=float) self.M = M * np.ones(self.init_state.shape[0]) self.size = size self.state = self.init_state.copy() self.time_elapsed = 0 self.bounds = bounds self.G = G def step(self, dt): """step once by dt seconds""" self.time_elapsed += dt # update positions self.state[:, :2] += dt * self.state[:, 2:] # find pairs of particles undergoing a collision D = squareform(pdist(self.state[:, :2])) ind1, ind2 = np.where(D < 2 * self.size) unique = (ind1 < ind2) ind1 = ind1[unique] ind2 = ind2[unique] # update velocities of colliding pairs for i1, i2 in zip(ind1, ind2): # mass m1 = self.M[i1] m2 = self.M[i2] # location vector r1 = self.state[i1, :2] r2 = self.state[i2, :2] # velocity vector v1 = self.state[i1, 2:] v2 = self.state[i2, 2:] # relative location & velocity vectors r_rel = r1 - r2 v_rel = v1 - v2 # momentum vector of the center of mass v_cm = (m1 * v1 + m2 * v2) / (m1 + m2) # collisions of spheres reflect v_rel over r_rel rr_rel = np.dot(r_rel, r_rel) vr_rel = np.dot(v_rel, r_rel) v_rel = 2 * r_rel * vr_rel / rr_rel - v_rel # assign new velocities self.state[i1, 2:] = v_cm + v_rel * m2 / (m1 + m2) self.state[i2, 2:] = v_cm - v_rel * m1 / (m1 + m2) # check for crossing boundary crossed_x1 = (self.state[:, 0] < self.bounds[0] + self.size) crossed_x2 = (self.state[:, 0] > self.bounds[1] - self.size) crossed_y1 = (self.state[:, 1] < self.bounds[2] + self.size) crossed_y2 = (self.state[:, 1] > self.bounds[3] - self.size) self.state[crossed_x1, 0] = self.bounds[0] + self.size self.state[crossed_x2, 0] = self.bounds[1] - self.size self.state[crossed_y1, 1] = self.bounds[2] + self.size self.state[crossed_y2, 1] = self.bounds[3] - self.size self.state[crossed_x1 | crossed_x2, 2] *= -1 self.state[crossed_y1 | crossed_y2, 3] *= -1 # add gravity self.state[:, 3] -= self.M * self.G * dt #------------------------------------------------------------ # set up initial state np.random.seed(0) init_state = -0.5 + np.random.random((50, 4)) init_state[:, :2] *= 3.9 box = ParticleBox(init_state, size=0.04) dt = 1. / 30 # 30fps #------------------------------------------------------------ # set up figure and animation fig = plt.figure() fig.subplots_adjust(left=0, right=1, bottom=0, top=1) ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-3.2, 3.2), ylim=(-2.4, 2.4)) # particles holds the locations of the particles particles, = ax.plot([], [], 'bo', ms=6) # rect is the box edge rect = plt.Rectangle(box.bounds[::2], box.bounds[1] - box.bounds[0], box.bounds[3] - box.bounds[2], ec='none', lw=2, fc='none') ax.add_patch(rect) def init(): """initialize animation""" global box, rect particles.set_data([], []) rect.set_edgecolor('none') return particles, rect def animate(i): """perform animation step""" global box, rect, dt, ax, fig box.step(dt) ms = int(fig.dpi * 2 * box.size * fig.get_figwidth() / np.diff(ax.get_xbound())[0]) # update pieces of the animation rect.set_edgecolor('k') particles.set_data(box.state[:, 0], box.state[:, 1]) particles.set_markersize(ms) return particles, rect ani = animation.FuncAnimation(fig, animate, frames=600, interval=10, blit=True, init_func=init) # save the animation as an mp4. This requires ffmpeg or mencoder to be # installed. The extra_args ensure that the x264 codec is used, so that # the video can be embedded in html5. You may need to adjust this for # your system: for more information, see # http://matplotlib.sourceforge.net/api/animation_api.html #ani.save('particle_box.mp4', fps=30, extra_args=['-vcodec', 'libx264']) plt.show() Toán học nên quen thuộc với bất cứ ai có nền tảng vật lý, và kết quả là khá mê hoặc. Tôi đã mã hóa nó trong một chuyến bay, và cuối cùng chỉ ngồi và xem nó trong khoảng mười phút. Đây chỉ là khởi đầu: nó có thể là một bài tập thú vị để thêm các yếu tố khác, như tính toán nhiệt độ và áp suất để chứng minh luật khí lý tưởng, hoặc vẽ sơ đồ thời gian thực của phân phối vận tốc để theo dõi nó tiếp cận phân phối Maxwellian dự kiến. Nó mở ra nhiều khả năng cho các bản demo vật lý ảo ... Tóm tắt nó lênMô -đun hoạt hình Matplotlib là một bổ sung tuyệt vời cho những gì đã là một gói tuyệt vời.Tôi nghĩ rằng tôi đã gãi bề mặt của những gì có thể với các công cụ này ... bạn có thể nghĩ ra ý tưởng hoạt hình tuyệt vời nào? Chỉnh sửa: Trong một bài đăng tiếp theo, tôi chỉ ra cách các công cụ này có thể được sử dụng để tạo hoạt hình của một hệ thống cơ học lượng tử đơn giản. |