Strange form of increasing when derivating sine waves

Hello everyone, I'm trying to get coordinates for a robot movement from a sine wave. Now when I change the duration of a period, the coordinates seem to drift, as shown in the picture figure 5. Can anyone explain? Or is there a mistake in my code? Thanks for your help, I'm thankful for every idea.

I'm using Matlab R2020b academic use.

clear;
b = 4;  %Factor for changing the duration of period                                                                                 
t = 0.05; % Frequenz
A = 50;  %Amplitude
x = 0:t:10*pi;
y = A*sin(b*x);
T = (1:length(x))*t;
figure(1)
plot(T,y,'.')
xlabel('time')
ylabel('distance')
%%
vel = gradient(y)./gradient(x);
figure(2)
plot(T,vel,'.')
xlabel('time')
ylabel('velocity')
%%
acc = gradient(vel)./gradient(x);
  
figure(3)
plot(T,acc,'.')
xlabel('time')
ylabel('acceleration')
acc = acc';
%xlswrite('Acceleration.xlsx',acc);
acc = acc';
%%
%velocity Cosinus Figure
v = zeros(1,length(acc));
v(1,1) = b*A * cos(0);
l_v = 1:length(acc);
for i = 2:length(acc)
    
    v(i) = ((acc(i)+acc(i-1))/2) * t + v(i-1);
    
end
figure(4)
scatter(T,v,'.')
xlabel('time')
ylabel('velocity')
v = v';
%xlswrite('Velocity.xlsx',v);
v = v';
%%
%distance Sinus Figure
s = zeros(1,length(v));
s(1,1) = b*sin(0);
for i = 2:length(v)
    
    s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);
    
end
figure(5)
scatter(T,s,'.')
xlabel('time')
ylabel('distance')
%%
s = s';
%xlswrite('idealwaves.xlsx',s)

[Scott MacKenzie]

I'm not really following the math here, but the drift disappears for the 5th plot if you remove the last term in the calculation for s. That is, by changing

s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);

to

s(i) = ((v(i)+v(i-1))/2) * t;

New Figure 5:

[Bruce Rogers]

The math is to derivate the first plot (distance) two times to get the acceleration with the analytic calculations.

acc = gradient(vel)./gradient(x);

Then I want to get back to the distance and I integrate the acceleration two times, but this time numeric with this formula:

s(1,1) = b*sin(0);
for i = 2:length(v)
    
    s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);
    
end

If I would cross out the +s(i-1) I'd only get the incremental distance, but to get the absolute, I have to add the distance from before.

So sadly, this is not a option for me, as long as you have another idea to get the absolut distance position.

[Jan]

The term (v(i)+v(i-1)) / 2 cannot be the maximum value of v. It is likely that this has more influence on the positive part than on the negative part (or vice versa). Therefore such drifts are typical for numerical integrations of accelerations or velocities. You are not the first person, who observes this.

There is no magic solution. To reduce the influence of drifting, you need absolute position values as a calibration.

[J. Alex Lee]

by the way you should be able to accomplish your numerical integration with cumtrapz, instead of manually summing

s = b*sin(0) + cumtrapz(T,v)

If you want to understand the "drift" you can systematically reduce your value of t (sample faster) to see how well your integrated position matches the analytical expectation.

Answer 1

[Paul]

I think you're seeing the effect of the discrete time step combined with numerical approximations to the differentiation and the integration. Try decreasing t. Also, I noticed few other things in the code that may be of interest.

clear;
b = 4;  %Factor for changing the duration of period                                                                                 
t = 0.05; % Frequenz  % time step
A = 50;  %Amplitude
x = 0:t:10*pi;
y = A*sin(b*x);
%T = (1:length(x))*t;
T = (0:length(x)-1)*t;  % corrected, why use T and not x, x is the independent variable?
figure(1)
plot(T,y)
xlabel('time')
ylabel('distance')
%%
%vel = gradient(y)./gradient(x);
vel = gradient(y,t);
% vel = b*A*cos(b*x);  % use this for better accuracy
figure(2)
plot(T,vel)
xlabel('time')
ylabel('velocity')
%%
% acc = gradient(vel)./gradient(x);
acc = gradient(vel,t);
% acc = -b^2*A*sin(b*x); % use this for better accuracy  
figure(3)
plot(T,acc)
xlabel('time')
ylabel('acceleration')
acc = acc';
%xlswrite('Acceleration.xlsx',acc);
acc = acc';
%%
%velocity Cosinus Figure
v = zeros(1,length(acc));
% v(1,1) = b*A * cos(0);
v(1,1) = b*A * cos(b*0);
l_v = 1:length(acc);
for i = 2:length(acc)
    
    v(i) = ((acc(i)+acc(i-1))/2) * t + v(i-1);
    
end
figure(4)
plot(T,v,T,vel,'o')
xlabel('time')
ylabel('velocity')
v = v';
%xlswrite('Velocity.xlsx',v);
v = v';
%%
%distance Sinus Figure
s = zeros(1,length(v));
% s(1,1) = b*sin(0);
s(1,1) = A*sin(b*0);  % amplitude is A
for i = 2:length(v)
    
    s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);
    
end
figure(5)
plot(T,s,T,y,'o')
xlabel('time')
ylabel('distance')
figure(6)
plot(T,s-y)  % error

[Bruce Rogers]

Hey Paul, that's great, thank you very much! I like the last figure, it shows the error when the drift appears