The accuracy of cumtrapz (numerical integration)????

I have a problem where I have to integrate a vector to get another vector. (this problem related to noncircular gear design).

I am trying to find phi 2 the output angle depending on the integration in equation 2.4.5.

I know that the last value of phi 2 should be pi=3.14 almost. but when I am using cumtrapz function as in the code down, the values of phi 2 is not ok.

I think that because of the approximation of the function (cumtrapz), isn't it?

is there any methodolgy or another function to numerically integrate the vector.

please help.

% this program is for calculating the NCG (non-circular gear) with no approximation
clear, clc

disp('This program generates NCG depending on the imported speed ratio form Excel file')
disp('and entered center distance')

Ratio=xlsread('speed_ratio.xlsx','A:A');
Ratio=Ratio';
Center_dis=89.305;
m=3;
disp('*')
disp('*')
disp('The solution depends on the main equation of NCG design and geanaration')
disp('*')
disp('*')

% the input angle is equal parts since the driving gear angular velocity is
% constant, so depending on the number of the speed ration elements, we
% can divide the one rotation (2*pi) to the same number of elements
Input_ang=0:0.0001*pi:2*pi;

% plot the speed ratio variations 
figure(1)
% plot the center
plot(0,0,'or','LineWidth',5)
hold on

plot(Input_ang,Ratio,'LineWidth',3)
title('The speed ratio variations of the NC-gearset')
xlabel('input angle')
ylabel('W_2/ W_1')
axis equal
grid on

% generating the first centrode with its evolute, base curve, tip and root curve. 
R_c1=(Center_dis)./(1+Ratio);
X_c1=R_c1.*cos(Input_ang);
Y_c1=R_c1.*sin(Input_ang);
       % plot it with clear zero point
       figure(2)
       % plot the center
       plot(0,0,'or','LineWidth',5)
         hold on
       plot(X_c1,Y_c1,'k')
        title('The first centrode')
        xlabel('X')
        ylabel('Y')
        axis equal
        grid on

% generating the second centrode with its evolute, base curve, tip and root curve.

% R_c2=Center_dis-R_c1;                    % first choice
R_c2=(Center_dis.*Ratio)./(1+Ratio);       % second choice


% finding the rotating angle of the second gear >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> here is my problem >>>>> the last value of this output angle must be equal to -pi
Ratio_2=1./Ratio;
out_ang_trap=cumtrapz(Input_ang,Ratio_2);
out_ang_trap=pi-out_ang_trap;

X_c2=R_c2.*cos(out_ang_trap);
Y_c2=R_c2.*sin(out_ang_trap);


       % plot it with clear zero point
       figure(3)
       % plot the center
       plot(0,0,'or','LineWidth',5)
         hold on
       plot(X_c2,Y_c2,'k')
        title('The second centrode')
        xlabel('X')
        ylabel('Y')
        axis equal
        grid on

Answer 1

[John D'Errico]

This is something covered in any basic class that talks about numerical methods. (And why are you using methods about which you seem to have not a clue? Time to do some reading!)

Trapezoidal rule is about the simplest integration rule you can use (bested in that respect perhaps only by methods such as rectangle rule, or arguably Monte Carlo, athough Monte Carlo falls into an entirely different class of methods, so a bad comparison.) As such, trapezoidal rule is not expected to be accurate. So there should be no surprise at all in this.

The classical way to decrease the error arising from trapezoidal rule (which is what cumtrapz is) is to decrease the step size, thus increase the number of points.

There are other methods you can use. Surely you can find at least one Simpson's rule tool on the file exchange. It will be more accurate, BUT you need to respect that no numerical integration method will be perfectly accurate in something like this! If you have not a clue as to what Simpson's rule is, then you really do need to go back and read your class notes from some long forgotten class.)

Or, you can use a higher order Newton-Cotes rule, but beware of too high an order, as I would bet that they too can introduce interesting artifacts. As well, higher order Newton-Cotes rules have an issue at the end point, since your data may not have the correct number of points to make up a complete panel. Even Simpson's rule suffers from this flaw, since it really wants an odd number of points to work correctly.

Other methods are things like fitting your curve with a spline, then integrating the spline. fnint will do the integration. (It is in the curve fitting toolbox, so you would need that. But there are simple enough ways to integrate a spline.)

I would suggest the simplest thing is a spline integral, then integrating the spline. But depending upon the shape of your integration kernel, a pchip integral may often be a better idea. At the same time, a pchip interpolant will be a lower order approximation, so it too will sometimes be problematic.

Sorry, but it seems to me you are trying to solve a problem without understanding the methods you are using. Would you recommend that someone play around with a table saw, who has never seen one, and has no idea how they work or what to do?

If this is a problem from a class, where you were told to use cumtrapz, then I would hope there would be discussion afterwards about the error that you see, and the source of that error.