# Error in Cumtrapz?

19 views (last 30 days)
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?

Walter Roberson on 24 Nov 2020
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.

Walter Roberson on 25 Nov 2020
When your delta_x is constant, then the code for cumulative trapazoid integration is:
my_cumtrapz = @(delta_x, y) (cumsum(y) - y/2 - y(1)/2) * delta_x
my_cumtrapz = function_handle with value:
@(delta_x,y)(cumsum(y)-y/2-y(1)/2)*delta_x
You can do experiments to confirm this numerically. For example:
deltax = rand();
y = randn(1,50);
cy = cumtrapz(deltax, y);
mcy = my_cumtrapz(deltax, y);
mcy2 = my_cumptrapz2(deltax, y);
max(abs(cy - mcy))
ans = 2.6645e-15
max(abs(cy - mcy2))
ans = 3.5527e-15
You can see that the differences between my optimized formula and a call to cumtrapz() and the long explicit way of writing the code, is down in the range of numeric round-off. Therefore, my optimized formula is correct. But if you have doubts, you can use this completely unoptimized my_cumptrapz2 code.
function cy = my_cumptrapz2(delta_x, y)
%and now an unoptimized way based upon the way the formula for trapazoid
%rule is typically expressed in terms of coefficients 2, 4, 4, 4, ... 4, 2
assert(isvector(y), 'we only handle vectors')
ny = length(y);
cy = zeros(1, ny);
if isempty(cy) || isscalar(cy); return; end
cy(1) = 0; %always
for K = 2 : ny
cy(K) = sum(addends) / 4 * delta_x;
end
end
Walter Roberson on 25 Nov 2020
Suppose that you have
A = sin(t)
then you would propose to calculate
V = int(A,t) -> -cos(t)
P = int(V,t) -> -sin(t)
So amplitude 1, frequency 1 cycle per 2*pi seconds, you would propose that the position would be strictly -sin(t)
I would then ask you: If the initial position was (say) 1 light year away from the origin, then how exactly is the object going to travel 1 light year in the infinitesimal time from 0 to dt ?
Your proposed formula does not take into account initial position or initial velocity.
KLETECH MOTORSPORTS on 29 Nov 2020
The first bit with the unoptimized formula went over my head. I've been trying to wrap my mind around it for days.
The second part, where A= sin(t), i understood the code, but i did in fact specify the initial displacement, and i took initial velocity = 0.
I am not quite able to follow the example "If the initial position was (say) 1 light year away from the origin, then how exactly is the object going to travel 1 light year in the infinitesimal time from 0 to dt ? "