在模拟弹跳粒子时，最大时间略有不同会产生截然不同的中间结果。我究竟做错了什么？

Question

我正在尝试模拟一个粒子在一个盒子里弹跳，盒子的壁以给定的速度定律（现在采用锯齿波）移动，并将粒子的能量记录为时间的函数。虽然我做了一些不完全准确的近似值（我将在我的代码示例中对此进行扩展），但我运行遇到了一个问题，即最大时间略有不同会在我预期的最大时间以下产生截然不同的结果结果一致（参见数字：在 this figure, the maximum time is 200 while in this one 中是 201。这些数字来自代码的二维版本，但即使在一维中问题仍然存在）。

在一维设置中，我的尝试过程如下；

首先，我设置了最小时间、最大时间和时间步长，然后以时间步长为间隔生成一个从最小时间到最大时间的时间数组：

TMin = 0;
TMax = 200;
TimeStep = .0001;
times = linspace(TMin,TMax,(TMax-TMin)/TimeStep);

据我所知，此步骤正在努力正确生成我想要的时间数组。接下来我生成与时间数组大小相同的粒子位置和速度数组；

x = zeros(1,TNum(2));
vx = zeros(1,TNum(2));

在这一步我还设置了粒子的初始条件。接下来我为墙壁的振荡设置参数，包括振幅和频率，以定义速度定律，例如

VLeft = ALeft*FLeft*sawtooth(FLeft.*times+SLeft)

其中 VLeft 是左壁的速度，ALeft 是位置振荡的幅度，FLeft 是振荡的频率，SLeft 是相对于另一壁的相移。据我所知，这一步也是在做我想做的。接下来，我为墙设置初始条件，并使用黎曼和近似根据速度定律生成墙的位置；

Left(1)=-1; %For example...
for i=2:TNum(2)
    Left(i)=Left(i-1)+VLeft(i).*(TimeStep);
end

只要我为 TimeStep 使用足够小的值（根据我的计算方式，这就是我所期望的），这似乎会生成作为时间函数的墙位置的正确形式).右墙也进行了类似的时间演化。

下一步是我尝试生成粒子的位置和速度作为时间的函数，我认为这一定是出了问题的地方。

for i=2:TNum(2)

    x(i)=x(i-1)+vx(i-1)*(TimeStep);%distance=velocity*time

    if x(i) < Right(i) && x(i) > Left(i) 
    vx(i)=vx(i-1); %If it doesn't hit a wall, velocity doesn't change
    end

    if x(i) < Left(i) || x(i) == Left(i) %If it hits the left wall, bounce off
    vx(i)=-vx(i-1)+VLeft(i);
    x(i) = Left(i)+abs(x(i)-Left(i));
    end

    if x(i) > Right(i) || x(i) == Right(i) %If it hits the right wall, bounce off
    vx(i)=-vx(i-1)+VRight(i);
    x(i) = Right(i)-abs(Right(i)-x(i));
    end
end

对于碰撞，我假设墙壁是无限大的。我知道的近似值是不正确的，如果粒子的位置超出了墙壁的位置，我可以简单地反转它相对于墙壁的位置并给它适当的速度 - 实际计算它的位置，我将不得不以某种方式计算确切的碰撞时间。这似乎是唯一可能发生问题的地方，但我不明白这将如何导致最长时间的变化以将结果更改为低于最长时间。

最后我计算能量为

Energy = zeros(1,TNum(2));

for i=1:TNum(2)
    Energy(i)=vx(i)^2+vy(i)^2;
end

我真的不明白为什么同一时间步长的最大时间变化会产生如此不同的结果。减少时间步长会增加结果匹配的区域，但如果我理解代码正确执行的操作，那么最大时间的增加不应该改变较低时间发生的任何事情，因为时间数组的值是直到那一刻都一样。

很抱歉这个问题冗长笼统，但我真的不确定自己做错了什么，希望提供尽可能多的信息。

编辑：我在下面添加了二维版本的完整代码。

%--------- THE TIME BLOCK ------------------------------------------------%

TMin = 0;
TMax = 201;
TimeStep = .001;
TDensity = (TMax-TMin)/TimeStep; %Declare instance variables that determine the time steps

times = linspace(TMin,TMax,TDensity); %Generate the array of time steps

 TNum = size(times); %For looping purposes and generating velocity, position arrays


%-------- THE OBJECT INITIAL CONDITIONS BLOCK-----------------------------%

Lx = 1; %Half the length of the box, which spans from -Lx to Lx
Ly = 1; %Half the width of the box, from -Ly to Ly

x0 = 0; %Initial x position
v0x = .1; %Initial x velocity
y0 = 0; %Initial y position
v0y = .1; %Initial y velocity

x = zeros(1,TNum(2));
vx = zeros(1,TNum(2)); %Initialize the arrays for position and velocity
y = zeros(1,TNum(2));
vy = zeros(1,TNum(2));

