1
$\begingroup$

I have a complicated function called f2[re]

f2[re_] := 
  1/(2 cep^(3/2) (cep - 2 cp) cp^(3/2) re)
    E^(-cp inf^2 - cep inf re - ce re^2 - (cep re \[Alpha]2)/cp - 
    2 (inf + re) \[Alpha]2)
    norm^2 \[Pi] (2 cep^(3/2) Sqrt[cp] E^((cep re \[Alpha]2)/cp) - 
     4 Sqrt[cep] cp^(3/2) E^((cep re \[Alpha]2)/cp) - 
     2 cep^(3/2) Sqrt[cp] E^(
      re (2 cep inf + 4 \[Alpha]2 + (cep \[Alpha]2)/cp)) + 
     4 Sqrt[cep] cp^(3/2) E^(
      re (2 cep inf + 4 \[Alpha]2 + (cep \[Alpha]2)/cp)) + 
     4 cp^(3/2) E^(
      2 (inf + re) \[Alpha]2 + (2 \[Alpha]2^2)/cep + (
       cep (cp (inf + re)^2 + 2 re \[Alpha]2))/(2 cp)) Sqrt[
      2 \[Pi]] \[Alpha]2 Erf[(cep (inf - re) + 2 \[Alpha]2)/(
       Sqrt[2] Sqrt[cep])] + 
     cep^(3/2) E^(
      cp inf^2 + cep inf re + 2 inf \[Alpha]2 + 
       4 re \[Alpha]2 + ((cep^2 re^2)/4 + \[Alpha]2^2)/cp)
       Sqrt[\[Pi]] (cep re - 2 (cp re + \[Alpha]2)) Erf[(
       2 cp inf - cep re + 2 \[Alpha]2)/(2 Sqrt[cp])] + 
     cep^(5/2) E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]]
       re Erf[(2 cp inf + cep re + 2 \[Alpha]2)/(2 Sqrt[cp])] - 
     2 cep^(3/2) cp E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]]
       re Erf[(2 cp inf + cep re + 2 \[Alpha]2)/(2 Sqrt[cp])] + 
     2 cep^(3/2) E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]] \[Alpha]2 Erf[(2 cp inf + cep re + 2 \[Alpha]2)/(
       2 Sqrt[cp])] - 
     cep^(5/2) E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]]
       re Erf[(cep re - 2 cp re + 2 \[Alpha]2)/(2 Sqrt[cp])] + 
     2 cep^(3/2) cp E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]]
       re Erf[(cep re - 2 cp re + 2 \[Alpha]2)/(2 Sqrt[cp])] - 
     2 cep^(3/2) E^(
      cp inf^2 + cep inf re + (cep^2 re^2)/(4 cp) + 
       2 inf \[Alpha]2 + (2 cep re \[Alpha]2)/cp + \[Alpha]2^2/cp)
       Sqrt[\[Pi]] \[Alpha]2 Erf[(cep re - 2 cp re + 2 \[Alpha]2)/(
       2 Sqrt[cp])] - 
     4 cp^(3/2) E^(
      2 (inf + re) \[Alpha]2 + (2 \[Alpha]2^2)/cep + (
       cep (cp (inf + re)^2 + 2 re \[Alpha]2))/(2 cp)) Sqrt[
      2 \[Pi]] \[Alpha]2 Erf[(cep (inf + re) + 2 \[Alpha]2)/(
       Sqrt[2] Sqrt[cep])] + 
     cep^(5/2) E^(
      cp inf^2 + cep inf re + 2 inf \[Alpha]2 + 
       4 re \[Alpha]2 + ((cep^2 re^2)/4 + \[Alpha]2^2)/cp) Sqrt[\[Pi]]
       re Erf[(cep re - 2 (cp re + \[Alpha]2))/(2 Sqrt[cp])] - 
     2 cep^(3/2) cp E^(
      cp inf^2 + cep inf re + 2 inf \[Alpha]2 + 
       4 re \[Alpha]2 + ((cep^2 re^2)/4 + \[Alpha]2^2)/cp) Sqrt[\[Pi]]
       re Erf[(cep re - 2 (cp re + \[Alpha]2))/(2 Sqrt[cp])] - 
     2 cep^(3/2) E^(
      cp inf^2 + cep inf re + 2 inf \[Alpha]2 + 
       4 re \[Alpha]2 + ((cep^2 re^2)/4 + \[Alpha]2^2)/cp)
       Sqrt[\[Pi]] \[Alpha]2 Erf[(cep re - 2 (cp re + \[Alpha]2))/(
       2 Sqrt[cp])]);

The values of function f2 must be equal to function f1 in different re. function f1 includes a numerical integration as follows

