Numerical analysis problem

[mmgh]

Hi everyone!

I have an issue regarding numerical analysis. I have written a c++ code to solve an equation where there is a term like this in the denominator:

1/pow(a-b, 11)

the problem is that if 0 < a-b << 1, then pow(a-b, 11) becomes zero with the double precision and thus the division is basically a NAN.

I implemented the code with the CLN library to have an arbitrary high precision and it works with 60 digits. However, my code has become extremely slow.

I was wondering if there was any other solution better than the CLN library.

Thank you in advance for you help.

Cheers,
M

[lastchance]

What problem are you trying to solve?

There is probably a better algorithm.

[mmgh]

I am trying to analytically solve an integral for two interacting particles. I solved the integral with Mathematica and then implemented it in C++.

Approximating the pow(a-b, 11) using series and expansions was with no luck so far.

[lastchance]

What integral? (Give integrand and limits).

If they are distances associated with subatomic particles then work in different length units (e.g. Angstrom units or smaller, rather than metres).

[mmgh]

That is the Coulomb integral between two interacting 3s Slater atomic orbitals in 3D Cartesian space.

[mmgh]

You can find more info about the integrand here http://www.theochem.ru.nl/~pwormer/Knowino/knowino.org/wiki/Slater_orbital.html

[lastchance]

If a and b are lengths then work in different length units (e.g. Angstroms for typical atomic sizes). Then they won't be so numerically small. Otherwise, non-dimensionalise your problem with respect to, say, the distance between particle centres; again, that would keep numerical factors to more sensible sizes.

If your limits are such that a is actually allowed to be equal to b then the integral won't exist.

[mmgh]

Lets say we have two 3s Slater orbitals i and j and r is their distances and x is the exponent (width) of the orbital:

a = r * x_i
b = r * x_j

I calculated the limit of the integral as r -> 0 which is implemented separately.

[lastchance]

a = r * x_i
b = r * x_j

OK, that will make x_i and x_j into non-dimensional variables with sensible sizes. Do the integral with x_i and x_j rather than a and b.

[mmgh]

r is in nm and x_i and x_j are in 1/nm so a and b are dimensionless.

[jonnin]

what computer are you on? Some computers have up to 128 bit double values available in the hardware. I don't know how many digits that is off the top of my head.

it comes down to 3 answers, always, for this type of question..
- do it in hardware if you can
- change the code to avoid the big/tiny number problem if you can somehow
- else do it with emulation, which is slow and your last resort.

[lastchance]

You still haven't stated what integral or expression you are trying to evaluate. It isn't in the link posted, and it's difficult to offer any advice without it.

For a simple electrostatic interaction between between two overlapping spherical orbitals it is hard to see where your power of -11 comes from, and you may also get near cancellation from two opposing terms. We need to see the full expression.

[Ganado]

In addition to what lastchance said, I would be very cautious of doing anything where you have (unit)*(1/unit) where unit is a very big or small number -- the floating-point rounding error can grow to be very inaccurate. I would suggest finding a way to refactor this or use different units.

[mmgh]

Yeah not so trivial. I am trying to find another way to implement the code.

Thank you all,
M