Strange form of increasing when derivating sine waves

2 次查看（过去 30 天）
Bruce Rogers
Bruce Rogers 2021-6-9
评论： Bruce Rogers 2021-6-11
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)
  4 个评论
显示 2更早的评论
Jan
Jan 2021-6-10
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
J. Alex Lee 2021-6-10
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2021-6-10
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
  1 个评论
Bruce Rogers
Bruce Rogers 2021-6-11
Hey Paul, that's great, thank you very much! I like the last figure, it shows the error when the drift appears

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numbers and Precision 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by