Error in Cumtrapz?

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?

[John]

Should you be taking into consideration the acceleration due to gravity ?

Answer 1

[Walter Roberson]

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.

[KLETECH MOTORSPORTS]

1)When you say "linear trend", exactly what do you mean?

2)How do i adjust the constant of integration? Can i evaluate it, perhaps?

3)Also, i do not understand what you said about trapz(x,y)

4) I also don't understand what you mean by ading half of the previous value aand half of the current value of y in:

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

What influence do you mean?

[Walter Roberson]

linear trend: the output you get is your expected output plus m times t, for some (negative) constant m.

[Walter Roberson]

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]

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]

(cumsum(y) - y/2 - y(1)/2) * delta_x is another way of calculating cumtrapz(x, y) with constant x differences.

[Walter Roberson]

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.

[Walter Roberson]

y = 1 : 10;
mat2str(y)

ans = '[1 2 3 4 5 6 7 8 9 10]'

y1 = cumtrapz(y);
strjoin(compose('%.4f', y1))

ans = '0.0000 1.5000 4.0000 7.5000 12.0000 17.5000 24.0000 31.5000 40.0000 49.5000'

syms Y b
Y1fun = int(Y, 1, b)

Y1fun = 

Y1 = subs(Y1fun, b, y);
strjoin(compose('%.4f', double(Y1)))

ans = '0.0000 1.5000 4.0000 7.5000 12.0000 17.5000 24.0000 31.5000 40.0000 49.5000'

%so that is an exact match between cumtrapz and a single level of int

Y2fun = int(Y1fun, 1, b)

Y2fun = 

y2 = cumtrapz(y1);
strjoin(compose('%.4f', y2))

ans = '0.0000 0.7500 3.5000 9.2500 19.0000 33.7500 54.5000 82.2500 118.0000 162.7500'

Y2 = subs(Y2fun, b, y);
strjoin(compose('%.4f', double(Y2)))

ans = '0.0000 0.6667 3.3333 9.0000 18.6667 33.3333 54.0000 81.6667 117.3333 162.0000'

hybrid_int_trap = cumtrapz(Y1);
strjoin(compose('%.4f', double(hybrid_int_trap)))

ans = '0.0000 0.7500 3.5000 9.2500 19.0000 33.7500 54.5000 82.2500 118.0000 162.7500'

%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]

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)

ans = '[0 -1882.69968752787 -3304.44883010182 -3917.15365238493 -3570.80264505458 -2350.1946141189 -554.176975744264 1377.52251389028 2971.95644934726 3838.75179539929 3765.6868207883 2770.65037952912 1097.2620955511 -844.774217771256 -2579.98034005237 -3683.51729512865 -3885.20074920062 -3135.6515587551 -1618.38550705507 295.117759939827 2136.36590670975 3454.55717115136 3926.95235820179 3437.89265071297 2107.11692163176 260.445481863336 -1649.99209521896 -3156.45406190557 -3890.10598905377 -3671.32429789221 -2553.67437669547 -810.795905576206 1130.59368072388 2795.17450314967 3775.39912208559 3831.2743642874 2949.1200217469 1344.91824372628 -588.566426022085 -2377.94950770852 -3585.12761601951 -3914.54144823213 -3285.53900951169 -1852.12203408393 34.7590103017299 1913.12983670266 3323.09975634313 3919.45895847509 3556.19791166236 2322.2555893366 519.744107246145 -1410.01885900761 -2994.56003244942 -3845.92847101535 -3755.67948843213 -2745.90918317943 -1063.84454295368 878.686344262665 2606.08416914591 3695.42169905935 3879.99111470566 3114.6033860512 1586.65212290082 -329.766916362992 -2165.44751347787 -3470.95103667109 -3926.64469243787 -3420.95878097462 -2077.70284982028 -225.752798609794 1681.46941110447 3177.0092656845 3894.70644995325 3658.84366263736 2527.16834007405 776.754069783487 -1163.83668703532 -2819.4796326444 -3784.81563139228 -3823.49676351554 -2926.05253881874 -1312.20860296923 622.909763762101 2405.51809558697 3599.17170223563 3911.62255067558 3266.37177610443 1821.39927204469 -69.515297331955 -1943.41009749099 -3341.49032698756 -3921.4571858883 -3541.31456008297 -2294.13462230628 -485.270518246932 1442.40473308026 3016.92900012568 3852.80382886305 3745.37790906367 2720.95285250435 1030.3436411037]'

