Multiplication of large number with small number

I'm trying to compute a rather ugly integral using MATLAB. What I'm having problem with though is a part where I multiply a very big number (>10^300) with a very small number (<10^-300). MATLAB returns 'inf' for this even though it should be in the range of 0-0.0005. This is what I have

    besselFunction = @(u)besseli(qb,2*sqrt(lambda*(theta + mu)).*u);
    exponentFuncion = @(u)exp(-u.*(lambda + theta + mu));

where qb = 5, lambda = 12, theta = 10, mu = 3. And what I want to find is

    besselFunction(u)*exponentFunction(u)

for all real values of u. The problem is that whenever u>28 it will be evaluated as 'inf'. This is because the besselFunction(29) is so large even though exponentFunction(29) is extremely small. I've heared of, and tried, to use MATLAB function 'vpa' but it doesn't seem to work well when I want to use functions...

Any tips will be appreciated at this point!

[laurie]

sorry i feel stupid for even suggesting this, but can't you divide your very large number bessel(u>28) by, say, 10^150, do your multiplication, and then multiply 10^150 again after that ? or is bessel (u>28) evaluated as 'inf' ?

[Filip Trönnberg]

Was my first thought too but unfortunately bessel(u>28) is evaluated as inf.

Answer 1

[Walter Roberson]

vpa() is for use with symbolic expressions only.

If you have the symbolic engine,

    besselFunction = @(u)besseli(qb,2*sqrt(lambda*(theta + mu)).*sym(u));
    exponentFuncion = @(u)exp(-sym(u).*(lambda + theta + mu));

Then besselFunction(u)*exponentFunction(u) will return a symbolic value. You can double() it if you want the double-precision representation.

Note: going symbolic will only postpone the problem: the symbolic engine is limited to about a billion decimal places.

[John D'Errico]

Walter - The number of digits is not the issue here, but the dynamic range of a double. That is virtually unlimited for either vpa or hpf, but of course there are always limits.

[Walter Roberson]

The dynamic range for symbolic numbers only goes up to somewhere around 10 to 1 billion.

[Filip Trönnberg]

Great suggestion, thanks! It is correct that the symbolic numbers could pose a limitation but for my needs it was more than enough. What John wrote do also work but requires some extra thinking!

Answer 2

[John D'Errico]

The simple answer is always to use logs. The log of these numbers will be well behaved. Add the logs, then exponentiate at the end.

You could also use my HPF toolbox. It is now on the FEX, although I've not yet implemented bessel functions.

And, of course, the symbolic toolbox is an option with vpa.

[Filip Trönnberg]

Nice toolbox, will check it out. Would really appreciate of you could be more specific when you say I could try to add the logs and then exponentiate at the end. Regards

[Sargondjani]

he means: take logs, add them, take exponential. this is equivalent to multipying 2 numbers:

a x b = exp (log (a x b)) = exp(log a + log b)

[John D'Errico]

Exactly as Sargondjani says. A problem may be you need the log of this bessel function, if it is itself too large for the dynamic range of a double. For that you might need to use some approximations. Perhaps a series approximation would be adequate for the log of the bessel function. Or, one day I'll get bessel functions written for HPF, or maybe someone else will do that.