Using the integral function for large upper limit

how do I know if integral(fun,x0,inf) is being accurate or not? Ive noticed that replacing inf by a very large number (given i know the fun decays sharply with x) gives a different and smaller answer? Why is this the case and what kind of magic is matlab doing when I give it an infinite limit?

Also, I notice a rescaling of my variables gives a different answer still. In SI units the lower bound x0 is of order 10^15. I have rescaled my units such that x0 is of order 1. Each gives a very different answer (yes i have accounted for converting back the units)...which can i trust?

[Torsten]

Play with the tolerances AbsTol and RelTol. Choose these parameters such that the value for the integral does not change much if you tighten them.

[will steel]

@Torsten This does seem to work for some numbers in the x0 = O(1) case - if i tighten tolerances sufficiently it reaches the inf limit value. However, for the SI unit case i run into this problem...

Warning: Minimum step size reached near x = 3.29e+15. There may be a singularity, or the tolerances may be too tight for this problem.

this is only for 'AbsTol',10^-6,'RelTol',10^-4 and there is certainly no singularity in the function (i have checked by plotting it).

Does integral just fail for large limits or something?

[Torsten]

I don't know how the error is estimated in "integral", but to request an absolute error of 1e-6 for an integral with a very high value (e.g. 10000 or something) is of course inadequate and not to accomplish. So scaling the function to be integrated is recommended.

[David Goodmanson]

Hi will,

could you supply an example integrand and limits of integration for the first paragraph in the question, and also for the SI case that does not work in the second paragraph? In the SI case it sounds, as Torsten mentioned, that 1e-6 is just too small a tolerance when the independent variable is as large as it its. As a hand-waving argument, 1e-6/3e15 is about 3e-22, which is far too small a relative number to be handled by double precision numbers which are good down to about a part in 10^16

[will steel]

@David Goodmanson@Torsten

Ah yes thanks i should have given a more lenient abstol. I have tried changing it to a large tolerance and though i lost the warning message i got the same result sadly. For more detail: this is how i am using integral:

fun = @(v) ((R/r)^2 * 2*pi*av0*v0^3 /(c^2))*(exp(-tv)./(v.*(exp(h*v/(k*T))-1)));

f = integral(fun, v0,inf,'AbsTol',10^-1)

The integral is over v and all the other symbols are physics constants (ill put them below in SI units). Worth noting is that the rescaling I used is seconds --> 10^-15 seconds (fs) and the integrand has no seconds unit - so the rescaling should increase the value of f by 10^15 (from dv scaling) I think?. However with the fs rescaling i get 6.65...e-11, and with the normal SI units i get 0.00531... which dont seem to share a simple relationship

v0 = 3.29*10^15; %this is the integral lower limit (3.26 if i rescale to fs)

c = 2.998*10^8; h = 6.626*10^-34; k = 1.381*10^-23;

av0 = 6.3*10^-22;

T = 20*10^3; nH = 10^6; L = 4*10^29; aB = 2.59*10^-19;

R = (15*L*h^3*c^2/(8*pi^6*k^4*T^4))^(1/2);

r = R;

tv = 0; %for checking integral

[David Goodmanson]

Hi will,

this all can be reduced to two specific situations with the integral function. Your integral falls off like 1 / ( x*exp(x) ), but x varies slowly enough that the effects show up with just the use of exp(-x).

[1] inf replaced by a large number

fun = @(x) exp(-x);
integral(fun,1,inf)
    ans = 0.3679        % 1/e, correct   
integral(fun,1,1e7)
    ans = 0.3679        % correct
integral(fun,1,1e8)
    ans = 1.8175e-22

I have wondered in the past why that is.

[2] Scaling the variable of integration. Should still give 1/e. Note that at the lower limit of integration, the argument of exp is always -1, so it's not like exp is immediately forced into overflow or underflow..

fun = @(x,A) A*exp(-A*x);

A = 1;    integral(@(x) fun(x,A),1/A,inf)
    ans = 0.3679        % correct
other values of A:
A = 1e-13
    ans = 0.3679        % correct
A = 1e-14
    Warning: Minimum step size reached near x = 1e+14, etc.
    ans = 1.0972e-07
A = 1e8
    ans = 0.3679        % correct
A = 1e9
    ans = 1.6570e-77

so a fair amount of scaling works. In your case the equivalent of A is h/(kT) = 2.3990e-15 which looks to be too small to work.

[will steel]

@David Goodmanson Ah okay i will use the scaled version and just note matlabs struggles for future reference. Many thanks

[Torsten]

@will steel You might also test MATLAB's "quad" for integration. "Old" does not always mean "Bad" :-)