y1 = cumtrapz(t, y);
strjoin(compose('%.4f', y1))

ans = '0.0000 -9.4135 -35.3492 -71.4573 -108.8970 -138.5020 -153.0239 -148.9072 -127.1598 -93.1062 -55.0840 -22.4023 -3.0628 -1.8003 -18.9241 -50.2416 -88.0852 -123.1894 -146.9596 -153.5760 -141.4186 -113.4639 -76.5564 -39.7322 -12.0071 -0.1693 -7.1170 -31.1493 -66.3821 -104.1892 -135.3142 -152.1366 -150.5376 -130.9087 -98.0559 -60.0225 -26.1205 -4.6503 -0.8686 -15.7012 -45.5165 -83.0149 -119.0153 -144.7036 -153.7904 -144.0510 -117.8698 -81.6570 -44.2787 -14.8865 -0.6765 -5.1278 -27.1507 -61.3532 -99.3612 -131.8692 -150.9179 -151.8437 -134.4199 -102.9123 -65.0353 -30.0623 -6.5560 -0.2716 -12.7477 -40.9297 -77.9176 -114.6557 -142.1490 -153.6663 -146.3877 -122.0953 -86.7367 -48.9689 -18.0389 -1.5193 -3.4547 -23.3713 -56.3927 -94.4343 -128.1821 -149.3734 -152.8199 -137.6777 -107.6543 -70.1003 -34.2103 -8.7715 -0.0120 -10.0767 -36.5012 -72.8159 -110.1298 -139.3070 -153.2040 -148.4184 -126.1217 -91.7730 -53.7821 -21.4505 -2.6940'

syms Y b T
Y1fun = int(-sym(A)*(sym(w).^2)*sin(w*T), T, 0, b)

Y1fun = 

Y1 = subs(Y1fun, b, t);
strjoin(compose('%.4f', double(Y1)))

ans = '0.0000 -9.6146 -36.1046 -72.9841 -111.2239 -141.4615 -156.2936 -152.0890 -129.8769 -95.0957 -56.2610 -22.8810 -3.1282 -1.8388 -19.3285 -51.3151 -89.9674 -125.8217 -150.0998 -156.8575 -144.4403 -115.8884 -78.1922 -40.5811 -12.2637 -0.1729 -7.2691 -31.8149 -67.8005 -106.4155 -138.2056 -155.3874 -153.7542 -133.7060 -100.1511 -61.3050 -26.6787 -4.7497 -0.8871 -16.0367 -46.4891 -84.7887 -121.5584 -147.7956 -157.0766 -147.1290 -120.3884 -83.4018 -45.2249 -15.2046 -0.6909 -5.2374 -27.7309 -62.6642 -101.4843 -134.6869 -154.1427 -155.0883 -137.2921 -105.1113 -66.4249 -30.7047 -6.6961 -0.2774 -13.0201 -41.8042 -79.5826 -117.1056 -145.1864 -156.9497 -149.5156 -124.7042 -88.5901 -50.0153 -18.4243 -1.5517 -3.5285 -23.8707 -57.5977 -96.4522 -130.9210 -152.5651 -156.0853 -140.6196 -109.9546 -71.5982 -34.9413 -8.9589 -0.0123 -10.2920 -37.2811 -74.3718 -112.4830 -142.2837 -156.4777 -151.5897 -128.8166 -93.7340 -54.9313 -21.9088 -2.7516'

%so that is NOT an exact match between cumtrapz and a single level of int
strjoin(compose('%.4f', double(Y1)-y1))

