Making Particles Move Randomly Within A Circle [cl

I am attempting to make plotted points move around within a circle of a known radius and center. At the moment I am able to generate points within the boundary, but now I need to make them move.

I have the following script to generate the initial locations of the particles.

function [particle_gen] = generate(n,centre,radius)
    %generates n particles in a circle with specified centre and radius

    angle = 2 * pi * rand(n,1);
    r = radius * sqrt(rand(n,1));
    X = r.*cos(angle) + centre(1);
    Y = r.*sin(angle) + centre(2);
    plot(X,Y,'.k')
end

I want to animate them so that the particles travel in a straight line at a constant velocity until they hit the circular boundary and bounce off. I need this to all happen within the same plot.

回答1:

Ok so there are a few facets to this question that I'm going to address separately.

Updating Plot Objects

As far as the logistics of creating a plot and updating the values of the plot, you can do all of this using MATLAB's handle graphics. If you create a plot object, you can use the returned handle to the plot to update the plot.

plt = plot(x, y, 'k.');
set(plt, 'XData', newx, 'YData', newy)

In your case, the newx and newy values will be computed in a loop and then you can assign them. By doing this, you won't constantly be creating new graphics objects which can get very expensive with MATLAB.

Since you're doing a lot of computation, you'll also want to use drawnow each time through your loop to force MATLAB to render the graphics.

Particle Trajectories

You will need to compute the xnew and ynew values that you are going to update your plot with. There are a number of ways that you could determine these trajectories, but I have presented one way to do that here.

Apply Displacements. Using the velocity and direction, you can apply the displacements to each of the particles.
Determine if point is outside of the circle. If the displacement places the point outside of the circle (basic equation here), then we need to get it back into the circle. If there are no points outside of the circle, then go back to step 1.
For each point outside of the circle, find the intersection of its path with the circle. This is a pretty trivial exercise using this solution as a reference.
Find the tangent line of the circle at this point. Again this is trivial as we can draw a line from the center of the circle to the intersection point on the circle. the negative inverse of the slope of this line is the slope of the tangent line to the curve at the intersection point. By combining this with the intersection point, we have the equation of a line.
Reflect the point that is outside of the circle across this tangent line to "bounce" it back into the circle. This reflection of a point across a line can be performed using the equations provided here.

Example Solution

I have implemented the above steps. Although this may not be the most computationally-efficient method it provides the desired result.

function bounce(npts, vmin, vmax, radius, center)

    % Initial direction/velocity of the points
    direction = rand(npts, 1) * 2 *pi;
    velocity = (rand(npts, 1) * (vmax - vmin)) + vmin;

    % Create random starting locations within the circle
    theta = rand(npts, 1) * 2*pi;
    r = radius * sqrt(rand(npts, 1));

    XY = [r .* cos(theta(:)) + center(1), ...
          r .* sin(theta(:)) + center(2)];

    % Initial plot objects
    hfig = figure('Color', 'w');
    hax = axes('Parent', hfig);

    % Plot the dots as black markers
    hdots = plot(XY(:,1), XY(:,2), ...
                'Parent', hax, ...
                'Marker', '.', ...
                'Color', 'k', ...
                'LineStyle', 'none', ...
                'MarkerSize', 12);

    hold(hax, 'on')
    axis(hax, 'equal')

    % Plot the circle as a reference
    t = linspace(0, 2*pi, 100);
    plot(radius * cos(t) + center(1), ...
         radius * sin(t) + center(2))

    % Keep simulating until we actually close the window
    while ishghandle(hfig);
        % Determine new dot locations
        [XY, direction] = step(XY, direction, velocity, radius, center);

        % Update the dot plot to reflect new locations
        set(hdots, 'XData', XY(:,1), 'YData', XY(:,2))

        % Force a redraw
        drawnow
    end
end

function [XYnew, direction] = step(XY, direction, velocity, radius, center)
    % Compute the next position of the points
    DX = [cos(direction(:)) .* velocity, ...
          sin(direction(:)) .* velocity];
    XYnew = XY + DX;

    % Now check that they are all inside circle
    isOutside = sum(bsxfun(@minus, XYnew, center).^2, 2) > radius^2;

    % The ones that are outside should "bounce" back into the circle
    if any(isOutside)        
        orig  = XY(isOutside,:);
        new   = XYnew(isOutside,:);
        delta = -DX(isOutside,:);

        % Find intersection of this path with the circle
        % Taken from: https://math.stackexchange.com/a/311956
        a = sum(delta.^2, 2);
        b = sum(2 .* delta .* bsxfun(@minus, orig, center), 2);
        c = sum(bsxfun(@minus, orig, center).^2, 2) - radius^2;

        t = (2 * c) ./ (-b + sqrt(b.^2 - 4 .* a .* c)); 
        xintersect = orig(:,1) + delta(:,1) .* t;
        yintersect = orig(:,2) + delta(:,2) .* t;

        % Get tangent at this intersection (slope/intercept form)
        m = - 1 ./ ((yintersect - center(2)) ./ (xintersect - center(1)));
        b = yintersect - m .* xintersect;

        % "Reflect" outside points across the tangent line to "bounce" them
        % Equations from: https://stackoverflow.com/a/3307181/670206
        d = (new(:,1) + (new(:,2) - b) .* m) ./ (1 + m.^2);

        XYnew(isOutside,1) = 2 * d - new(:,1);
        XYnew(isOutside,2) = 2 .* d .* m - new(:,2) + 2 .* b;

        % Recompute the direction of the particles that "bounced"
        direction(isOutside) = atan2(XYnew(isOutside,2) - yintersect, ...
                                     XYnew(isOutside,1) - xintersect);
    end
end

The Result

By running the following command I was able to obtain the following result.

bounce(100, 0.01, 0.2, 5, [0 0]);

回答2:

This is not a trivial problem. Your current method gives you initial starting positions. You will also want to generate starting 'velocities' which I will call VX and VY.

That is, store X,Y position arrays as well as VX and VY direction arrays where each 'line' of VX and VY has a constant norm. You will then increment each position X,Y by their corresponding VX and VY each iteration, then check if you are against a wall, and then change the VX and VY for points hitting the wall in such a way that they 'bounce' off of it.

Hopefully, this gets you started! Good luck.