Projectielbaan

[OOOVincentOOO]

Een laatste:
[ax.set_xlim] en [ax.set_ylim] weglaten uit [def init()] werkt bij mij beter, Slechts eenmaal lopen laten bij programma start:

Zonder handmatig schermafmeting veranderen werkt het op deze manier. Nog enkel een hoogte text label wat niet ververst.

Code: Selecteer alles

# Set up a new Figure, with equal aspect ratio so the ball appears round.
fig, ax = plt.subplots()
ax.set_xlim(0, XMAX)
ax.set_ylim(0, YMAX)
ax.set_aspect('equal')

Excuses voor mijn bemoeienissen. Nu laat ik het erbij.

[edit:]
echter super!!! met sympy en animation!!!

Tekts ook oke:

Code: Selecteer alles

height_text = ax.text(XMAX*0.5, y0*0.8, f'Height: {y0:.1f} m')
height_text.set_text("")

[wnvl1]

Dit is ook nog interessant. Zo kan je meerdere ballen tegelijk laten botsen op de ondergrond van ukster en ze tegelijk ook nog tegen elkaar laten botsen. Wel best vaste wanden introduceren op 0 en 10 dan anders vliegen ze direct uit het spel.

https://scipython.com/blog/two-dimensional-collisions/

[wnvl1]

Zelfde bodem maar veel ballen tegelijk.

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib import animation
from itertools import combinations
import sympy as sym


class Particle:
    """A class representing a two-dimensional particle."""

    def __init__(self, x, y, vx, vy, radius=0.01, color='blue'):
        """Initialize the particle's position, velocity, and radius.

        Any key-value pairs passed in the styles dictionary will be passed
        as arguments to Matplotlib's Circle patch constructor.

        """

        self.r = np.array((x, y))
        self.v = np.array((vx, vy))
        self.radius = radius
        self.color = color


        x = sym.Symbol('x')
        y_expr = sym.sin(x) + (0.2 * sym.cos(4 * x + sym.sin(4 * x))) - 0.2 * x + 3
        y_expr_diff = sym.diff(y_expr, x)

        self.f = sym.lambdify(x, y_expr, "numpy")
        self.f_diff = sym.lambdify(x, y_expr_diff, "numpy")

    # For convenience, map the components of the particle's position and
    # velocity vector onto the attributes x, y, vx and vy.
    @property
    def x(self):
        return self.r[0]
    @x.setter
    def x(self, value):
        self.r[0] = value
    @property
    def y(self):
        return self.r[1]
    @y.setter
    def y(self, value):
        self.r[1] = value
    @property
    def vx(self):
        return self.v[0]
    @vx.setter
    def vx(self, value):
        self.v[0] = value
    @property
    def vy(self):
        return self.v[1]
    @vy.setter
    def vy(self, value):
        self.v[1] = value

    def overlaps(self, other):
        """Does the circle of this Particle overlap that of other?"""

        return np.hypot(*(self.r - other.r)) < self.radius + other.radius

    def draw(self, ax):
        """Add this Particle's Circle patch to the Matplotlib Axes ax."""

        circle = Circle(xy=self.r, radius=self.radius, color=self.color)
        ax.add_patch(circle)
        return circle

    def advance(self, dt):
        """Advance the Particle's position forward in time by dt."""

        self.r += self.v * dt
        self.vy -= 9.81 * dt

        # Make the Particles bounce off the walls
        if self.x - self.radius < 0:
            self.x = self.radius
            self.vx = -self.vx
        if self.x + self.radius > 10:
            self.x = 10-self.radius
            self.vx = -self.vx
        if (self.y - self.radius) < self.f(self.x):
            # bounce!
            self.y = self.radius + self.f(self.x)
            vector_v = np.array([self.vx, self.vy])
            vector_normal = np.array([-self.f_diff(self.x), 1])
            vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
            v_reflected = -(1 + 1) * np.dot(vector_normal, vector_v) * vector_normal + vector_v
            self.vx = v_reflected[0]
            self.vy = v_reflected[1]

        if self.y + self.radius > 10:
            self.y = 10-self.radius
            self.vy = -self.vy

