Value of trapz shrinking to zero

Question

Allan Andre do Nascimento 2016-8-1

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/297954-value-of-trapz-shrinking-to-zero

评论： Allan Andre do Nascimento 2016-8-1

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?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Allan Andre do Nascimento 2016-8-1

在 MATLAB Online 中打开

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!

John D'Errico 2016-8-1

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 2016-8-1

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

请先登录，再进行评论。

请先登录，再回答此问题。