ans = '0.0000 -0.2011 -0.7553 -1.5269 -2.3269 -2.9595 -3.2698 -3.1818 -2.7171 -1.9895 -1.1770 -0.4787 -0.0654 -0.0385 -0.4044 -1.0735 -1.8822 -2.6323 -3.1402 -3.2816 -3.0218 -2.4245 -1.6358 -0.8490 -0.2566 -0.0036 -0.1521 -0.6656 -1.4184 -2.2263 -2.8914 -3.2508 -3.2166 -2.7972 -2.0952 -1.2825 -0.5581 -0.0994 -0.0186 -0.3355 -0.9726 -1.7738 -2.5431 -3.0920 -3.2861 -3.0780 -2.5186 -1.7448 -0.9461 -0.3181 -0.0145 -0.1096 -0.5801 -1.3110 -2.1231 -2.8177 -3.2248 -3.2445 -2.8722 -2.1990 -1.3897 -0.6424 -0.1401 -0.0058 -0.2724 -0.8746 -1.6649 -2.4499 -3.0374 -3.2835 -3.1280 -2.6089 -1.8534 -1.0464 -0.3854 -0.0325 -0.0738 -0.4994 -1.2050 -2.0178 -2.7390 -3.1918 -3.2654 -2.9419 -2.3003 -1.4979 -0.7310 -0.1874 -0.0003 -0.2153 -0.7799 -1.5559 -2.3532 -2.9767 -3.2736 -3.1714 -2.6949 -1.9610 -1.1492 -0.4583 -0.0576'

%and it is not even a constant difference

Y2fun = int(subs(Y1fun, b, T), T, 0, b)

Y2fun = 

y2 = cumtrapz(t, y1);
strjoin(compose('%.4f', y2))

ans = '0.0000 -0.0471 -0.2709 -0.8049 -1.7067 -2.9437 -4.4013 -5.9110 -7.2913 -8.3926 -9.1336 -9.5210 -9.6483 -9.6727 -9.7763 -10.1221 -10.8137 -11.8701 -13.2209 -14.7235 -16.1985 -17.4729 -18.4230 -19.0045 -19.2632 -19.3240 -19.3605 -19.5518 -20.0395 -20.8923 -22.0898 -23.5271 -25.0405 -26.4477 -27.5925 -28.3829 -28.8136 -28.9675 -28.9951 -29.0779 -29.3840 -30.0267 -31.0368 -32.3554 -33.8479 -35.3371 -36.6467 -37.6443 -38.2740 -38.5698 -38.6476 -38.6767 -38.8381 -39.2806 -40.0842 -41.2403 -42.6542 -44.1680 -45.5994 -46.7860 -47.6258 -48.1013 -48.2843 -48.3185 -48.3836 -48.6520 -49.2462 -50.2091 -51.4931 -52.9722 -54.4724 -55.8149 -56.8590 -57.5375 -57.8726 -57.9704 -57.9952 -58.1294 -58.5282 -59.2823 -60.3954 -61.7832 -63.2942 -64.7466 -65.9733 -66.8621 -67.3836 -67.5985 -67.6425 -67.6929 -67.9258 -68.4724 -69.3871 -70.6343 -72.0968 -73.6049 -74.9777 -76.0671 -76.7949 -77.1711 -77.2918'

Y2 = subs(Y2fun, b, t);
strjoin(compose('%.4f', double(Y2)))

ans = '0.0000 -0.0323 -0.2490 -0.7893 -1.7133 -2.9869 -4.4907 -6.0488 -7.4720 -8.6041 -9.3603 -9.7476 -9.8637 -9.8723 -9.9636 -10.3076 -11.0123 -12.0975 -13.4898 -15.0406 -16.5625 -17.8752 -18.8495 -19.4393 -19.6924 -19.7391 -19.7604 -19.9432 -20.4351 -21.3080 -22.5405 -24.0230 -25.5850 -27.0362 -28.2137 -29.0214 -29.4540 -29.5977 -29.6097 -29.6793 -29.9819 -30.6355 -31.6725 -33.0313 -34.5714 -36.1082 -37.4576 -38.4815 -39.1216 -39.4134 -39.4778 -39.4913 -39.6429 -40.0877 -40.9092 -42.0985 -43.5568 -45.1192 -46.5955 -47.8167 -48.6759 -49.1551 -49.3293 -49.3482 -49.3993 -49.6625 -50.2656 -51.2533 -52.5760 -54.1022 -55.6505 -57.0341 -58.1066 -58.7976 -59.1303 -59.2156 -59.2247 -59.3479 -59.7471 -60.5171 -61.6614 -63.0924 -64.6518 -66.1503 -67.4131 -68.3235 -68.8508 -69.0582 -69.0872 -69.1231 -69.3492 -69.9026 -70.8401 -72.1244 -73.6333 -75.1898 -76.6050 -77.7247 -78.4672 -78.8428 -78.9520'

