Error in Cumtrapz?
Show older comments
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?
1 Comment
John
on 19 Oct 2022
Should you be taking into consideration the acceleration due to gravity ?
Accepted Answer
Walter Roberson
on 24 Nov 2020
0 votes
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.
15 Comments
KLETECH MOTORSPORTS
on 24 Nov 2020
Walter Roberson
on 24 Nov 2020
linear trend: the output you get is your expected output plus m times t, for some (negative) constant m.
Walter Roberson
on 24 Nov 2020
The formula for trapezoid rule with constant differences in x is
(2*f(x1) + 4*f(x2) + 4*f(x3) + ... + 4*f(x(n-1)) + 2*f(xn))/4 * delta_x
Divide through by 4 to get
(f(x1)/2 + f(x2) + f(x3) + ... + f(x(n-1)) + f(xn)/2) * delta_x
But f(x1)/2 = f(x1)-f(x1)/2 and f(xn)/2 = f(xn)- f(xn)/2 so this can be rewritten as
(f(x1)+f(x2)+f(x3)...+f(xn)- f(x1)/2-f(x2)/2) * delta_x
which is (sum(f(x)) -(f(x1)+f(xn))/2) * delta_x
or (sum(y) - (y(1)+y(end))/2) * delta_x
Walter Roberson
on 24 Nov 2020
Now cumtrapz is the sequence
[trapz(x(1:1),y(1:1)), trapz(x(1:2),y(1:2)), trapz(x(1:3),y(1:3))... trapz(x(1:end), y(1:end))]
consider two adjacent entries:
(sum(y(1:n)) - y(1)/2 - y(n)/2) * delta_x
(sum(y(1:n+1)) - y(1)/2 - y(n+1)/2) * delta_x
subtract to find the difference:
(sum(y(1:n+1)) - sum(y(1:n)) - y(1)/2 - (-y(1)/2) - y(n+1)/2 - (-y(n)/2)) * delta_x
which is
(y(n+1) - 0 - y(n+1)/2 + y(n) /2) * delta_x
which is
(y(n+1)/2 + y(n)/2) * delta_x
So if you know cumtrapz up to a particular point, you can find the next value by adding half of the previous y and half of the current y, all times delta_x, as I indicated earlier.
Walter Roberson
on 24 Nov 2020
(cumsum(y) - y/2 - y(1)/2) * delta_x is another way of calculating cumtrapz(x, y) with constant x differences.
Walter Roberson
on 24 Nov 2020
You are expecting cumtrapz of cumtrapz to be like double integral, but you can see from the above that the second cumtrapz is going to be mathematically subtracting off values.
What does this all mean? It means that cumtrapz of cumtrapz is mathematically not the same as a double integral that assumes that the constants of integration are 0.
y = 1 : 10;
mat2str(y)
y1 = cumtrapz(y);
strjoin(compose('%.4f', y1))
syms Y b
Y1fun = int(Y, 1, b)
Y1 = subs(Y1fun, b, y);
strjoin(compose('%.4f', double(Y1)))
%so that is an exact match between cumtrapz and a single level of int
Y2fun = int(Y1fun, 1, b)
y2 = cumtrapz(y1);
strjoin(compose('%.4f', y2))
Y2 = subs(Y2fun, b, y);
strjoin(compose('%.4f', double(Y2)))
hybrid_int_trap = cumtrapz(Y1);
strjoin(compose('%.4f', double(hybrid_int_trap)))
%but not an exact match between second cumtrapz versus double integral
%but there is an exact match between second cumtrapz and cumtrapz of integral
plot(y, y2, 'k*-', y, Y2, 'b^-', y, hybrid_int_trap, 'gv-' );
legend({'traptrap', 'intint', 'inttrap'})
Walter Roberson
on 24 Nov 2020
Edited: Walter Roberson
on 24 Nov 2020
I am still working on it; here are some investigations
t = 0:0.01:1;
w = 50;
A = pi/2;
y = -A*(w.^2)*sin(w*t); %here y1 is diff(x,t,2) , i.e, accln
mat2str(y)
y1 = cumtrapz(t, y);
strjoin(compose('%.4f', y1))
syms Y b T
Y1fun = int(-sym(A)*(sym(w).^2)*sin(w*T), T, 0, b)
Y1 = subs(Y1fun, b, t);
strjoin(compose('%.4f', double(Y1)))
%so that is NOT an exact match between cumtrapz and a single level of int
strjoin(compose('%.4f', double(Y1)-y1))
%and it is not even a constant difference
Y2fun = int(subs(Y1fun, b, T), T, 0, b)
y2 = cumtrapz(t, y1);
strjoin(compose('%.4f', y2))
Y2 = subs(Y2fun, b, t);
strjoin(compose('%.4f', double(Y2)))
hybrid_int_trap = cumtrapz(t, Y1);
strjoin(compose('%.4f', double(hybrid_int_trap)))
%not an exact match between second cumtrapz versus double integral
%and not an exact match between second cumtrapz and cumtrapz of integral
strjoin(compose('%.4f', double(hybrid_int_trap)-y2))
%that is a growing difference between the two but comparatively small to the values
plot(t, y2, 'k*-', t, Y2, 'b^-', t, hybrid_int_trap, 'gv-' );
legend({'traptrap', 'intint', 'inttrap'})
So, where are we? Well look at Y2fun above, and notice that it is not 
and is instead that minus 
where b is the current time . The true double integral reflects that the first integral introduces a boundary condition, and the true double integral might not match exactly with the double cumtrapz but you can see from above that the difference is small compared to the value of the integral. The double cumtrapz() is thus a comparatively good approximation of the double integral.
and is instead that minus where b is the current time . The true double integral reflects that the first integral introduces a boundary condition, and the true double integral might not match exactly with the double cumtrapz but you can see from above that the difference is small compared to the value of the integral. The double cumtrapz() is thus a comparatively good approximation of the double integral.The error in your logic was in neglecting that constant of integration.
You might notice that the symbolic integral of the original sin() was in terms of sin()^2, but that has to do with formulas for cos(theta/2) and has to do with sin^2 + cos^2 = 1 so cos^2 can be rewritten in terms of sin^2
syms w t T
P = int(sin(w*t), t, 0, T)
rewrite(P, 'cos')
Q = int(P, T, 0, T)
and there the double integral is in the expected sin() form, but with a linear dependence on time.
Paul
on 25 Nov 2020
"The error in your logic was in neglecting that constant of integration."
Are you saying this was the only error? Or are you saying that there is an additional issue such that double cumtrapz is not a valid approximation a double integral for a fundamental reason (rather than a numerical reason, such as accumulation of error)?
Walter Roberson
on 25 Nov 2020
Double cumtrapz approaches the double integral as the distance between points gets smaller, except possibly for a boundary value adjustment. For continuous finite functions the boundary differences should tend to zero, but if the function has discontinuities such as Dirac, or Heaviside, or derivative of abs(), then the boundary differences could potentially be substantial.
However when I say double integral here the integral of the constant of integration must be included.
KLETECH MOTORSPORTS
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
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))
max(abs(cy - mcy2))
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
addends = zeros(1, K);
addends(1) = y(1) * 2;
addends(2:K-1) = y(2:K-1) * 4;
addends(K) = y(K) * 2;
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
More Answers (0)
Categories
Find more on Numerical Integration and Differentiation in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
