# General:munfl error problem

I have the following function, that is the result of The integration between x and some upper limit of a positive function, thus should be positive

$$f(x,vE) = \frac{e^{-1140.09 x} \left(e^{67945.5 \text{vE}^2} \theta (\text{1.859514170346336\grave{ }*{}^{\wedge}-9}-x) \left(e^{1140.09 x} \left(-\text{3.3591048468065713\grave{ }*{}^{\wedge}-37} \text{vE}^8+\text{1.2371349122599874\grave{ }*{}^{\wedge}-16} x^8+\text{4.6259292692714846\grave{ }*{}^{\wedge}-18} x^7-\text{5.720218174055417\grave{ }*{}^{\wedge}-21} x^6-\text{3.830505707541381\grave{ }*{}^{\wedge}-19} x^5-\text{6.875632678803749\grave{ }*{}^{\wedge}-22} x^4-\text{2.592004365421134\grave{ }*{}^{\wedge}-25} x^3+\text{1.1383704954727588\grave{ }*{}^{\wedge}-27} x^2-\text{1.0946897704176435\grave{ }*{}^{\wedge}-14} x+\text{2.0355911402247853\grave{ }*{}^{\wedge}-23}\right)+\text{1.1896770509204967\grave{ }*{}^{\wedge}-8} x^5-\text{8.30796052725145\grave{ }*{}^{\wedge}-9} x^4+\text{1.4721183020268512\grave{ }*{}^{\wedge}-9} x^3-\text{7.569314750088174\grave{ }*{}^{\wedge}-12} x^2+\text{1.0946949062266079\grave{ }*{}^{\wedge}-14} x-\text{2.0355980730229598\grave{ }*{}^{\wedge}-23}\right)+e^{67945.5 \text{vE}^2} \left(-\text{1.1896770509204967\grave{ }*{}^{\wedge}-8} x^5+\text{8.30796052725145\grave{ }*{}^{\wedge}-9} x^4-\text{1.4721183020268512\grave{ }*{}^{\wedge}-9} x^3+\text{7.569314750088174\grave{ }*{}^{\wedge}-12} x^2-\text{1.0946949062266079\grave{ }*{}^{\wedge}-14} x+\text{2.0355980730229598\grave{ }*{}^{\wedge}-23}\right)+e^{67945.5 \text{vE}^2+1140.09 x} \left(\text{vE}^{16} (1. x-\text{1.8595141703463354\grave{ }*{}^{\wedge}-9})+\text{vE}^{14} (\text{2.0083576467270298\grave{ }*{}^{\wedge}-14}-0.0000108004 x)+\text{vE}^{12} (\text{2.663670413232032\grave{ }*{}^{\wedge}-6} x-\text{4.953132878537244\grave{ }*{}^{\wedge}-15})+\text{vE}^{10} (\text{5.166994397269747\grave{ }*{}^{\wedge}-10} x-\text{9.608099299823222\grave{ }*{}^{\wedge}-19})+\text{vE}^8 (\text{2.571855800860126\grave{ }*{}^{\wedge}-21}-\text{1.383079431108135\grave{ }*{}^{\wedge}-12} x)+\text{vE}^6 (\text{1.7552515711539956\grave{ }*{}^{\wedge}-25}-\text{9.439301937812491\grave{ }*{}^{\wedge}-17} x)+\text{vE}^4 (\text{1.4313862961095703\grave{ }*{}^{\wedge}-21} x-\text{2.6616831008553027\grave{ }*{}^{\wedge}-30})+\text{vE}^2 (\text{2.9153325110971034\grave{ }*{}^{\wedge}-34}-\text{1.5677925759254207\grave{ }*{}^{\wedge}-25} x)-\text{1.2371349122599874\grave{ }*{}^{\wedge}-16} x^8-\text{4.6259292692714846\grave{ }*{}^{\wedge}-18} x^7+\text{5.720218174055417\grave{ }*{}^{\wedge}-21} x^6+\text{3.830505707541381\grave{ }*{}^{\wedge}-19} x^5+\text{6.875632678803749\grave{ }*{}^{\wedge}-22} x^4+\text{2.592004365421134\grave{ }*{}^{\wedge}-25} x^3-\text{1.1383704954727588\grave{ }*{}^{\wedge}-27} x^2+\text{1.2253923290809874\grave{ }*{}^{\wedge}-30} x-\text{1.6264357714319357\grave{ }*{}^{\wedge}-39}\right)+e^{33972.8 \text{vE}^2+1140.09 x} \left(\text{vE}^8 (0.00937982 x-\text{1.744191593358305\grave{ }*{}^{\wedge}-11})+\text{vE}^6 (\text{4.087595291558347\grave{ }*{}^{\wedge}-13}-0.000219821 x)+\text{vE}^4 (\text{1.3071488339510708\grave{ }*{}^{\wedge}-6} x-\text{2.430661779483706\grave{ }*{}^{\wedge}-15})+\text{vE}^2 (\text{4.1941725957388944\grave{ }*{}^{\wedge}-19}-\text{2.255520642231904\grave{ }*{}^{\wedge}-10} x)-\text{6.639200573944934\grave{ }*{}^{\wedge}-15} x+\text{1.2345687547022132\grave{ }*{}^{\wedge}-23}\right)+(\text{1.7586098278121368\grave{ }*{}^{\wedge}-14} x-\text{3.270159894926998\grave{ }*{}^{\wedge}-23}) e^{33972.8 \text{vE}^2+1140.09 x}\right) \theta \left(29.7983 \text{vE}^2-x+\text{1.859514170346336\grave{ }*{}^{\wedge}-9}\right)}{\left(e^{33972.8 \text{vE}^2} \left(0.00937982 \text{vE}^8-0.000219821 \text{vE}^6+\text{1.3071488339510708\grave{ }*{}^{\wedge}-6} \text{vE}^4-\text{2.2555206422319046\grave{ }*{}^{\wedge}-10} \text{vE}^2-\text{6.639200573944936\grave{ }*{}^{\wedge}-15}\right)+\text{1.758609827812137\grave{ }*{}^{\wedge}-14} e^{33972.8 \text{vE}^2}+e^{67945.5 \text{vE}^2} \left(1. \text{vE}^{16}-0.0000108004 \text{vE}^{14}+\text{2.6636704132320325\grave{ }*{}^{\wedge}-6} \text{vE}^{12}+\text{5.166994397269748\grave{ }*{}^{\wedge}-10} \text{vE}^{10}-\text{1.3830794311081354\grave{ }*{}^{\wedge}-12} \text{vE}^8-\text{9.439301937812495\grave{ }*{}^{\wedge}-17} \text{vE}^6+\text{1.4313862961095709\grave{ }*{}^{\wedge}-21} \text{vE}^4-\text{1.567792575925421\grave{ }*{}^{\wedge}-25} \text{vE}^2-\text{1.0946897704176435\grave{ }*{}^{\wedge}-14}\right)\right) (x-\text{1.859514170346336\grave{ }*{}^{\wedge}-9})}$$

