Cumtrapz integration drift with time

Hello, I'm trying to get a displacement from a velocity. For a constant velocity, the formula is:

x_h = x_c + V*t 

Where x_c is an initial value, V is the velocity magnitude and t is a time vector. Please note: the time vector t is linearly increasing. I now have a velocity u that is as below.

Theta is an azimuthal position, it is essentially related to time.

It starts from an undisturbed value (equal to V) then drops and recovers twice and at the end of the period (theta = 360 deg) it goes back to the undisturbed value. To compute the displacement, I do:

x_h = x_c + cumtrapz(u,t) 

Below the displacements for constant velocity (formula with no integration) and non-constant velocity are plotted:

The issue is that, as you can see, cumtrapz drifts more and more away from the undisturbed value as time goes by. Consider that there's also another period preceding theta = 0 deg, so x_h is not starting at the undisturbed value at theta=0 deg, although it should. This would only happen if there was no drift introduced by cumtrapz. If I try to do:

x_h(t) = x_h(t-1) + u(t-1) dt(t-1) 

I get basically the same. I checked that x_h(t) - x_h(t-1) behaves as the above-plotted u velocity, which is the case, so both the integral formulas are correct. However, I get the above-mentioned drift issue. Please note: u velocity is supposed to behave like that, its trend is not due to some sort of bias or noise in the data. How can I fix this, preferably keeping x_h(t) - x_h(t-1) behaving as the u velocity? Thanks for your help.

[David Goodmanson]

Hi Andrea,

This behavior is expected. Without the complications of angles, suppose you have an object moving at v1 m/sec but there is a portion tau seconds long where the object is is moving at a slower velocity v2 m/sec. The distance traveled in that portion will be v2*tau rather than v1*tau, so there will be a shortfall in distance of (v2-v1)*tau, and after the object is once again traveling at velocity v1, the shortfall will persist in a plot of distance vs time. That's the effect you are seeing with the red line dropping down between 200 and 300 degrees and then staying below the blue line by what should be a constant distance after that.

[Andrea Improta]

Hi @David Goodmanson, thank you for your answer. So you're saying it is something purely mathematical, and reducing the time step of the integration won't fix it.

Answer 1

[John D'Errico]

This is expected behavior. And reducing the time step WILL reduce the problem. And it IS something purely mathematical, yes, you have that right.

How does cumtrapz work? It is purely a cumulative trapezoidal integration. Put trapezoidal pieces under your function. Sum of the areas in a cumulative way.

Look at a simple example. Integrate a function with positive curvature. For example, x^2, from 0 to 1. We know the answer already, the integral will be (done in MATLAB, since this is a MATLAB forum)

syms x
int(x^2), disp(char(ans))

x^3/3

Yeah, I know. calc 101. But this way we have some MATLAB content in this answer. ;-)

But now, look at what cumtrapz does.

X = 0:.5:1;
Y = @(x) x.^2;
fplot(Y,[0,1],'k')
hold on
plot(X,Y(X),'-ro')
hold off
grid on

So a trapezoidal integration tool will overestimate the area. And cumtrapz, since it just sums the pieces, will do the same.

vpa([int(x^2,[0,1/2]),int(x^2,[0,1])],5), disp(char(ans))

[0.041667, 0.33333]

cumtrapz(X,Y(X))

ans = 1×3

0 0.0625 0.3750

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So instead of 1/3 for the integral, we will see 0.375 from trapz. And each segment over that domain will similarly over-estimate the area. We see that from cumtrapz.

Next, plot the integral, comparing to to the known value.

fplot(x^3/3,[0,1])
hold on
plot(X,cumtrapz(X,Y(X)),'ro')
hold off
grid on

That is a very coarse interval, but we see a steady and consistent bias arises.

But what happens when we reduce the interval?

X = 0:.25:1;
fplot(Y,[0,1],'k')
hold on
plot(X,Y(X),'-ro')
hold off
grid on

Do you see the two curves are now much closer? So the trapezoidal segments fit much more tightly against the curve. This is of course expected. And of course, the integral a trapezoidal rule wil yield is now much closer, but below we see the red circles always lie above the true curve.

fplot(x^3/3,[0,1])
hold on
plot(X,cumtrapz(X,Y(X)),'ro')
hold off
grid on

You should see the decrease in the bias, even though I only cut the interval in half. It will ALWAYS over-estimate the integral of this particular function, even if we used a very tiny increment.

trapz(X,Y(X))

ans = 0.3438

X = linspace(0,1,10000);
format long G
trapz(X,Y(X))

ans = 
         0.333333335000333

Again, the over-estimate arises because of the positive curvature of the integrand. Again, this is completely expected, as is the reduction in the bias if we take a smaller time step.

Finally, you can show this is the expected behaviour for a trapezoidal integration tool, as is cumtrapz, when applied to such an integrand. But that would begin to turn into a class in numerical analysis.

[Andrea Improta]

Wow, what a complete answer! Thank you very much, John