class Simulation:
    """A class for a simple hard-circle molecular dynamics simulation.

    The simulation is carried out on a square domain: 0 <= x < 1, 0 <= y < 1.

    """

    def __init__(self, n, radius=0.01, styles=None):
        """Initialize the simulation with n Particles with radii radius.

        radius can be a single value or a sequence with n values.

        Any key-value pairs passed in the styles dictionary will be passed
        as arguments to Matplotlib's Circle patch constructor when drawing
        the Particles.

        """

        x = sym.Symbol('x')
        y_expr = sym.sin(x) + (0.2 * sym.cos(4 * x + sym.sin(4 * x))) - 0.2 * x + 3
        y_expr_diff = sym.diff(y_expr, x)

        self.f = sym.lambdify(x, y_expr, "numpy")
        self.f_diff = sym.lambdify(x, y_expr_diff, "numpy")

        self.init_particles(n, radius, styles)

    def init_particles(self, n, radius, styles=None):
        """Initialize the n Particles of the simulation.

        Positions and velocities are chosen randomly; radius can be a single
        value or a sequence with n values.

        """

        try:
            iterator = iter(radius)
            assert n == len(radius)
        except TypeError:
            # r isn't iterable: turn it into a generator that returns the
            # same value n times.
            def r_gen(n, radius):
                for i in range(n):
                    yield radius
            radius = r_gen(n, radius)

        self.n = n
        self.particles = []
        for i, rad in enumerate(radius):
            # Try to find a random initial position for this particle.
            while True:
                # Choose x, y so that the Particle is entirely inside the
                # domain of the simulation.
                x = float(i+1)
                y = float(8)

                # Choose a random velocity (within some reasonable range of
                # values) for the Particle.
                vx = float(0)
                vy = float(-1)
                colors = ['black', 'blue', 'yellow', 'orange', 'green', 'purple', 'maroon', 'teal', 'brown']
                particle = Particle(x, y, vx, vy, rad, colors[i])
                # Check that the Particle doesn't overlap one that's already
                # been placed.
                for p2 in self.particles:
                    if p2.overlaps(particle):
                        break
                else:
                    self.particles.append(particle)
                    break

    def handle_collisions(self):
        """Detect and handle any collisions between the Particles.

        When two Particles collide, they do so elastically: their velocities
        change such that both energy and momentum are conserved.

        """

        def change_velocities(p1, p2):
            """
            Particles p1 and p2 have collided elastically: update their
            velocities.

            """

            m1, m2 = p1.radius**2, p2.radius**2
            M = m1 + m2
            r1, r2 = p1.r, p2.r
            d = np.linalg.norm(r1 - r2)**2
            v1, v2 = p1.v, p2.v
            u1 = v1 - 2*m2 / M * np.dot(v1-v2, r1-r2) / d * (r1 - r2)
            u2 = v2 - 2*m1 / M * np.dot(v2-v1, r2-r1) / d * (r2 - r1)
            p1.v = u1
            p2.v = u2

        # We're going to need a sequence of all of the pairs of particles when
        # we are detecting collisions. combinations generates pairs of indexes
        # into the self.particles list of Particles on the fly.
        pairs = combinations(range(self.n), 2)
        for i,j in pairs:
            if self.particles[i].overlaps(self.particles[j]):
                change_velocities(self.particles[i], self.particles[j])

    def advance_animation(self, dt):
        """Advance the animation by dt, returning the updated Circles list."""

        for i, p in enumerate(self.particles):
            p.advance(dt)
            self.circles[i].center = p.r
        self.handle_collisions()
        return self.circles

    def advance(self, dt):
        """Advance the animation by dt."""
        for i, p in enumerate(self.particles):
            p.advance(dt)
        self.handle_collisions()

    def init(self):
        """Initialize the Matplotlib animation."""

        self.circles = []
        for particle in self.particles:
            self.circles.append(particle.draw(self.ax))
        return self.circles

    def animate(self, i):
        """The function passed to Matplotlib's FuncAnimation routine."""

        self.advance_animation(0.01)
        return self.circles

    def do_animation(self, save=False):
        """Set up and carry out the animation of the molecular dynamics.

        To save the animation as a MP4 movie, set save=True.
        """

        fig, self.ax = plt.subplots()

        x1 = np.arange(0, 10, 0.01)
        y1 = self.f(x1)

        self.ax.plot(x1, y1, lw=2)

        for s in ['top','bottom','left','right']:
            self.ax.spines[s].set_linewidth(2)
        self.ax.set_aspect('equal', 'box')
        self.ax.set_xlim(0, 10)
        self.ax.set_ylim(0, 10)
        self.ax.xaxis.set_ticks([])
        self.ax.yaxis.set_ticks([])

        anim = animation.FuncAnimation(fig, self.animate, init_func=self.init,
                               frames=50, interval=80, blit=True)
        if save:
            Writer = animation.writers['ffmpeg']
            writer = Writer(fps=100, bitrate=1800)
            anim.save('collision.mp4', writer=writer)
        else:
            plt.show()


if __name__ == '__main__':
    nparticles = 9
    radii = np.random.random(nparticles)*0.12+0.08
    styles = {'edgecolor': 'C0', 'linewidth': 2, 'fill': None}
    sim = Simulation(nparticles, radii, styles)
    sim.do_animation(save=False)

[wnvl1]

Volgende stap zou zijn afstappen van puntmassa's en een vierkantje laten botsen op deze bodem.

[wnvl1]

Op basis van vergelijking 11 gecombineerd met vergelijkingen 7 tot 10 (onderstaande link) moet dat wel gaan.

https://www.myphysicslab.com/engine2D/collision-en.html

Daarnaast moet nog uitgerekend worden waar de plaats is van de impact met de bodem. Want dat is niet noodzakelijk op een van de hoekpunten van het vierkant. dat is ook telkens het snijpunt van een rechte met de bodem berekenen.

[wnvl1]

Vraag is nu of er een library bestaat in python die zulke dingen al implementeert? Op het eerste gezicht vind ik niets.

[wnvl1]

Hier een foutieve poging van mij voor de projectielbaan van een vierkant (zonder zwaartekracht) op dezelfde ondergrond. Het algoritme zal echter verfijnd moeten worden. Ik heb het gevoel dat wanneer twee knooppunten kort na mekaar botsen het vierkant achter de wand geraakt en het niet meer goed komt. De stoten worden berekend met behoud van impuls en energie zoals beschreven in https://www.myphysicslab.com/engine2D/collision-en.html



Het is zoiezo veel complexer dan met puntmassa's. In python vind ik ook geen voorbeeldcode van anderen.

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle
from matplotlib import animation
from itertools import combinations
import sympy as sym


class SingletonMeta(type):
    """
    The Singleton class can be implemented in different ways in Python. Some
    possible methods include: base class, decorator, metaclass. We will use the
    metaclass because it is best suited for this purpose.
    """

    _instances = {}

    def __call__(cls, *args, **kwargs):
        """
        Possible changes to the value of the `__init__` argument do not affect
        the returned instance.
        """
        if cls not in cls._instances:
            instance = super().__call__(*args, **kwargs)
            cls._instances[cls] = instance
        return cls._instances[cls]


class Bottom(metaclass=SingletonMeta):

    def __init__(self):
        x = sym.Symbol('x')
        y_expr = sym.sin(x) + (0.2 * sym.cos(4 * x + sym.sin(4 * x))) - 0.2 * x + 3
        y_expr_diff = sym.diff(y_expr, x)

        self.__f = sym.lambdify(x, y_expr, "numpy")
        self.__f_diff = sym.lambdify(x, y_expr_diff, "numpy")

    def normal3(self, x):
        vector_normal = np.array([-self.f_diff(x), 1, 0])
        vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
        return vector_normal

    def normal(self, x):
        vector_normal = np.array([-self.f_diff(x), 1])
        vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
        return vector_normal

    @property
    def f(self):
        return self.__f

    @property
    def f_diff(self):
        return self.__f_diff


