AbsTol and RelTol in vpaintegral

Hi

I am integrating some complex equations using vpaintegral in matlab. it is triple integration and the equations are too complex, so the calculation time is too high

when the reltol and abstol is set to 1e-1, it takes about 3-4 hours to compute the result, but I need more accuracy.

when I set reltol and abstol to 1e-2 it takes about 12 hours!! and yet I need more accuracy!!

The questions is that is it necessary to give both reltol and abstol values to get accurate result?? or one of them are enough in my work? and does it make any difference in calculation time to set only one of them (or it doesn't affect)? if so, which one should I use?

(the value of my final result is between 990 to 1001)

thank you

Answer 1

[Walter Roberson]

Internally, vpaintegral() always works with tolerances. It has default values that it uses. It is not the case that providing a tolerance as an option results in the code proceeding differently -- no

if relative_tolerance_was_specified
    do some extra work
end

So you should not worry about only providing one instead of the other to have it do less work: it is always going to be doing the work, just with default values if you did not specify the option.

Typically reltol is more important in figuring out when to stop. Typically abstol is there as a fall-back for places where the relative values get small. You know that your final result is roughly 1000, but there could still be ranges of your integral where the values are (say) near 1e-15, and over that stretch vpaintegral() might be doing a lot of fine work to try to be more precise about the small numbers, but you in your knowledge of the overall range of results might say that it is pointless to be precise down below 1e-15 about a value that is effectively going to contribute nothing to the overall result, so you might want to put in an absolute tolerance of (say) 1e-5.

Absolute tolerance of 1e-5 does not mean that it will persue the integration until it is certain that the result over the entire integration is within 1e-5 of the "correct" value: instead it means that values less than 1e-5 would more or less filtered out.

There are some cases where the absolute values might get quite small, and yet the contribution could be very important. For example the integral of 1/x, as x goes to infinity, is infinite, ln(infinity), and yet the contributions up around realmax are only about 5e-309 each, with the total contribution by realmax only being about 709 -- the contributions beyond realmax are all negligibly small from the perspective of IEEE 754 double precision, and yet in theory they contribute an infinite amount out to infinity.

Therefore, it is not always appropriate to use an absolute tolerance: you have to know enough about the function to be able to establish that the integral of the missed values will not matter to your accuracy.

In the meantime, I recommand that you use matlabFunction() to turn the vpaintegral() with symbolic upper bound into an anonymous function that calls integral()

I was about to suggest you consider using the matlabFunction() options 'File' and 'optimize', but then I remembered that I have open cases against bugs in the optimization function that can result in the optimized function being quite wrong. It might still be worth using the 'File' option, but if you do, then explicitly set 'optimize' to false.

[Pooyan Jafari]

Hello again,

Thank you for your answer

I am not familiar with matlabFunction and I don't know how to use it. I also read matlab help about this function but I didn't understand how to apply it to my code. would you please take a berif look at my code and tell me where to use matlabfunction()?

should I use vpaintegral() or integral() when I use matlabfunction? and what is File option?

Thank you so much for your kind help

clc,clear,close all;
opengl('save', 'software')
syms beta theta r ;

B = 10;
H = 10;
C = 10;
L1 = 4;

Ri = B/2 + L1;
Rf = C + B/2 + L1;
Rbeta = Ri + (2*beta*(Rf-Ri)/pi);
Rbeta0 = Ri;
Rbetapi2 = Rf;
theta1 = atan((H/2)/(B/2+L1));
theta2 = pi - atan((H/2)/(B/2-L1));
Rt0 = piecewise((0<=theta) & (theta<theta1),(B/2+L1)/cos(theta),(theta1<=theta) & (theta<theta2),(H/2)/sin(theta),(theta2<=theta) & (theta<=pi),(B/2-L1)/cos(pi-theta));

Rmaxtb = (Rt0*Rbeta)/Ri;
Rmaxtb0 = (Rt0*Rbeta0)/Ri;
Rmaxtbpi2 = (Rt0*Rbetapi2)/Ri;

Vmbeta = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B)+(4*beta*C))^2));
Vmbeta0 = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B))^2));
Vmbetapi2 = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B)+(4*(pi/2)*C))^2));

Vbeta = Vmbeta * (1 - ((r^2) / (Rmaxtb^2)));
Vbeta0 = Vmbeta0 * (1 - ((r^2) / (Rmaxtb0^2)));
Vbetapi2 = Vmbetapi2 * (1 - ((r^2) / (Rmaxtbpi2^2)));

a_coeff = 1 / ( (Rmaxtb^3) * ((int(Rmaxtb^2,theta,0,pi))^2) * r * (r*cos(theta)-Rf));
b_coeff = (((r^2) * Rmaxtb) - (Rmaxtb^3)) * (int(Rmaxtb * diff(Rmaxtb,beta),theta,0,pi));
c_coeff = ( int(Rmaxtb^2,theta,0,pi) ) * ( diff(Rmaxtb,beta) ) * ((r^2) + (Rmaxtb^2));
d_coeff = (H*(r-Rmaxtb)) * (B*(r+Rmaxtb)) * (b_coeff + c_coeff);
Vr = a_coeff * d_coeff;

epsrr = diff(Vr,r);
epsrtheta = (1/(2*r)) * (diff(Vr,theta));
epsrbeta = (1/(2*(Rf-(r*cos(theta))))) * ( (Vbeta*cos(theta)) + (diff(Vr,beta)) + ( (Rf-(r*cos(theta)))*(diff(Vbeta,r)) ) );
epsthetatheta = Vr/r;
epsthetabeta = (1/(2*r*(Rf-(r*cos(theta))))) * ( (-1*Vbeta*r*sin(theta)) + ( (Rf-(r*cos(theta)))*(diff(Vbeta,theta)) ) );
epsbetabeta = ( 1/(Rf-(r*cos(theta))) ) * ( (-1*Vr*cos(theta)) + (diff(Vbeta,beta)) );

g_coeff = (epsrtheta^2) + (epsthetabeta^2) + (epsrbeta^2);
h_coeff = (epsrr * epsthetatheta) + (epsthetatheta * epsbetabeta) + (epsbetabeta * epsrr);
p = g_coeff - h_coeff;

i_coeff = (epsrr * epsthetatheta * epsbetabeta) + (2*epsrtheta*epsthetabeta*epsrbeta);
j_coeff = (epsbetabeta * (epsrtheta^2)) + (epsrr * (epsthetabeta^2)) + (epsthetatheta * (epsrbeta^2));
q = i_coeff - j_coeff;

epsmax = 2 * sqrt(p/3) * cos((1/3) * acos( (3*sqrt(3)*abs(q)) / (2*sqrt(p^3)) ));

Dc = 2 * vpaintegral(vpaintegral(vpaintegral( (2*epsmax) *r*(Rf-(r*cos(theta))), r, [0 Rmaxtb], 'AbsTol',1e-1, 'RelTol',1e-1), theta, [0 pi], 'AbsTol',1e-1, 'RelTol',1e-1),beta ,[0 pi/2], 'AbsTol',1e-1, 'RelTol',1e-1)