f1[re_?NumericQ] := (2*\[Pi]*norm^2)/
   re*(Exp[-ce re^2]) (1/cep^(3/2)) NIntegrate[
    rp *Exp[-cp rp^2]*(E^(-cep re rp) (Sqrt[
          cep] (-E^(-2 (re + rp) \[Alpha]2) + E^(
            2 cep re rp - 2  \[Alpha]2  Abs[re - rp])) - 
         E^(1/2 cep (re + rp)^2 + (2 \[Alpha]2^2)/cep) Sqrt[
          2 \[Pi]] \[Alpha]2 Erf[(cep (re + rp) + 2 \[Alpha]2)/(
           Sqrt[2] Sqrt[cep])] + 
         E^(1/2 cep (re + rp)^2 + (2 \[Alpha]2^2)/cep) Sqrt[
          2 \[Pi]] \[Alpha]2 Erf[(2 \[Alpha]2 + cep Abs[re - rp])/(
           Sqrt[2] Sqrt[cep])])), {rp, 0, inf}];

So I specify the needed constants as follows

ce = 0.0006533997022026266`;
cp = 2.8846533997022026`;
cep = -0.026693200595594747`;
inf = 15;
norm = 0.5293704593281444`;
{\[Alpha]2, \[Beta]2} = {0.9948987042561277`, 0.0002930460049474671`};

when I try

f1[0.1]
f2[0.1]

I get

0.0968872575537 + 0. I
0.0968872575541 + 0. I

which is OK as I expect. But as soon as I set inf=16 it returns

General::munfl: Exp[-770.463385382] is too small to represent as a normalized machine number; precision may be lost.
0. + 0. I

I googled this error and it said it relates to MachinePrecision, It's obvious if inf>15 some parts of f2[re] become so small so that Mathematica makes them zero, so the whole function becomes zero. However I don't any idea to address this problem and increase the precision of calculations. Any idea?

Note: inf is finite number (its maximum is 20) which I used as upper limit of integration to obtain f2[re]

$\endgroup$
1
  • 1
    $\begingroup$ You will need to work in higher precision. This can be done by making the constants into bignums (say, precision of 30), setting NIntegrate to use WorkingPrecision->25 (or anything modestly lower than the precision of the constants), and using a bignum instead of .1 in the inputs. $\endgroup$ Apr 27, 2021 at 15:36

2 Answers 2

1
$\begingroup$

Your function f1 is not relevant to your issue, and your definition of f2 could be simplified a little using Mathematica's Simplify. You need to ensure none of the numeric values are machine precision numbers as I do below.

ce=0.0006533997022026266`25;
cp=2.8846533997022026`25;
cep=-0.026693200595594747`25;
inf=16;
norm=0.5293704593281444`25;
{\[Alpha]2,\[Beta]2}={0.9948987042561277`25,0.0002930460049474671`25};
f2[1/10]
(*  0.09688725755409674223410719883437690442`20.5674697906684  *)

In the above, Mathematica's tracking of precision predicts the result is correct to about 20.6 decimal places. The advantage of Machine Precision numerics is tat it runs faster than arbitrary precision numerics. You can learn about arbitrary precision numbers in Mathematica here and here.

$\endgroup$
1
  • $\begingroup$ Thanks a lot, I wonder how you did this?!! I read your references, but can you tell me what does25 at the end of each constant value do? and why f2[0.1] does not the correct answer but f2[1/10] does? $\endgroup$
    – Wisdom
    Apr 27, 2021 at 15:47
2
$\begingroup$

This is too long for a comment. On my computer Mathematica treats 3.141592653589793 as a machine precision number. In Mathematica computing with machine precision numbers works like most computing software. Hence, I get the floating point approximation of 0 from the following.

1+Cos[3.141592653589793]

However, 3.14159265358979300000000000 is not a machine precision number, and Mathematica is much more careful (and slower) when computing with a number such as this, so the result is not 0.

1+Cos[3.14159265358979300000000000]
(* 2.843221614*10^-32 *)

Instead of typing 3.14159265358979300000000000 we can type 3.141592653589793`25 which means this number should be considered correct to 25 digits to the right of the decimal point. Now consider the following:

0.1 + 2.123456789123456789123456789

In the above the machine precision number 0.1 is added to an arbitrary precision number. Except when a result is exact, computation with machine precision numbers returns machine precision numbers. By default Mathematica only displays the first six digits of machine numbers. You can copy the output of the last result into a new cell and you will see the result is the machine number 2.223456789123457`. When you need to compute with arbitrary precision numbers, you normally need to avoid using machine precision numbers for the reason just illustrated. Now look at the the following:

1/10+2.123456789123456789123456789

That involves no machine numbers, so the result will have many digits of precision. I used 1/10 instead of 0.1 in my first answer because I needed to avoid results with machine precision.

We also have the following:

2.345 (* A machine precision number. *)
2.345`9 (* Not a machine precision number *)

For the second example above, Mathematica will treat it as correct to 9 digits to the right of the decimal point. The documentation points out the following:

123.45678901234567890 (* A machine precision number on some computers. *)
123.45678901234567890` (* A machine precision number on all computers. *) 

We also have:

123.45678901234567891234567891234` (* A machine precision number *) 

However, with the last example, digits too far to the right for a machine precision number are be ignored.

It also helps to understand the output of the following:

 1+Cos[3.14159265358979323846264338328]

The output is 0.*10^-59. This is approximately zero, but is not a machine precision number. It is between -10^-59 and 10^-59.

$\endgroup$
1
  • $\begingroup$ Thanks a lot for your detailed answer. $\endgroup$
    – Wisdom
    Apr 28, 2021 at 3:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.