MATLAB: Error in Cumtrapz

cumtrapzerrorinitial conditionsintegrationMATLABnonlinearnumerical integrationodependulumtrapz

Hi, i'm trying to integrate acceleration of a simple pendulum twice to get the displacement, and i've plotted it against time.
i.e,
On the same graph, i've plotted the displacement versus time
i.e,
with certain initial conditions as follows:
%this part integrates the acceleration(y1 = d^2(x)/dt^2) twice
% to get both velocity(y2) and displacement(y3)
t=0:0.01:1    
w=50;      
A=pi/2;       
y1=-A*(w.^2)*sin(w*t);   %here y1 is diff(x,t,2) , i.e, accln
y2=cumtrapz(t',y1');     %here y2 is diff(x,t) , i.e, velocity
y3=cumtrapz(t',y2);      %here y3 is the double integrated acceleration
y3a= A*sin(w*t);         %direct displacement expression
plot(t,y3, t,y3a)
legend('y3', 'displ')
There is an error in the blue curve obtained by double integration, as both curves should be similar.
I think it might have to do with an error in trapezoidal integration.
Any idea what the error is and how i can correct it?

Best Answer

Remember after one integration you get a formula with a constant of integration, and when you integrate that, you get the double integral plus the integral of the constant of integration, plus a second constant of integration. And the integral of the first constant of integration is the constant times the variable of integration.
Your plot shows a linear trend. That is the constant of integration multiplied by time.
Your blue plot is therefore one of the mathematically correct solutions.
If it is not the solution you want then you need to adjust the constant of integration.
Note that mathematically, trapz(x, y) with equidistant x, is equivalent to (sum(y) - (y(1)+y(end))/2) * delta x
cumtrapz with constant x difference effectively proceeds step by step, taking the proceeding result and adding half of the previous y value and half of the current y value. Those halfs have an influence.
Related Question