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; 

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);


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

title('Restoring force surface')

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);
        xx(i,1) = 999; % Just to give some value to be discarded later on
        fx(i,1) = 999;
xx(xx == 999) = []; % Elimination of points outside the range [-delta_v +delta_v]]
fx(fx == 999) = [];

x_sort = sort(xx);
fx_sort = sort(fx);


[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);        
        xx_dot(i,1) = 999;
        fv(i,1) = 999;

xx_dot(xx_dot == 999) = [];
fv(fv == 999) = [];

x_dot_sort = sort(xx_dot);
fv_sort = sort(fv);


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?

  • 2
    $\begingroup$ You may get more help if you change the title of your question to something more descriptive. $\endgroup$ Jan 17, 2016 at 10:56
  • $\begingroup$ As far as I remember, the best way to obtain $y=0$ or any other hyperplane from ode45 is to define event functions which correct triggers. Have a look at event location. This should give you the corresponding points in phase space with the necessary accuracy. $\endgroup$
    – Bort
    Jan 9, 2017 at 13:23

1 Answer 1


I think you do not need to use sort function and your delta_s is too high, I would suggest that reduce the delta_s to 0.000005 and remove all your sort functions and you will get right slope for damping.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.