hybrid_int_trap = cumtrapz(t, Y1);
strjoin(compose('%.4f', double(hybrid_int_trap)))

ans = '0.0000 -0.0481 -0.2767 -0.8221 -1.7432 -3.0066 -4.4954 -6.0373 -7.4471 -8.5720 -9.3287 -9.7245 -9.8545 -9.8793 -9.9852 -10.3384 -11.0448 -12.1237 -13.5034 -15.0381 -16.5446 -17.8463 -18.8167 -19.4105 -19.6748 -19.7370 -19.7742 -19.9696 -20.4677 -21.3387 -22.5618 -24.0298 -25.5755 -27.0128 -28.1821 -28.9894 -29.4293 -29.5864 -29.6146 -29.6992 -30.0119 -30.6683 -31.7000 -33.0468 -34.5711 -36.0922 -37.4297 -38.4487 -39.0918 -39.3940 -39.4735 -39.5031 -39.6679 -40.1199 -40.9407 -42.1215 -43.5657 -45.1118 -46.5737 -47.7857 -48.6434 -49.1291 -49.3161 -49.3509 -49.4174 -49.6915 -50.2985 -51.2819 -52.5934 -54.1041 -55.6364 -57.0075 -58.0740 -58.7670 -59.1092 -59.2091 -59.2345 -59.3715 -59.7788 -60.5491 -61.6859 -63.1033 -64.6466 -66.1301 -67.3830 -68.2908 -68.8235 -69.0430 -69.0878 -69.1393 -69.3772 -69.9355 -70.8697 -72.1436 -73.6374 -75.1777 -76.5797 -77.6925 -78.4358 -78.8200 -78.9433'

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

ans = '0.0000 -0.0010 -0.0058 -0.0172 -0.0365 -0.0629 -0.0940 -0.1263 -0.1558 -0.1793 -0.1952 -0.2034 -0.2062 -0.2067 -0.2089 -0.2163 -0.2311 -0.2536 -0.2825 -0.3146 -0.3461 -0.3734 -0.3937 -0.4061 -0.4116 -0.4129 -0.4137 -0.4178 -0.4282 -0.4464 -0.4720 -0.5027 -0.5351 -0.5651 -0.5896 -0.6065 -0.6157 -0.6190 -0.6196 -0.6213 -0.6279 -0.6416 -0.6632 -0.6914 -0.7233 -0.7551 -0.7831 -0.8044 -0.8178 -0.8241 -0.8258 -0.8264 -0.8299 -0.8393 -0.8565 -0.8812 -0.9114 -0.9438 -0.9744 -0.9997 -1.0177 -1.0278 -1.0317 -1.0325 -1.0338 -1.0396 -1.0523 -1.0729 -1.1003 -1.1319 -1.1639 -1.1926 -1.2149 -1.2294 -1.2366 -1.2387 -1.2392 -1.2421 -1.2506 -1.2667 -1.2905 -1.3202 -1.3524 -1.3835 -1.4097 -1.4287 -1.4398 -1.4444 -1.4454 -1.4464 -1.4514 -1.4631 -1.4826 -1.5093 -1.5405 -1.5728 -1.6021 -1.6254 -1.6409 -1.6490 -1.6515'

%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'})

[Walter Roberson]

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.

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)

P = 

rewrite(P, 'cos')

ans = 

Q = int(P, T, 0, T)

ans = 

and there the double integral is in the expected sin() form, but with a linear dependence on time.

[Paul]

"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]

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]

Thanks Walter, for that in depth explanation, but in all honesty, i think it would be simpler to find another way to double integrate my ode. I'm quite new to matlab, so the simpler the method, the easier, at least for now.

Do you have any suggestions as to how i could integrate it another way? perhaps if i write a function with the actual steps involved in tarpezoidal integration?

Thank you though for all the time you put into answering my questions.

[Walter Roberson]

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
       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]

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]

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 ? "