Моделирование отскока шарика в дискретном времени

[Tatarin]

Привет,
Я написал код в матлабе чтобы смоделировать отскок шарика при ударении о поверхность. Начальная высота падения 20м. Вес шарика 0.1кг. Радиус шарика 2м. Упругость я смоделировал как пружину которая сжимается при позиции когда высота pt меньше радиуса r шарика. Почему то график высоты не выглыдит как а полусинусоида, а выглядит как пилообразная кривая. Может кто поскажет где я ошибся с физикой? Трение о воздух или нагрев от удара я не учитываю, т.к. хотел отладить простую модель сначала.

Code: Select all

clear all;

% Bouncing ball definition
r=2;      % ball radius in m
k=500;    % spring stiffness
m=0.1;      % mass of the ball in kg
dt=0.1;   % time increment in seconds
N=1000; % simulation length
pt=zeros(1,N); % position at time t in m
pt(1)=20;      % initial position is 20m
at=zeros(1,N); % acceleration at time t m^2/sec
at(1)=-9.8;    %  initial acceleration

for j=2:N
    pt(j)=pt(j-1)+(0.5*dt^2)*at(j-1); % integration of a to get location
    if pt(j)>r
        at(j)=-9.8;
    else
        at(j)=at(j-1)+k*(r-pt(j-1))/m; 
    end;
    
end

%plot(at);
plot(pt);

[perasperaadastra]

Первая часть (pt>r) у вас правильная: delta pt = 0.5*a*t^2
Проблема со второй частью, где at меняется, и интеграция усложняется.
dp/dt = v
dv/dt = a
a = функция позиции

Кто тут дифуры численно решать любит?

[Tatarin]

Я переписал решение дифура и все заработало
График приложил

Code: Select all

for i=2:N
    t(i)=(i-1)*dt;
    y(i)=y(i-1)+dt*z(i-1);
    
    if y(i-1)>r
        z(i)=z(i-1)+dt*(-G);
    else
        z(i)=z(i-1)+dt*(K*(R-y(i-1))-G);
    end
    
end

otskok.jpg