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

3 ビュー (過去 30 日間)
Joe
Joe 2023 年 1 月 11 日
コメント済み: Mathieu NOE 2023 年 1 月 31 日
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

回答 (1 件)

Mathieu NOE
Mathieu NOE 2023 年 1 月 11 日
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 件のコメント
Mathieu NOE
Mathieu NOE 2023 年 1 月 11 日
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

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeGraphics Object Programming についてさらに検索

タグ

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by