Value of trapz shrinking to zero

Dear Matlab community,

I am calculating an index with trapz and I am having some strange results. This index gets a solution from an ode and makes use of trapz. The strange thing is that depending on the number of points generated by the ode, the trapz value decreases indefinetly (the bigger the number of points the smaller the value of trapz). A decrease could be expected if the trapezoids approximating the function were "bigger" that the real functon, but as the number of points increased, trapz value should approximate a fixed value (which is the real value of trapz). Furthermore trapz is calculating the area based only on positive numbers, thus the area should not go to zero.... Could someone please help me with that?

[Steven Lord]

Please post a SMALL sample of code that people can run that demonstrates the behavior you're describing.

[John D'Errico]

Ok I'll help. You have made a mistake. Or perhaps you misunderstand how to use trapz or how to solve whatever problem you want to solve.

The thing is, there are infinitely many ways you could have done something wrong. We cannot guess which of them it is, unless you show us what you did. The crystal ball is so foggy today.

So please attach your code to your question, or to a comment (NOT to an answer.) Make it easy to let us help you.

[Allan Andre do Nascimento]

This is the code:

function [t,x]=call_passive_quarter()
%%%%%%%%%%%%%%%%%Parameter and Input Definition %%%%%%%%%%%%%%%%%%%%%%%

%%%If you have changed anything on passive quarter, do not forget to
%%%change it here as well (be it input function or variable values such as
%%%V and amp

V=30;lam=4;
ome=2*pi*(V/(lam*3.6));
amp=0.025;

%%%%%%%%%%%%%%%Constant declaration and ode15 simulation %%%%%%%%%%%%%%%%%
k2 = 190000;
g=9.81;
m1 = 290;
m2 = 59;
x0 = [0 0 0 0]';

tspan=[0 5]; 
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,x] = ode15s(@passive_quarter, tspan, x0,opts); % Had to use ode 45 because ode15s does not get enough points

%%%Bump input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp above as desired 

%w=amp*(1-cos(8*pi*t)).*((heaviside(t-0.5)-heaviside(t-0.75))+(heaviside(t-1.5)-heaviside(t-1.75)));

%%%Sinusoidal input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp
%%%and V above as desired 

%w = amp*(heaviside(t)-heaviside(t-1.2)).*(sin(ome*t));

%%%Random input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp
%%%above and the value of V below as desired

w=amp*(-heaviside(t-0.1)+heaviside(t-1.6));%.*(1+0.1*randn(length(t),1));

%%%%%%%%%%%%%%%%%%%Renaming Variables %%%%%%%%%%%%%%%%%%%%%%%%%%
x1=x(:,1);
x1dot=x(:,2);
x2=x(:,3);
x2dot=x(:,4);
tmp=t;
%%%%%%%%%%%%%Trajectory and Suspension Deflection Plot %%%%%%%%%%%%%%%%%
figure(1);
subplot(2,1,1);hold on; plot(t,x1,'b');plot(t,w,'r');legend('Sprung Mass','Road Input');title('Sprung Mass Displacement x Road Input');xlabel('Time (s)');ylabel('Displacement Xs (m)');
subplot(2,1,2);hold on;plot(t,x2,'b');plot(t,w,'r');legend('Unsprung Mass','Road Input');title('Unprung Mass Displacement x Road Input');xlabel('Time (s)');ylabel('Displacement Xu (m)');
figure(2);plot(t,x1-x2);title('Suspension displacement');xlabel('Time (s)');ylabel('Displacement Xs-Xu (m)');

%%%%%%%%%%%%%%%%%%Calculating Comfort Index %%%%%%%%%%%%%%%%%%%%%%%%%%%
ttest=0:0.0001:5; % numero de pontos influencia no resultado. Logo decidiu-se pegar o mesmo numero de pontos de interpolaçao em todos os casos para um resultado mais justo.
xstest = interp1(t,x1(:,1),ttest);
xutest = interp1(t,x2(:,1),ttest);
wtest  = interp1(t,w(:,1),ttest);
acci=diff(diff(xstest));
%acci=diff(diff((x(:,1))'));
quos=(acci.^2);
tempo2= ttest(1,2:(length(ttest)-1));
normacc2=(10^9)*sqrt((1/max(t))*(trapz(tempo2,quos))); % Trapz is a numerical integration command
normacc3=(10^10)*sqrt(sum(quos(:)));
%tempo2= tmp(2:(length(tmp)-1),1)';
%normacc2=(10^9)*sqrt((1/max(t))*(trapz(tempo2,quos)));
%%%%%%%%%%%%%%%%%%%Calculating Safety Index %%%%%%%%%%%%%%%%%%%%%%%%
Fdinp=(k2*(wtest-xutest));
Fstat=((m1+m2)*g);
Fdinpnumi=(Fdinp/Fstat).^2;
Fdinpnum=trapz(ttest,Fdinpnumi); % Numerical diferentiation loses 2 points, thus here ttest must be used
Fdin=sqrt(Fdinpnum/(max(t)));

%%%%%%%%%%%%%%%%%%%% Plotting Indexes  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
subplot(2,1,1);plot(max(t),normacc2,'o');title('Comfort Index');
subplot(2,1,2);plot(max(t),Fdin,'o');title('Safety Index');

%%%confirm values
figure(4)
plot(tempo2,acci)
figure(5)
plot(ttest,(Fdinp/Fstat))

which calls from this code:

function dxdt = passive_quarter(t,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Parameter Definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
amp=0.025;
%%%Bump input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp above as desired 

%w=amp*(1-cos(8*pi*t)).*((heaviside(t-0.5)-heaviside(t-0.75))+(heaviside(t-1.5)-heaviside(t-1.75)));

%%%Sinusoidal input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp
%%%above and the value of V below as desired 

V=30;lam=4;
ome=2*pi*(V/(lam*3.6));
%w = amp*(heaviside(t)-heaviside(t-1.2)).*(sin(ome*t));

%%%Random input - uncomment to use it as input or if it is uncommented, this
%%%is the input that will be utilized for the current simulation. Change the value of amp
%%%above and the value of V below as desired 

w=amp*(-heaviside(t-0.1)+heaviside(t-1.6));%.*(1+0.1*randn(length(t),1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Constant Definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=9.81;
m1 = 290;
m2 = 59;
k1 = 16812;
k2 = 190000;
b1 = 1000;
alfa=4.515*(10^13);
beta=1;
gama=1.545*(10^9);
tau=(1/30);
Ps=10342500;
A=3.35*(10^(-4));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%State Space Equations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dxdt=zeros(4,1);
dxdt(1) = x(2);
dxdt(2) = -(k1/m1)*x(1)-(b1/m1)*x(2)+(k1/m1)*x(3)+(b1/m1)*x(4);
dxdt(3) = x(4);
dxdt(4) = (k1/m2)*x(1)+(b1/m2)*x(2)-((k1+k2)/m2)*x(3)-(b1/m2)*x(4)+(k2/m2)*w;

Thank you very much for helping!

[Allan Andre do Nascimento]

I think I got why it is behaving this way. I am interpolating between points which are, in general, approximate points. As ode returns more points in the solution, such answer is a more accurate representation of the real graph. As I interpolate between those more accurate points and the acceleration is very small, it gets indeed closer to zero...

[John D'Errico]

This is always the best way to solve a problem - by thinking it out, making a conjecture about the problem. When I have a problem like this, once I have some conjecture as to the reason, there is usually a way I can test the validity of my conjecture.

For example, use a finer interval. Or use a better interpolant. The default for interp1 is a piecewise linear interpolation. However you can also use various spline options for interp1, which will generally yield a higher accuracy to approximate a smooth function.

[Allan Andre do Nascimento]

Thank you very much for the help John! I am thinking about defining the internal step in tspan. Like this, I will have the internal answers calculated by ode in a pre-defined interval and number of points, which can give me a better accuracy than interp1 ,which is just a linear interpolation performed from the internal data provided by ode (please correct me if I am wrong).