# 2D cross section from 3D surface

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?

• You may get more help if you change the title of your question to something more descriptive. – David Ketcheson Jan 17 '16 at 10:56
• 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. – Bort Jan 9 '17 at 13:23

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.