我正在尝试将“恢复力面”方法应用于动态线性系统。这种方法背后的想法是，在知道加速度、位移、速度和输入力的情况下，可以通过在相空间中绘制一个表面来确定所研究的系统是否是线性的。

为了验证我的代码，我将该方法应用于我已经知道的系统，因此我生成了所需的数据（$\ddot{y},\dot{y},y$) 求解以下方程：

$\ddot{y} + 40 \dot{y} + 10^{4} y = f(t)$

在哪里$f(t) = 30 \sin(10 t)$.

这是我为此编写的代码：

[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

在哪里$odefun$是：

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

对于该方法的应用，我必须将所谓的恢复力计算为

$z(y,\dot{y}) = f(t) - \ddot{y}$

因此我可以绘制三元组$(y,\dot{y},z)$作为一个表面$y$和$\dot{y}$表示相空间中的一个点，并且$z$高于该点的高度。

% 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')    

一旦我获得了表面，我需要从中提取横截面：一个用于$y = 0$另一个用于$\dot{y} = 0$. 横截面的斜率为$\dot{y} = 0$应该代表系统的刚度（因此它应该是$k = 10^4$)，而横截面的斜率在$y = 0$应该代表阻尼（$= 40$）。

但这里我有问题，事实上，由于系统已知是线性的（方程中的系数是常数），我希望有直线作为横截面，而$velocity-force$不是一条直线。

为了提取横截面，我保存了$y$和$z$对应于$|\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

相同的$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')

最后一个图甚至不是一条直线，这意味着阻尼不是线性的，但我知道它必须是线性的，因为我从线性方程开始。

任何想法？