CODEWITHSUNDEEP
matlabmathnumberslargenumber
Groot
Groot

Reputation: 14251

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'. I've heared, 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!

Upvotes: 1

Views: 1666

Answers (2)

craigim
craigim

Reputation: 3914

If you are doing an integration, consider using Gauss–Laguerre quadrature instead. The basic idea is that for equations of the form exp(-x)*f(x), the integral from 0 to inf can be approximated as sum(w(X).*f(X)) where the values of X are the zeros of a Laguerre polynomial and W(X) are specific weights (see the Wikipedia article). Sort of like a very advanced Simpson's rule. Since your equation already has an exp(-x) part, it is particularly suited.

To find the roots of the polynomial, there is a function on MATLAB Central called LaguerrePoly, and from there it is pretty straightforward to compute the weights.

Upvotes: 0

duffymo
duffymo

Reputation: 308763

I'd use logarithms.

Let x = Bessel function of u and y = x*exp(-u) (simpler than your equation, but similar).

Since log(v*w) = log(v) + log(w), then log(y) = log(x) + log(exp(-u))

This simplifies to 

log(y) = log(x) - u

This will be better behaved numerically.

The other key will be to not evaluate that Bessel function that turns into a large number and passing it to a math function to get the log. Better to write your own that returns the logarithm of the Bessel function directly. Look at a reference like Abramowitz and Stegun to try and find one.

Upvotes: 5

Related Questions