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

1 visualización (últimos 30 días)
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

Respuestas (1)

Mathieu NOE
Mathieu NOE el 11 de En. de 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 comentarios
Mathieu NOE
Mathieu NOE el 11 de En. de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by