I am trying to get values of this function for vE up to 0.141421, however I start getting problems at vE=0.10. This is probably due to underflow errors. I have tried wrapping the function inside a SetPrecision, as well as using

$$f[SetPrecision[x,$$MachienPrecision],SetPrecision[vE,$$MachienPrecision]]$$

with no luck. I am afraid that the problem is the precision of the numeric coefficients of the function, that need to be set to arbitrary precision to fix the underflow.

Any suggestion on how to fix this problem without having to rewrite the function (its generated automatically by my program, chaging the part that generates it would be cubersome)

EDIT: Inputform expression

((E^(33972.7745398073*vE^2 + 1140.0906489396789*x)*(1.2345687547022132*^-23 +
vE^6*(4.087595291558347*^-13 - 0.00021982060458281043*x) +
vE^2*(4.1941725957388944*^-19 - 2.255520642231904*^-10*x) +
vE^4*(-2.430661779483706*^-15 + 1.3071488339510708*^-6*x) + vE^8*(-1.744191593358305*^-11 + 0.00937982415607754*x) -
6.639200573944934*^-15*x) + E^(33972.774539807324*vE^2 + 1140.0906489396789*x)*
(-3.270159894926998*^-23 + 1.7586098278121368*^-14*x) + E^(67945.54907961463*vE^2)*
(2.0355980730229598*^-23 - 1.0946949062266079*^-14*x + 7.569314750088174*^-12*x^2 - 1.4721183020268512*^-9*x^3 +
8.30796052725145*^-9*x^4 - 1.1896770509204967*^-8*x^5) + E^(67945.54907961463*vE^2 + 1140.0906489396789*x)*
(-1.6264357714319357*^-39 + vE^14*(2.0083576467270298*^-14 - 0.000010800442818636717*x) +
vE^8*(2.571855800860126*^-21 - 1.383079431108135*^-12*x) +
vE^6*(1.7552515711539956*^-25 - 9.439301937812491*^-17*x) +
vE^2*(2.9153325110971034*^-34 - 1.5677925759254207*^-25*x) +
vE^4*(-2.6616831008553027*^-30 + 1.4313862961095703*^-21*x) +
vE^10*(-9.608099299823222*^-19 + 5.166994397269747*^-10*x) +
vE^12*(-4.953132878537244*^-15 + 2.663670413232032*^-6*x) + vE^16*(-1.8595141703463354*^-9 + 0.9999999999999998*x) +
1.2253923290809874*^-30*x - 1.1383704954727588*^-27*x^2 + 2.592004365421134*^-25*x^3 + 6.875632678803749*^-22*x^4 +
3.830505707541381*^-19*x^5 + 5.720218174055417*^-21*x^6 - 4.6259292692714846*^-18*x^7 -
1.2371349122599874*^-16*x^8) + E^(67945.54907961463*vE^2)*(-2.0355980730229598*^-23 + 1.0946949062266079*^-14*x -
7.569314750088174*^-12*x^2 + 1.4721183020268512*^-9*x^3 - 8.30796052725145*^-9*x^4 + 1.1896770509204967*^-8*x^5 +
E^(1140.0906489396789*x)*(2.0355911402247853*^-23 - 3.3591048468065713*^-37*vE^8 - 1.0946897704176435*^-14*x +
1.1383704954727588*^-27*x^2 - 2.592004365421134*^-25*x^3 - 6.875632678803749*^-22*x^4 -
3.830505707541381*^-19*x^5 - 5.720218174055417*^-21*x^6 + 4.6259292692714846*^-18*x^7 +
1.2371349122599874*^-16*x^8))*HeavisideTheta[1.859514170346336*^-9 - x])*
HeavisideTheta[1.859514170346336*^-9 + 29.79830996018*vE^2 - x])/
(E^(1140.0906489396789*x)*(1.758609827812137*^-14*E^(33972.774539807324*vE^2) +
E^(33972.7745398073*vE^2)*(-6.639200573944936*^-15 - 2.2555206422319046*^-10*vE^2 + 1.3071488339510708*^-6*vE^4 -
0.00021982060458281048*vE^6 + 0.009379824156077542*vE^8) + E^(67945.54907961463*vE^2)*
(-1.0946897704176435*^-14 - 1.567792575925421*^-25*vE^2 + 1.4313862961095709*^-21*vE^4 -
9.439301937812495*^-17*vE^6 - 1.3830794311081354*^-12*vE^8 + 5.166994397269748*^-10*vE^10 +
2.6636704132320325*^-6*vE^12 - 0.000010800442818636719*vE^14 + 1.*vE^16))*(-1.859514170346336*^-9 + x))