class Particle:
    """A class representing a two-dimensional particle."""

    def __init__(self, m, x, y, vx, vy, side_length=0.01, color='blue'):
        """Initialize the particle's position, velocity, and radius.

        Any key-value pairs passed in the styles dictionary will be passed
        as arguments to Matplotlib's Circle patch constructor.

        """
        self.center_position = np.array([x, y, 0], dtype=float)
        self.center_velocity = np.array([vx, vy, 0], dtype=float)
        self.theta = float(0)
        self.omega = np.array([0, 0, 0], dtype=float)
        self.side_length = side_length
        self.color = color
        self.elasticity = 1
        self.mass = m
        self.I = m * side_length * side_length / 6
        self.bottom = Bottom()

        delta_x = np.abs(np.cos(self.theta + np.pi/4)*(self.side_length/np.sqrt(2)))
        delta_y = np.abs(np.sin(self.theta + np.pi/4)*(self.side_length/np.sqrt(2)))

        self.corner1 = np.array([self.center_position[0]-delta_x , self.center_position[1]-delta_y, 0])
        self.corner2 = np.array([self.center_position[0]+delta_x , self.center_position[1]-delta_y, 0])
        self.corner3 = np.array([self.center_position[0]+delta_x , self.center_position[1]+delta_y, 0])
        self.corner4 = np.array([self.center_position[0]-delta_x , self.center_position[1]+delta_y, 0])

    @property
    def calculate_corners(self):
        delta_x = np.cos(self.theta + np.pi/4)*(self.side_length/np.sqrt(2))
        delta_y = np.sin(self.theta + np.pi/4)*(self.side_length/np.sqrt(2))

        self.corner1 = np.array([self.center_position[0]-delta_x, self.center_position[1]-delta_y, 0])
        self.corner2 = np.array([self.center_position[0]+delta_y , self.center_position[1]-delta_x, 0])
        self.corner3 = np.array([self.center_position[0]+delta_x , self.center_position[1]+delta_y, 0])
        self.corner4 = np.array([self.center_position[0]-delta_y , self.center_position[1]+delta_x, 0])

    def draw(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""

        square = Rectangle((self.corner1[0], self.corner1[1]), self.side_length, self.side_length, angle=self.theta*180/np.pi, color=self.color)
        return square

    def draw_center(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""
        circle_center = Circle((self.center_position[0], self.center_position[1]), 0.1)
        return circle_center

    def draw_corners(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""
        circle_corner1 = Circle((self.corner1[0], self.corner1[1]), 0.1, color='green')
        circle_corner2 = Circle((self.corner2[0], self.corner2[1]), 0.1, color='yellow')
        circle_corner3 = Circle((self.corner3[0], self.corner3[1]), 0.1, color='red')
        circle_corner4 = Circle((self.corner4[0], self.corner4[1]), 0.1, color='black')
        return circle_corner1, circle_corner2, circle_corner3, circle_corner4


    def intersectionBottom(self):
        division = 10
        for i in range(division):
            x = self.corner1[0] + i*(self.corner2[0]-self.corner1[0])/division

            y = self.corner1[1] + i*(self.corner2[1]-self.corner1[1])/division
            if y < self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner2[0] + i*(self.corner3[0]-self.corner2[0])/division
            y = self.corner2[1] + i*(self.corner3[1]-self.corner2[1])/division
            if y <  self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner3[0] + i * (self.corner4[0] - self.corner3[0]) / division
            y = self.corner3[1] + i * (self.corner4[1] - self.corner3[1]) / division
            if y <  self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner4[0] + i * (self.corner1[0] - self.corner4[0]) / division
            y = self.corner4[1] + i * (self.corner1[1] - self.corner4[1]) / division
            if y <  self.bottom.f(x):
                return True,  np.array([x, y, 0])
        return False, np.array([0, 0, 0])

    def intersectionTop(self):
        if self.corner1[1] > 10:
            return True, self.corner1
        if self.corner2[1] > 10:
            return True, self.corner2
        if self.corner3[1] > 10:
            return True, self.corner3
        if self.corner4[1] > 10:
            return True, self.corner4
        return False, (0, 0)

    def intersectionLeft(self):
        minimum = np.min(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if minimum < 0:
            if self.corner1[0] == minimum:
                return True, self.corner1
            if self.corner2[0] == minimum:
                return True, self.corner2
            if self.corner3[0] == minimum:
                return True, self.corner3
            if self.corner4[0] == minimum:
                return True, self.corner4
        return False, (0, 0)

    def intersectionRight(self):
        maximum = np.max(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if maximum > 10:
            corners = []
            if self.corner1[0] == maximum:
                corners.append(self.corner1)
            if self.corner2[0] == maximum:
                corners.append(self.corner2)
            if self.corner3[0] == maximum:
                corners.append(self.corner3)
            if self.corner4[0] == maximum:
                corners.append(self.corner4)
            x=float(0)
            y=float(0)
            for corner in corners:
                x += corner[0]
                y += corner[1]
            return True, (x/len(corners), y/len(corners), 0)
        return False, (0, 0, 0)

    def velocity3(self, position):
        return self.center_velocity + np.cross(self.omega, position)

    def calculateJ(self, velocity3, normal3, position3 ):
        crossproduct_position_normal = np.cross(position3, normal3)
        j = -(1+self.elasticity) * np.dot(velocity3, normal3) / (1/self.mass + np.dot(crossproduct_position_normal, crossproduct_position_normal)/self.I)
        print("j=", j)
        return j

    def calculate_new_velocity(self, velocity3, normal3, j ):
        # print("auto")
        # print(velocity3)
        # print(normal3)
        # print(j/self.mass)
        return velocity3 + np.multiply(normal3, j/self.mass)

    def calculate_new_angular_velocity3(self, omega3, position3, normal3, j ):
        return omega3 + (np.cross(position3, normal3)*j/self.I)

    def advance(self, dt):
        """Advance the Particle's position forward in time by dt."""

        self.center_position += self.center_velocity * dt
        print(self.center_velocity * dt)
        print(self.center_position)
        print('')

        self.theta += self.omega[2] * dt

        self.calculate_corners

        intersection, coordinates_intersection = self.intersectionBottom()
        if intersection :
            print("botsing bottom")
            normal3 = self.bottom.normal3(coordinates_intersection[0])
            position3 = coordinates_intersection - self.center_position
            velocity3 = self.velocity3(position3)

            print("Normal: ", normal3);
            print("Intersection", coordinates_intersection);
            print("Center position", self.center_position)
            print("Center velocity", self.center_velocity)
            print("Omega", self.omega)

            j = self.calculateJ(velocity3, normal3, position3 )
            new_velocity3 = self.calculate_new_velocity(velocity3, normal3, j )
            self.omega = self.calculate_new_angular_velocity3(self.omega, position3, normal3, j)
            self.center_velocity = new_velocity3 - np.cross(self.omega, position3)

            print("New center velocity", self.center_velocity)
            print("New omega", self.omega)

            print();

        intersection, coordinates_intersection = self.intersectionTop()
        if intersection :
            print("botsing top")

            normal3 = np.array([0, -1, 0],  dtype=float)
            position3 = coordinates_intersection - self.center_position
            velocity3 = self.velocity3(position3)

            j = self.calculateJ(velocity3, normal3, position3 )
            new_velocity3 = self.calculate_new_velocity(velocity3, normal3, j )
            self.omega = self.calculate_new_angular_velocity3(self.omega, position3, normal3, j)
            self.center_velocity = new_velocity3 - np.cross(self.omega, position3)

        intersection, coordinates_intersection = self.intersectionLeft()
        if intersection :
            print("botsing left")
            normal3 = np.array([1, 0, 0])
            position3 = coordinates_intersection - self.center_position
            velocity3 = self.velocity3(position3)
            j = self.calculateJ(velocity3, normal3, position3)
            new_velocity3 = self.calculate_new_velocity(velocity3, normal3, j)
            self.omega = self.calculate_new_angular_velocity3(self.omega, position3, normal3, j)
            self.center_velocity = new_velocity3 - np.cross(self.omega, position3)

        intersection, coordinates_intersection = self.intersectionRight()
        if intersection :
            print("botsing right")
            normal3 = np.array([-1, 0, 0], )
            position3 = coordinates_intersection - self.center_position
            velocity3 = self.velocity3(position3)

            print("position3", position3)
            print("velocity3", velocity3)

            print("Normal: ", normal3);
            print("Intersection", coordinates_intersection);
            print("Center position", self.center_position)
            print("Center velocity", self.center_velocity)
            print("Omega", self.omega)

            j = self.calculateJ(velocity3, normal3, position3)
            new_velocity3 = self.calculate_new_velocity(velocity3, normal3, j)
            self.omega = self.calculate_new_angular_velocity3(self.omega, position3, normal3, j)
            self.center_velocity = new_velocity3 - np.cross(self.omega, position3)

            print("New center velocity", self.center_velocity)
            print("New omega", self.omega)

            print()



class Simulation:
    def __init__(self, side_length=0.01):
        self.bottom = Bottom()
        self.particle = Particle(1, 2, 8, 0, -1, side_length, 'blue')

    def init(self):
        """Initialize the Matplotlib animation."""
        objects_to_draw = []
        objects_to_draw.append(self.particle.draw())
        ax.add_patch(objects_to_draw[0])

        self.velocity_x_text = ax.text(6, 9, '')
        self.velocity_y_text = ax.text(6, 8, '')
        self.theta_text = ax.text(6, 7, '')
        self.omega_text = ax.text(6, 6, '')

        return objects_to_draw

    def advance_animation(self, dt):
        self.particle.advance(dt)
        return

    def animate(self, i):
        global fig, ax
        self.advance_animation(0.01)
        objects_to_draw = []
        objects_to_draw.append(self.particle.draw())
        ax.add_patch(objects_to_draw[0])

        objects_to_draw.append(self.particle.draw_center())
        ax.add_patch(objects_to_draw[1])

        corner1, corner2, corner3, corner4 = self.particle.draw_corners()
        objects_to_draw.append(corner1)
        objects_to_draw.append(corner2)
        objects_to_draw.append(corner3)
        objects_to_draw.append(corner4)
        ax.add_patch(objects_to_draw[2])
        ax.add_patch(objects_to_draw[3])
        ax.add_patch(objects_to_draw[4])
        ax.add_patch(objects_to_draw[5])

        self.velocity_x_text.set_text(f'Vx: {self.particle.center_velocity[0]:.1f} m')
        self.velocity_y_text.set_text(f'Vy: {self.particle.center_velocity[1]:.1f} m')
        theta_degrees = self.particle.theta * 180 / np.pi
        self.theta_text.set_text(f'Theta: {theta_degrees:.1f} graden')
        self.omega_text.set_text(f'Omega: {self.particle.omega[2]:.1f} rad/s')

        objects_to_draw.append(self.velocity_x_text)
        objects_to_draw.append(self.velocity_y_text)
        objects_to_draw.append(self.theta_text)
        objects_to_draw.append(self.omega_text)

        return objects_to_draw

    def do_animation(self):
        global fig, ax
        anim = animation.FuncAnimation(fig, self.animate, init_func=self.init,
                               frames=50, interval=1, blit=True)
        plt.show()

if __name__ == '__main__':
    global fig, ax
    fig, ax = plt.subplots()

    x1 = np.arange(0, 10, 0.01)
    bottom = Bottom()
    y1 = bottom.f(x1)
    ax.plot(x1, y1, lw=2)

    for s in ['top', 'bottom', 'left', 'right']:
        ax.spines[s].set_linewidth(2)

    ax.set_aspect('equal', 'box')
    ax.set_xlim(0, 10)
    ax.set_ylim(0, 10)
    ax.xaxis.set_ticks([])
    ax.yaxis.set_ticks([])

    side_length = 1
    sim = Simulation(side_length)
    sim.do_animation()

[wnvl1]

Nu een veel betere simulatie van een vierkant botsend in een zwaartekracht veld nog altijd met dezelfde ondergrond. Om te debuggen druk ik ook altijd op het scherm de kinetische energie (translatie + rotatie), de potentiele energie en de totale energie af. Opvallend is dat de totale energie langzaam toeneemt. Als ik het zwaartekracht veld gelijk stel aan nul verdwijnt het probleem. Het is dus niet gerelateerd aan de botsingen, want die zijn er ook zonder zwaartekrachtveld. Het probleem lijkt mij mogelijk een numeriek probleem te zijn gelinkt aan de voorwaartse Euler methode die gebruikt wordt.

Code: Selecteer alles

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle
from matplotlib import animation
from itertools import combinations
import sympy as sym


class SingletonMeta(type):
    """
    The Singleton class can be implemented in different ways in Python. Some
    possible methods include: base class, decorator, metaclass. We will use the
    metaclass because it is best suited for this purpose.
    """

    _instances = {}

    def __call__(cls, *args, **kwargs):
        """
        Possible changes to the value of the `__init__` argument do not affect
        the returned instance.
        """
        if cls not in cls._instances:
            instance = super().__call__(*args, **kwargs)
            cls._instances[cls] = instance
        return cls._instances[cls]


class Bottom(metaclass=SingletonMeta):

    def __init__(self):
        x = sym.Symbol('x')
        y_expr = sym.sin(x) + (0.2 * sym.cos(4 * x + sym.sin(4 * x))) - 0.2 * x + 3
        #y_expr = 0.1 * x
        #y_expr = 1+sym.sin(x*np.pi/4)
        y_expr_diff = sym.diff(y_expr, x)

        self.__f = sym.lambdify(x, y_expr, "numpy")
        self.__f_diff = sym.lambdify(x, y_expr_diff, "numpy")

    def normal(self, x):
        vector_normal = np.array([-self.f_diff(x), 1, 0])
        vector_normal = vector_normal / np.sqrt(np.sum(vector_normal ** 2))
        return vector_normal

    @property
    def f(self):
        return self.__f

    @property
    def f_diff(self):
        return self.__f_diff


class Particle:
    """A class representing a two-dimensional particle."""

    def __init__(self, m, x, y, vx, vy, side_length=0.01, color='blue'):
        """Initialize the particle's position, velocity, and radius.

        Any key-value pairs passed in the styles dictionary will be passed
        as arguments to Matplotlib's Circle patch constructor.

        """
        self.center_position = np.array([x, y, 0], dtype=float)
        self.center_velocity = np.array([vx, vy, 0], dtype=float)
        self.theta = float(0)
        self.omega = np.array([0, 0, 0], dtype=float)
        self.side_length = side_length
        self.color = color
        self.elasticity = 1
        self.g = 9.81
        self.mass = m
        self.I = m * (side_length**2) / 6
        self.bottom = Bottom()
        self.energy_kin = float(0)
        self.energy_pot = float(0)
        self.collision_position = np.array([0, 0, 0], dtype=float)
        self.collision = False

        self.calculate_corners

    @property
    def calculate_corners(self):
        delta_x = np.cos(self.theta + np.pi/4)*(self.side_length/np.sqrt(2))
        delta_y = np.sin(self.theta + np.pi/4)*(self.side_length/np.sqrt(2))

        self.corner1 = np.array([self.center_position[0]-delta_x, self.center_position[1]-delta_y, 0])
        self.corner2 = np.array([self.center_position[0]+delta_y , self.center_position[1]-delta_x, 0])
        self.corner3 = np.array([self.center_position[0]+delta_x , self.center_position[1]+delta_y, 0])
        self.corner4 = np.array([self.center_position[0]-delta_y , self.center_position[1]+delta_x, 0])

    def draw(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""

        square = Rectangle((self.corner1[0], self.corner1[1]), self.side_length, self.side_length, angle=self.theta*180/np.pi, color=self.color)
        return square

    def draw_center(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""
        circle_center = Circle((self.center_position[0], self.center_position[1]), 0.1)
        return circle_center

    def draw_corners(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""
        circle_corner1 = Circle((self.corner1[0], self.corner1[1]), 0.1, color='green')
        circle_corner2 = Circle((self.corner2[0], self.corner2[1]), 0.1, color='yellow')
        circle_corner3 = Circle((self.corner3[0], self.corner3[1]), 0.1, color='red')
        circle_corner4 = Circle((self.corner4[0], self.corner4[1]), 0.1, color='black')
        return circle_corner1, circle_corner2, circle_corner3, circle_corner4

    def draw_collision(self):
        """Add this Particle's Rectangle patch to the Matplotlib Axes ax."""
        circle_collision = Circle((self.collision_position[0], self.collision_position[1]), 0.1, color='pink')
        return circle_collision

    def intersectionBottom(self):
        division = 10
        for i in range(division):
            x = self.corner1[0] + i*(self.corner2[0]-self.corner1[0])/division

            y = self.corner1[1] + i*(self.corner2[1]-self.corner1[1])/division
            if y <= self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner2[0] + i*(self.corner3[0]-self.corner2[0])/division
            y = self.corner2[1] + i*(self.corner3[1]-self.corner2[1])/division
            if y <=  self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner3[0] + i * (self.corner4[0] - self.corner3[0]) / division
            y = self.corner3[1] + i * (self.corner4[1] - self.corner3[1]) / division
            if y <=  self.bottom.f(x):
                return True,  np.array([x, y, 0])
            x = self.corner4[0] + i * (self.corner1[0] - self.corner4[0]) / division
            y = self.corner4[1] + i * (self.corner1[1] - self.corner4[1]) / division
            if y <=  self.bottom.f(x):
                return True,  np.array([x, y, 0])
        return False, np.array([0, 0, 0])

    def intersectionTop(self):
        if self.corner1[1] >= 10:
            return True, self.corner1
        if self.corner2[1] >= 10:
            return True, self.corner2
        if self.corner3[1] >= 10:
            return True, self.corner3
        if self.corner4[1] >= 10:
            return True, self.corner4
        return False, (0, 0)

    def intersectionLeft(self):
        minimum = np.min(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if minimum <= 0:
            if self.corner1[0] == minimum:
                return True, self.corner1
            if self.corner2[0] == minimum:
                return True, self.corner2
            if self.corner3[0] == minimum:
                return True, self.corner3
            if self.corner4[0] == minimum:
                return True, self.corner4
        return False, (0, 0)

    def intersectionRight(self):
        maximum = np.max(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if maximum >= 10:
            corners = []
            if self.corner1[0] == maximum:
                corners.append(self.corner1)
            if self.corner2[0] == maximum:
                corners.append(self.corner2)
            if self.corner3[0] == maximum:
                corners.append(self.corner3)
            if self.corner4[0] == maximum:
                corners.append(self.corner4)
            x=float(0)
            y=float(0)
            for corner in corners:
                x += corner[0]
                y += corner[1]
            return True, (x/len(corners), y/len(corners), 0)
        return False, (0, 0, 0)

    def calculate_velocity_contact(self, position_relative_contact):
        return self.center_velocity + np.cross(self.omega, position_relative_contact)

    def calculate_impulse_J(self, position_relative_contact, velocity_contact, normal_contact):
        crossproduct_position_normal = np.cross(position_relative_contact, normal_contact)
        j = -(1+self.elasticity) * np.dot(velocity_contact, normal_contact) / ((1/self.mass) + (np.dot(crossproduct_position_normal, crossproduct_position_normal)/self.I))
        return j

    def calculate_new_velocity_center(self, normal_contact, j ):
        print("vermenigvuldiging")
        print(normal_contact)
        print(j/self.mass)
        print(np.multiply(normal_contact, j / self.mass))

        return self.center_velocity + np.multiply(normal_contact, j/self.mass)

    def calculate_new_angular_velocity(self, position_relative_contact, normal_contact, j ):
        return self.omega + (np.cross(position_relative_contact, normal_contact)*j/self.I)

    def calculate_collision(self, coordinates_intersection, normal_contact):
        position_relative_contact = coordinates_intersection - self.center_position
        velocity_contact = self.calculate_velocity_contact(position_relative_contact)


        print("Massa", self.mass)
        print("Traagheidsmoment", self.I)

        print("VOOR botsing:")
        print("-------------")
        print("Normal contact: ", normal_contact)
        print("Intersection: ", coordinates_intersection)
        print("Center position: ", self.center_position)
        print("Relative position contact w.r.t. center: ", position_relative_contact)
        print("Velocity contact: ", velocity_contact)
        print("Velocity contact normal: ", np.dot(velocity_contact, normal_contact) * normal_contact)
        print("Velocity contact tangential: ", velocity_contact - np.dot(velocity_contact, normal_contact) * normal_contact)
        print("Center velocity: ", self.center_velocity)
        print("Omega: ", self.omega)

        print("Energy translation: ", (self.mass * np.dot(self.center_velocity, self.center_velocity))/2)
        print("Energy rotation: ", (self.I * (self.omega[2]**2))/2)
        print("Energy potential", self.mass * self.g * self.center_position[1])
        print("Energy total", (self.mass * np.dot(self.center_velocity, self.center_velocity))/2 + (self.I * (self.omega[2]**2))/2 + self.mass * self.g * self.center_position[1])


        j = self.calculate_impulse_J(position_relative_contact, velocity_contact, normal_contact)
        self.center_velocity = self.calculate_new_velocity_center(normal_contact, j)
        self.omega = self.calculate_new_angular_velocity(position_relative_contact, normal_contact, j)
        new_velocity_contact = self.center_velocity + np.cross(self.omega, position_relative_contact)

        print();
        print("NA botsing:")
        print("-----------")
        print("New velocity contact", new_velocity_contact)
        print("New velocity contact normal", np.dot(new_velocity_contact, normal_contact) * normal_contact)
        print("New velocity contact tangential", new_velocity_contact - np.dot(new_velocity_contact, normal_contact) * normal_contact)
        print("New velocity center", self.center_velocity)
        print("New omega", self.omega)
        print("Energy translation: ", (self.mass * np.dot(self.center_velocity, self.center_velocity))/2)
        print("Energy rotation: ", (self.I * (self.omega[2]**2))/2)
        print("Energy potential", self.mass * self.g * self.center_position[1])
        print("Energy total", (self.mass * np.dot(self.center_velocity, self.center_velocity))/2 + (self.I * (self.omega[2]**2))/2 + self.mass * self.g * self.center_position[1])
        print();

    def advance(self, dt):
        """Advance the Particle's position forward in time by dt."""
        self.center_position += self.center_velocity * dt
        self.theta += self.omega[2] * dt
        self.center_velocity[1] += -self.g * dt
        self.calculate_corners
        self.collision = False

        self.energy_kin = (self.mass * (np.dot(self.center_velocity, self.center_velocity)) + self.I * (self.omega[2]**2))/2
        self.energy_pot = self.mass * self.g * self.center_position[1]

        maximum = np.max(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if (maximum > 10):
            self.center_position[0] -= maximum-10

        minimum = np.min(np.array([self.corner1[0], self.corner2[0], self.corner3[0], self.corner4[0]]))
        if (minimum < 0):
            self.center_position[0] -= minimum

        maximum = np.max(np.array([self.corner1[1], self.corner2[1], self.corner3[1], self.corner4[1]]))
        if (maximum > 10):
            self.center_position[1] -= maximum-10

        self.calculate_corners

        intersection, coordinates_intersection = self.intersectionBottom()
        if intersection :
            print("botsing bottom")
            self.collision_position = np.array([coordinates_intersection[0], coordinates_intersection[1], 0], dtype=float)
            self.collision = True

            if (self.collision_position[1] < self.bottom.f(self.collision_position[0])):
                self.center_position[1] += self.bottom.f(self.collision_position[0]) - self.collision_position[1]
                self.calculate_corners

            normal_contact = self.bottom.normal(coordinates_intersection[0])
            self.calculate_collision(coordinates_intersection, normal_contact)


        intersection, coordinates_intersection = self.intersectionTop()
        if intersection :
            print("botsing top")
            self.collision_position = np.array([coordinates_intersection[0], coordinates_intersection[1], 0], dtype=float)
            self.collision = True

            normal_contact = np.array([0, -1, 0],  dtype=float)
            self.calculate_collision(coordinates_intersection, normal_contact)

        intersection, coordinates_intersection = self.intersectionLeft()
        if intersection :
            print("botsing left")
            self.collision_position = np.array([coordinates_intersection[0], coordinates_intersection[1], 0], dtype=float)
            self.collision = True
            normal_contact = np.array([1, 0, 0],  dtype=float)
            self.calculate_collision(coordinates_intersection, normal_contact)

        intersection, coordinates_intersection = self.intersectionRight()
        if intersection :
            print("botsing right")
            self.collision_position = np.array([coordinates_intersection[0], coordinates_intersection[1], 0], dtype=float)
            self.collision = True
            normal_contact = np.array([-1, 0, 0])
            self.calculate_collision(coordinates_intersection, normal_contact)



class Simulation:
    def __init__(self, side_length=0.01):
        self.bottom = Bottom()
        self.particle = Particle(1, 2, 8, 0, -1, side_length, 'blue')

    def init(self):
        """Initialize the Matplotlib animation."""
        objects_to_draw = []
        objects_to_draw.append(self.particle.draw())
        ax.add_patch(objects_to_draw[0])

        self.velocity_x_text = ax.text(1, 9, '')
        self.velocity_y_text = ax.text(1, 8, '')
        self.energy_kin_text = ax.text(1, 7, '')
        self.energy_pot_text = ax.text(1, 6, '')
        self.theta_text = ax.text(6, 9, '')
        self.omega_text = ax.text(6, 8, '')
        self.energy_tot_text = ax.text(6, 7, '')
        return objects_to_draw

    def advance_animation(self, dt):
        self.particle.advance(dt)
        return

    def animate(self, i):
        global fig, ax
        self.advance_animation(0.005)
        objects_to_draw = []
        objects_to_draw.append(self.particle.draw())
        ax.add_patch(objects_to_draw[0])

        objects_to_draw.append(self.particle.draw_center())
        ax.add_patch(objects_to_draw[1])

        corner1, corner2, corner3, corner4 = self.particle.draw_corners()
        objects_to_draw.append(corner1)
        objects_to_draw.append(corner2)
        objects_to_draw.append(corner3)
        objects_to_draw.append(corner4)

        ax.add_patch(objects_to_draw[2])
        ax.add_patch(objects_to_draw[3])
        ax.add_patch(objects_to_draw[4])
        ax.add_patch(objects_to_draw[5])

        if True:
            collision = self.particle.draw_collision()
            objects_to_draw.append(collision)
            ax.add_patch(collision)

        self.velocity_x_text.set_text(f'Vx: {self.particle.center_velocity[0]:.1f} m')
        self.velocity_y_text.set_text(f'Vy: {self.particle.center_velocity[1]:.1f} m')

        self.energy_kin_text.set_text(f'E Kin: {self.particle.energy_kin:.1f} J')
        self.energy_pot_text.set_text(f'E Pot: {self.particle.energy_pot:.1f} J')
        self.energy_tot_text.set_text(f'E Tot: {self.particle.energy_pot + self.particle.energy_kin:.1f} J')
        theta_degrees = self.particle.theta * 180 / np.pi
        self.theta_text.set_text(f'Theta: {theta_degrees:.1f} °')
        self.omega_text.set_text(f'Omega: {self.particle.omega[2]:.1f} rad/s')

        objects_to_draw.append(self.velocity_x_text)
        objects_to_draw.append(self.velocity_y_text)
        objects_to_draw.append(self.energy_pot_text)
        objects_to_draw.append(self.energy_kin_text)
        objects_to_draw.append(self.energy_tot_text)
        objects_to_draw.append(self.theta_text)
        objects_to_draw.append(self.omega_text)

        return objects_to_draw

    def do_animation(self):
        global fig, ax
        anim = animation.FuncAnimation(fig, self.animate, init_func=self.init,
                               frames=50, interval=1, blit=True)
        plt.show()

if __name__ == '__main__':
    global fig, ax
    fig = plt.figure()
    ax = fig.add_subplot()

    x1 = np.arange(0, 10, 0.01)
    bottom = Bottom()
    y1 = bottom.f(x1)
    ax.plot(x1, y1, lw=2)

    for s in ['top', 'bottom', 'left', 'right']:
        ax.spines[s].set_linewidth(2)

    ax.set_aspect('equal', 'box')
    ax.set_xlim(0, 10)
    ax.set_ylim(0, 10)
    ax.xaxis.set_ticks([])
    ax.yaxis.set_ticks([])

    side_length = 1
    sim = Simulation(side_length)
    sim.do_animation()

[physicalattraction]

Mooie simulatie, wnvl1!

[ReinW]

Dit commentaar komt als mosterd na de maaltijd, want de discussie is op 16-11-2021 al beëindigd. Misschien toch interessant voor degenen die een “echte” projectiel baan willen berekenen.

Wetenschapsforum

(16.4 KiB) 75 keer gedownload

[wnvl1]

In de eerste vergelijkingen die ik in dit topic heb gezet ontbrak de m in mg, maar in de latere berekeningen hebben ik en de anderen dat wel correct meegenomen.

Jouw berekeningen kloppen alleen maar bij een verticale beweging, lijkt mij. Als er ook een horizontale compenent is, kom je niet tot de verticale krachtenbalans zoals in jouw document staat. Je moet dan echt de differentiaalvergelijkingen uit mijn posts gebruiken. Je gebruikt ook een andere definitie van Cd, lijkt mij. Dat leidt tot een andere interpretatie. Daarom niet verkeerd, maar was niet de bedoeling van ukster, lijkt mij.

[ReinW]

Dat er in je latere berekeningen mg is meegenomen is me ontgaan, sorry. In ieder geval is het handig om door m te delen, dan staat deze alleen in de noemer van de drag force en is de vergelijking wat makkelijker op te lossen.
Mijn berekeningen gelden voor de totale beweging, de vergelijking is de ontbondene in de y-richting. De vergelijking in x-richting is simpel: du/dt = -Fdx. Dat kun je inderdaad combineren tot één differentiaalvergelijking. Oplossing hiervan levert de genoemde getallen.
De drag force coëfficiënt Cd is een dimensieloos getal, welbekend in de stromingsleer voor verschillende vormen van vallende voorwerpen. Ik ken geen andere definities.

[Xilvo]

Je kunt de berekening niet onafhankelijk voor x- en y-richting doen bij een kracht die van de snelheid in het kwadraat afhangt.

Ook kun je niet stellen dat de k met g in de noemer fout is. Het is een andere definitie, en die kwam uit het artikel waar Ukster de vraag uit had. Maar daar is al uitgebreid over gediscussieerd.

Voor een balletje met een straal van ongeveer 4 cm en een dichtheid van ongeveer 900 kg/m³ is die gegeven k een realistische waarde.

Tenslotte, je schrijft

De snelheid neemt af van 50.sin(40) = 32,139 m/s tot 0 m/s
//
De bereikte maximale hoogte is dan 0,5.9,81.t2 = 52,65 m. De Fd is niet nul, dus de bereikte hoogte is altijd lager dan 52,65 m. Waarschijnlijk is het teken van Fd niet juist

Die F_D is niet fout, er wordt daar gerekend met een initiële hoek van 45 graden.

[wnvl1]

De drag force coëfficiënt Cd is een dimensieloos getal, welbekend in de stromingsleer voor verschillende vormen van vallende voorwerpen. Ik ken geen andere definities.

Je hebt wel gelijk, maar deze oefening was waarschijnlijk door ukster geplaatst vanuit het standpunt om wat te leren (of als amusement) omtrent het oplossen van systemen van differentiaal vergelijkingen of omtrent mechanica in het algemeen. Het zou dan vervelend zijn om de Cd zoals die gekend is in de fluidummechanica te gebruiken. Dan moet je die parameters van lucht allemaal invullen terwijl dat niet de kern is van de oefening. Je kan dat perfect incorporeren in een nieuwe "C" zoals we gedaan hebben.

[ReinW]

Je kunt de berekening niet onafhankelijk voor x- en y-richting doen bij een kracht die van de snelheid in het kwadraat afhangt.

Dat doe ik ook niet! Sorry voor het misverstand, dat blijkbaar ontstaat doordat ik de vergelijking in x-richting even buiten beschouwing heb gelaten.
Mijn aanpak is gewoon hetzelfde als in het Python programma: v is ontbonden in vx en vy, a in ax en ay, en daarmee worden vx.append, vy.append, enz. uitgerekend. Alleen staat mijn programma in C# en gebruik ik een tijdstap van 0,00001 s in plaats van 0,05 of 0,01 s. Bij een kleinere tijdstap neemt de nauwkeurigheid behoorlijk toe.

De oorspronkelijke de vraag was heel duidelijk. Kan de analytische berekening geverifieerd worden met een numerieke berekening. Het antwoord is ja en is uitgebreid aangetoond. Dan mag je kiezen wat je wilt voor Fdrag, Cd, m of k, als iedereen maar hetzelfde doet. Er komen waarden langs van m = 207 gram, 1 kg, 4,08 kg; de lanceerhoek is van 40 naar 45 naar 80 graden gegaan en de lanceersnelheid is 50 of 450 m/s. Prima, ik heb geen enkel commentaar op de oefening als zodanig!

Wat betreft de drag force waarvan ik veronderstelde dat het teken niet juist is in de vergelijking dvy/dt = -g - k.vy^2. Die vergelijking geldt alleen voor de beweging omhoog; zodra de snelheid omdraait wordt het dvy/dt = -g + k.vy^2. Het is een veelgemaakte fout, vandaar mijn opmerking. Pas in de Python code van Xilvo (30-10, 19:05) is dat netjes opgelost door hoek van de baan t.o.v. de x-as te gebruiken. Als deze negatief wordt, draait het teken van Fd om.

Het enige punt dat ik wilde maken is dat er geen fysische betekenis aan de resultaten valt toe te kennen. Wat voor voorwerp wordt er omhoog geschoten (wat voor oppervlak); beweegt het door water of een ander medium (de dichtheid), enz. Dat is allemaal samengevat in een k. Maar als de termen in de vergelijkingen verschillende dimensies hebben, wordt ik achterdochtig, want dat leidt vaak tot betekenisloze resultaten.

Volgens mijn bescheiden mening zijn we het (voor een groot deel) eens.

PS.

Voor een balletje met een straal van ongeveer 4 cm en een dichtheid van ongeveer 900 kg/m3 is die gegeven k een realistische waarde.

Welke definitie van k en welke waarde gebruik je?