%------ WALL INITIAL CONDITIONS BLOCK ------------------------------------%

ALeft = 0.15; %Determines the size of the oscillations
ARight = 0.2;
ATop = 0.25;
ABottom = 0.12;

FLeft = 1; %Frequencies
FRight = 0.2;
FTop = 3;
FBottom = 0.2;


SLeft = pi/2; %Phase _S_hifts
SRight = 0;

STop = pi/2;
SBottom = 0;

Left = zeros(1,TNum(2));
Right = zeros(1,TNum(2));
Top = zeros(1,TNum(2));
Bottom = zeros(1,TNum(2));

VLeft = ALeft*FLeft*sawtooth(FLeft.*times+SLeft);
VRight = ARight*FRight*sawtooth(FRight.*times+SRight);
VTop = ATop*FTop*sawtooth(FTop.*times+STop);
VBottom = ABottom*FBottom*sawtooth(FBottom.*times+SBottom);

%---- TIME TO MOVE -------------------------------------------------------%

x(1)=x0; %Set the initial conditions for the variable arrays
y(1)=y0;
vx(1)=v0x;
vy(1)=v0y;
Left(1)=-Lx;
Right(1)=Lx;
Top(1)=Ly;
Bottom(1)=-Ly;

for i=2:TNum(2)
    Left(i)=Left(i-1)+(1/2)*TimeStep*(VLeft(i-1)+VLeft(i)); 
    Right(i)=Right(i-1)+(1/2)*TimeStep*(VRight(i-1)+VRight(i));
    Top(i)=Top(i-1)+(1/2)*TimeStep*(VTop(i-1)+VTop(i));
    Bottom(i)=Bottom(i-1)+(1/2)*TimeStep*(VBottom(i-1)+VBottom(i));
end

for i=2:TNum(2)

    x(i)=x(i-1)+vx(i-1)*(TimeStep);%distance=velocity*time
    y(i)=y(i-1)+vy(i-1)*(TimeStep);

    if x(i) < Right(i) && x(i) > Left(i) %If it doesn't hit a wall, velocity doesn't change
        vx(i)=vx(i-1);
    end

    if y(i) < Top(i) && y(i) > Bottom(i) %If it doesn't hit a wall, velocity doesn't change
        vy(i)=vy(i-1);
    end

    if x(i) < Left(i) %If it hits the left wall, bounce off
        vx(i)=-vx(i-1)+VLeft(i);
        x(i) = Left(i)+abs(x(i)-Left(i));
    end

    if x(i) > Right(i) %If it hits the right wall, bounce off
        vx(i)=-vx(i-1)+VRight(i);
        x(i) = Right(i)-abs(Right(i)-x(i));
    end

    if y(i) < Bottom(i) %If it hits the bottom wall, bounce off
        vy(i)=-vy(i-1)+VBottom(i);
        y(i) = Bottom(i)+abs(y(i)-Bottom(i));
    end

    if y(i) > Top(i) %If it hits the top wall, bounce off
        vy(i)=-vy(i-1)+VTop(i);
        y(i) = Top(i)-abs(Top(i)-y(i));
    end
end

Energy = zeros(1,TNum(2));

for i=1:TNum(2)
    Energy(i)=vx(i)^2+vy(i)^2;
end

% figure
% plot(x,y);
% title('Particle Trajectory')
% xlabel('X Position')
% ylabel('Y Position')

% figure
% subplot(2,1,1)
% plot(times,x);
% hold on
% plot(times,Right);
% plot(times,Left);
% ylabel('X Position')
% title('One-Dimensional Positions of Particle and Walls')
% subplot(2,1,2)
% plot(times, y);
% hold on
% plot(times,Top);
% plot(times,Bottom);
% xlabel('Time')
% ylabel('Y Position')



figure
plot(times,Energy);
axis([0 200 0 80])
title('Particle Energy (v^2)')
xlabel('Time')
ylabel('Energy')

Answer 1

TMin = 0;
TMax = 200;
TimeStep = .0001;
times = linspace(TMin,TMax,(TMax-TMin)/TimeStep);

TMin = 0;
TMax = 201;
TimeStep = .0001;
times_2 = linspace(TMin,TMax,(TMax-TMin)/TimeStep);

all(times(1:2000000) == times_2(1:2000000))

ans =
0

因此，由于舍入误差，时间的前2,000,000个值并不完全等于times_2对应的值。差异从大约 10^-13（对于开始附近的值）到最后几个值的 10^-4 不等。

现在，由于您没有提供壁振荡频率的确切值，因此很难判断发生了什么，但我怀疑这是问题运动学的一部分。弄清楚应该很有趣。也许你可以 post 完整的代码。

在模拟弹跳粒子时，最大时间略有不同会产生截然不同的中间结果。我究竟做错了什么？

While simulating a bouncing particle, slightly different maximum times yield wildly different intermediate results. What am I doing wrong?

simulation

matlab