Sidenote: the function should indicate a probability, and is the integral of a probability density defined in [0,x_max], therefore it should hold $$f(0,vE)=1$$ for any $$vE$$

• Please provide copy & paste-able code (InputForm) so that we have something to work with. Commented Aug 3, 2021 at 23:36
• you are right, done! Commented Aug 4, 2021 at 7:34

You might try to use exact numbers like

x = 1/100; Table[
Rationalize[((E^(33972.7745398073*vE^2 + 1140.0906489396789*x)*
(1.2345687547022132*^-23 +
vE^6*(4.087595291558347*^-13 - 0.00021982060458281043*
x) +
vE^2*(4.1941725957388944*^-19 -
2.255520642231904*^-10*x) +

vE^4*(-2.430661779483706*^-15 +
1.3071488339510708*^-6*x) +

vE^8*(-1.744191593358305*^-11 + 0.00937982415607754*x) -
6.639200573944934*^-15*x) +
E^(33972.774539807324*vE^2 + 1140.0906489396789*x)*
(-3.270159894926998*^-23 +
1.7586098278121368*^-14*x) +

E^(67945.54907961463*vE^2)*(2.0355980730229598*^-23 -
1.0946949062266079*^-14*x +
7.569314750088174*^-12*x^2 -
1.4721183020268512*^-9*x^3 +
8.30796052725145*^-9*x^4 -
1.1896770509204967*^-8*x^5) +

E^(67945.54907961463*vE^2 +
1140.0906489396789*x)*(-1.6264357714319357*^-39 +

vE^14*(2.0083576467270298*^-14 -
0.000010800442818636717*x) +

vE^8*(2.571855800860126*^-21 - 1.383079431108135*^-12*x) +

vE^6*(1.7552515711539956*^-25 -
9.439301937812491*^-17*x) +

vE^2*(2.9153325110971034*^-34 - 1.5677925759254207*^-25*x) +

vE^4*(-2.6616831008553027*^-30 +
1.4313862961095703*^-21*x) +

vE^10*(-9.608099299823222*^-19 + 5.166994397269747*^-10*x) +

vE^12*(-4.953132878537244*^-15 +
2.663670413232032*^-6*x) +

vE^16*(-1.8595141703463354*^-9 + 0.9999999999999998*x) +
1.2253923290809874*^-30*x -
1.1383704954727588*^-27*x^2 +
2.592004365421134*^-25*x^3 +
6.875632678803749*^-22*x^4 +
3.830505707541381*^-19*x^5 +
5.720218174055417*^-21*x^6 -
4.6259292692714846*^-18*x^7 -
1.2371349122599874*^-16*x^8) +

E^(67945.54907961463*vE^2)*(-2.0355980730229598*^-23 +
1.0946949062266079*^-14*
x - 7.569314750088174*^-12*x^2 +
1.4721183020268512*^-9*x^3 -
8.30796052725145*^-9*x^4 +
1.1896770509204967*^-8*x^5 +

E^(1140.0906489396789*x)*(2.0355911402247853*^-23 -
3.3591048468065713*^-37*
vE^8 - 1.0946897704176435*^-14*x +
1.1383704954727588*^-27*x^2 -
2.592004365421134*^-25*x^3 -
6.875632678803749*^-22*x^4 -
3.830505707541381*^-19*x^5 -
5.720218174055417*^-21*x^6 +
4.6259292692714846*^-18*x^7 +
1.2371349122599874*^-16*x^8))*
HeavisideTheta[1.859514170346336*^-9 - x])*
HeavisideTheta[

1.859514170346336*^-9 + 29.79830996018*vE^2 -
x])/(E^(1140.0906489396789*x)*
(1.758609827812137*^-14*E^(33972.774539807324*vE^2) +
E^(33972.7745398073*vE^2)*
(-6.639200573944936*^-15 -
2.2555206422319046*^-10*vE^2 +
1.3071488339510708*^-6*vE^4 -
0.00021982060458281048*vE^6 +
0.009379824156077542*vE^8) +
E^(67945.54907961463*vE^2)*
(-1.0946897704176435*^-14 -
1.567792575925421*^-25*vE^2 +
1.4313862961095709*^-21*vE^4 -
9.439301937812495*^-17*vE^6 -
1.3830794311081354*^-12*vE^8 +
5.166994397269748*^-10*vE^10 +
2.6636704132320325*^-6*vE^12 -
0.000010800442818636719*vE^14 + 1.*vE^16))*
(-1.859514170346336*^-9 + x)), 10^(-100)], {vE, 1/100,
14/100, 1/100}]
(* {0, 430648/6879181385, 698481/9203324930, 1805038/23785891251,
3126065/41297181986,
1576001/21651692414,
313326/7428814151, -(450748/2278033113), -(4060782/2406568931),
-(1242769/132435345), -(99750955/2239724446), -(49680141/
241024459),
-(50095271/31604284), 66016483/33048707} *)
`

Then your expression seems to evaluate .

Edit: If you change the parameter of Rationalize to 10^-500 then your sidenote appears to hold.

• thanks for the suggestion. I see, however, negative numbers in your result. The result should be, as I anticipated, a probability between 0 and 1. Maybe there is a problem with the coefficients of the function itself, and the function really turns a bit negative? Commented Aug 4, 2021 at 9:52
• Another strange thing: If I run x = 1/1000; Table[ Rationalize[FD1[x, vE]["WD_O", 1], 10^(-500)], {vE, 10/100, 14/100, 1/100}] Clear[x]; only the first result is not zero (FD1[x, vE]["WD_O", 1] is the function I gave), but if I run FD1[x, vE]["WD_O", 1]//Inputform and copy/paste the output to substitute FD1[x, vE]["WD_O", 1] in the code before I get nonzero results for all entries Commented Aug 4, 2021 at 10:07
• Maybe machine precision of the coefficients is not enough. try to obtain them as exact numbers right away from your integration Commented Aug 4, 2021 at 10:36
• There is something I still don-t get, if I run FD1[1/1000, 14/100]["WD_O", 1] I get 0, while if I run FD1[x, vE]["WD_O", 1] /. x -> 1/1000 /. vE -> 14/100 I get 1.85 Commented Aug 4, 2021 at 11:19