I am trying to apply the "restoring force surface" method to a dynamic linear system. The idea behind this method is that, knowing acceleration, displacement, velocity and input force it is possible to establish if the system under study is linear or not by plotting a surface in the phase-space.
To verify my code I applied the method to a system I already know, therefore I generated the needed data ($\ddot{y},\dot{y},y$) solving the following equation:
$ \ddot{y} + 40 \dot{y} + 10^{4} y = f(t) $
Where $f(t) = 30 \sin(10 t)$.
This is the code I wrote to do that:
[t,u] = ode45('odefun',[0 60],[0 0]);
vel = u(:,2);
spost = u(:,1);
f_input = 30*sin(10*t);
acc = zeros(size(vel)); % obtain the acceleration from the velocity
acc(1) = 0;
for i = 2:length(vel)
h = t(i)-t(i-1);
acc(i) = (vel(i)-vel(i-1))/h;
end
Where the $odefun$ is:
function [F] = odefun(t,y)
F = zeros(size(y));
F(1) = y(2);
F(2) = -10^4*y(1)-40*y(2)+30*sin(10*t);
end
For the application of the method I have to compute the so called restoring force as
$ z(y,\dot{y}) = f(t) - \ddot{y} $
therefore I can plot the triplet $(y,\dot{y},z)$ as a surface where $y$ and $\dot{y}$ represent a point in the phase-space and $z$ the height above that point.
% Surface
x_plot = linspace(min(spost),max(spost),200);
y_plot = linspace(min(vel),max(vel),200);
zgrid = gridfit(sort(spost),sort(vel),sort(z), x_plot,y_plot);
% trisurf should work as well instead of gridfit
% sort() used not to have oscillatory data
figure()
surf(x_plot,y_plot,zgrid)
title('Restoring force surface')
xlabel('displacement')
ylabel('velocity')
zlabel('force')
Once I have obtained the surface I need to extract from it the cross sections: one for $y = 0$ and another one for $\dot{y} = 0$. The slope of the cross section at $\dot{y} = 0$ should represent the stiffness of the system (therefore it should be $k = 10^4$), whereas the slope of the cross section obatined at $y = 0$ should represent the damping ($ = 40$).
But here I have problems, infact, since the system is known to be linear (the coefficients in the equation are constants), I expect to have straight lines as cross sections whereas the plot of $velocity-force$ is not a straight line as it should be.
To exctract a cross-section I saved the values of $y$ and $z$ corresponding to $|\dot{y}| < \delta$
% Surface section to plot force-displacement
delta_v = 0.0001;
xx = zeros(length(spost),1);
fx = zeros(length(z),1);
for i = 1:length(vel)
if abs(vel(i)) <= delta_v
xx(i,1) = spost(i);
fx(i,1) = z(i);
else
xx(i,1) = 999; % Just to give some value to be discarded later on
fx(i,1) = 999;
end
end
xx(xx == 999) = []; % Elimination of points outside the range [-delta_v +delta_v]]
fx(fx == 999) = [];
x_sort = sort(xx);
fx_sort = sort(fx);
figure()
plot(x_sort,fx_sort)
xlabel('displacement')
ylabel('force')
[b,bint] = polyfit(x_sort,fx_sort,1);
err_K = 10^4-b(1); % It should be 0
The same for the $y = 0$:
delta_s = 0.01;
xx_dot = zeros(length(vel),1);
fv = zeros(length(z),1);
for i = 1:length(spost)
if abs(spost(i)) <= delta_s
xx_dot(i,1) = vel(i);
fv(i,1) = z(i);
else
xx_dot(i,1) = 999;
fv(i,1) = 999;
end
end
xx_dot(xx_dot == 999) = [];
fv(fv == 999) = [];
x_dot_sort = sort(xx_dot);
fv_sort = sort(fv);
figure()
plot(x_dot_sort,fv_sort)
xlabel('velocity')
ylabel('force')
This last plot is not even a straight line and this means the damping is not linear but I know it must be linear since I started from a linear equation.
Any idea?
