- Forum
- General C++ Programming
- Numerical analysis problem
Numerical analysis problem
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
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
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.
Approximating the pow(a-b, 11) using series and expansions was with no luck so far.
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).
If they are distances associated with subatomic particles then work in different length units (e.g. Angstrom units or smaller, rather than metres).
Last edited on
That is the Coulomb integral between two interacting 3s Slater atomic orbitals in 3D Cartesian space.
You can find more info about the integrand here http://www.theochem.ru.nl/~pwormer/Knowino/knowino.org/wiki/Slater_orbital.html
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.
If your limits are such that a is actually allowed to be equal to b then the integral won't exist.
Last edited on
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.
a = r * x_i
b = r * x_j
I calculated the limit of the integral as r -> 0 which is implemented separately.
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.
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.
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.
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.
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.
Last edited on
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.
Last edited on
Yeah not so trivial. I am trying to find another way to implement the code.
Thank you all,
M
Thank you all,
M
Topic archived. No new replies allowed.