Why is there a diiference in state-space velocity and derivative of state space displacement?

1 view (last 30 days)
I stumbled on this matlab code online from mathworks. However, I am surprised that the velocity (yvector(:,1)) is different from the velocity obtained from diff(yvector(:,2))./diff(t). Why des this occur? https://www.mathworks.com/matlabcentral/answers/1613980-how-can-i-obtain-the-derivative-of-a-vector-displacement-velocity-but-what-about-acceleration-i-am
clc; %command window
clear;%clear workspace
close all; % close figures
global m k
m = 80; %Defining mass of SDOF system
k = 10000;%Defining stiffness of SDOF system
dt = .02; %Defining Time Step:
t = 0:dt:4;% Defining time vector.
y0 = [0.085; 0.1]; %initial vel and disp [vel disp] %velocity.
[tsol, yvectorl] = ode45(@testode1,t,y0);
%%plot
figure()
plot(tsol,yvectorl(:,2)), title("Displacement") %Disp.(2)
figure()
plot(tsol,yvectorl(:,1)), title("Velocity") %Vel. (1)
% THIS NEEDS TO GO AT THE END OF THE FILE
%%function statespace
function dy = testode1(t,y)
global m k
dy=[-k*y(2)/m;y(1)];
end

Answers (1)

Mathieu NOE
Mathieu NOE on 11 Jan 2023
hello
where do you see a difference ?
diff or gradient based dispkacement derivatives do pretty well overlay with computed velocity
clc; %command window
clear;%clear workspace
close all; % close figures
global m k
m = 80; %Defining mass of SDOF system
k = 10000;%Defining stiffness of SDOF system
dt = .02; %Defining Time Step:
t = (0:dt:4)';% Defining time vector.
y0 = [0.085; 0.1]; %initial vel and disp [vel disp] %velocity.
[tsol, yvectorl] = ode45(@testode1,t,y0);
%%plot
figure()
plot(tsol,yvectorl(:,2)), title("Displacement") %Disp.(2)
% derivative of disp
% ddisp = [0;diff(yvectorl(:,2))./diff(t)];
ddisp = gradient(yvectorl(:,2),dt);
figure()
plot(tsol,yvectorl(:,1),tsol,ddisp,'*r'), title("Velocity") %Vel. (1)
% THIS NEEDS TO GO AT THE END OF THE FILE
%%function statespace
function dy = testode1(t,y)
global m k
dy=[-k*y(2)/m;y(1)];
end
  4 Comments
Mathieu NOE
Mathieu NOE on 11 Jan 2023
Here we can test 3 methods
the 3rd option is the best and gives less than 1% error
also if you further reduce dt by factor 5 you can reduce that error to less than 0.1%
the difference between true and estimated velocity is plotted in the second subplot of fig 2
clc; %command window
clear;%clear workspace
close all; % close figures
global m k
m = 80; %Defining mass of SDOF system
k = 10000;%Defining stiffness of SDOF system
dt = .02/5; %Defining Time Step:
t = (0:dt:4)';% Defining time vector.
y0 = [0.085; 0.1]; %initial vel and disp [vel disp] %velocity.
[tsol, yvectorl] = ode45(@testode1,t,y0);
%%plot
figure(1)
plot(tsol,yvectorl(:,2)), title("Displacement") %Disp.(2)
% derivative of disp
% ddisp = [0;diff(yvectorl(:,2))./diff(t)]; % worst
% ddisp = gradient(yvectorl(:,2),dt); % better in average, except first and last point
[ddisp, ~] = firstsecondderivatives(t,yvectorl(:,2)); % best
figure(2)
subplot(2,1,1),plot(tsol,yvectorl(:,1),tsol,ddisp,'*r'), title("Velocity") %Vel. (1)
subplot(2,1,2),plot(tsol,yvectorl(:,1)-ddisp(:)) %Vel. (1)
% THIS NEEDS TO GO AT THE END OF THE FILE
%%function statespace
function dy = testode1(t,y)
global m k
dy=[-k*y(2)/m;y(1)];
end
%%%%%%%%%%%%%%%%%%
function [dy, ddy] = firstsecondderivatives(x,y)
% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.
n = length (x);
dy = zeros;
ddy = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1) = (-3*y(1) + 4*y(2) - y(3)) / (2*(x(2) - x(1))); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i = 2:n-1
dy(i) = (y(i+1) - y(i-1)) / (x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / (2*(x(n) - x(n-1)));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end

Sign in to comment.

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!