2
$\begingroup$

I'm trying to plot a function such as:

fwC[k1_, tau_, FE_, COH_, X_, t_] = 1 + (Exp[-k1 t] FE tau (-1 + Exp[k1 t] X (-1 + k1 tau) + Exp[t (k1 - 1/tau)] (1 + X - k1 X tau)))/(COH (-1 + k1 tau))

When I try to plot the function with the values:

Plot[fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, x],{x, 0, 40.}, PlotRange -> All, Frame -> True]

I get the plot:

enter image description here

For values of x > 36. I have the warning "General::munfl: Exp[-803.6] is too small to represent as a normalized machine number; precision may be lost." But writting the equation in 'numerical format' I have:

fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, x] = 1 + 497.869 Exp[-20.09 x] (-1 - 455.546 Exp[20.0856 x] + 456.546 Exp[20.09 x])

Once simplified gives:

f(x)= 227301. - 497.869 Exp[-20.09 x] - 226802. Exp[-0.00439947 x]

Which can be plotted in all the range without any precision problem:

Plot[{227301. - 497.869 Exp[-20.09 x] - 226802.130 Exp[-0.00439947 x],fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, x]}, {x, 0, 1000.}, PlotRange -> All, Frame -> True]

enter image description here

In orange is the function, in blue the numerical simplified expression.

Any help to overcome this kind of ploblems?

Best regards

$\endgroup$
2
  • $\begingroup$ E^(-20.09 x) will underflow when x is greater than Solve[-20.09` x == Log[$MinMachineNumber]]. In general Exp[-k1 t] will underflow when k1 t > -Log[$MinMachineNumber]. $\endgroup$
    – Michael E2
    Aug 15, 2020 at 12:01
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/197758/… $\endgroup$
    – Michael E2
    Aug 15, 2020 at 12:07

2 Answers 2

1
$\begingroup$

It seems to me the OP already has a solution to the problem at hand in the question. It is described in words, but here's an approach to the idea:

Plot[
 fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, x] // Expand // Evaluate,
 {x, 0, 1000.}, PlotRange -> All, Frame -> True]

enter image description here

The problem is a factor that underflows to zero in machine precision. In this case it is the factor Exp[-k1 t] of the second term, which underflows when k1 t is greater than -Log[$MinMachineNumber] == 708.396. When it underflows, the second term will be zero, no matter how big the remaining factors.

Expand distributes the factor and transforms the function expression to a sum of terms, some of which may underflow. Those that underflow are negligible in this form.

$\endgroup$
1
  • $\begingroup$ Thanks @MichaelE2, yes I identified what the problem was, but I don't found how to proceed with Mathematica to distribute the exponential terms in the calculations. The problem was how the calculus proceeds.' Expand' is the shortest way to solve this issue. $\endgroup$ Aug 15, 2020 at 17:30
1
$\begingroup$
Clear["Global`*"]

fwC[k1_, tau_, FE_, COH_, X_, t_] = 
  1 + (Exp[-k1 t] FE tau (-1 + Exp[k1 t] X (-1 + k1 tau) + 
        Exp[t (k1 - 1/tau)] (1 + X - k1 X tau)))/(COH (-1 + k1 tau));

It is a precision issue. To support high precision, Rationalize the function's arguments. Also specify a WorkingPrecision to cause the calculations to be done with arbitrary-precision rather than machine precision.

Plot[Evaluate[
  fwC[k1, tau, FE, COH, X, t] /.
    Thread[{k1, tau, FE, COH, X, t} ->
       {20.09, 227.3, 1000. 10^-8, 
        10^-9, 0.1, x} //
      Rationalize] // FullSimplify],
 {x, 0, 40},
 PlotRange -> All,
 Frame -> True,
 WorkingPrecision -> 25]

enter image description here

In the same way,

Plot[Evaluate[
  fwC[k1, tau, FE, COH, X, t] /.
    Thread[{k1, tau, FE, COH, X, t} ->
       {20.09, 227.3, 1000. 10^-8, 
        10^-9, 0.1, x} //
      Rationalize] // FullSimplify],
 {x, 0, 1000},
 PlotRange -> All,
 Frame -> True,
 WorkingPrecision -> 25]

enter image description here

EDIT: To use this approach more generally, redefine fwC with an optional argument to specify a working precision.

Clear["Global`*"]

fwC[k1_, tau_, FE_, COH_, X_, t_,
   wp_ : MachinePrecision] := Module[{k1p, taup, FEp, COHp, Xp, tp},
   {k1p, taup, FEp, COHp, Xp, tp} =
    If[wp === MachinePrecision,
     {k1, tau, FE, COH, X, t} 
     (* use arguments as given *),
     SetPrecision[{k1, tau, FE, COH, X, t}, wp]
     (* set precision to that specified *)];
   1 + (Exp[-k1p tp] FEp taup (-1 + Exp[k1p tp] Xp (-1 + k1p taup) + 
          Exp[tp (k1p - 1/taup)] (1 + Xp - k1p Xp taup)))/(COHp (-1 + 
          k1p taup)) // Simplify];

Without specifying a working precision (default value of wp, i.e., use precision of arguments as given)

fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, 100.]

(* General::munfl: Exp[-2009.] is too small to represent as a normalized machine number; precision may be lost.

1. *)

% // Precision

(* MachinePrecision *)

With machine precision numbers there is no attempt to track or control precision; you get whatever the machine operations produce.

If the inputs have specified precision or are exact,

fwC[20.09`10, 227.3`20, 1000.0`25 10^-8, 10^-9, 0.1`15, 100.0`15]

(* 81224.5 *)

% // Precision

(* 5.94886 *)

Note that the complexity of the calculation resulted in a loss of precision of about 4.1 digits from the argument with the lowest arbitrary-precision (10).

Specifying a working precision (e.g., wp == 25)

fwC[20.09, 227.3, 1000. 10^-8, 10^-9, 0.1, 100., 25]

(* 81224.455613146224781 *)

% // Precision

(* 20.6477 *)

Note that the complexity of the calculation resulted in a loss of precision of about 4.4 digits from the specified precision (25).

$\endgroup$
3
  • $\begingroup$ Thank you Bob for your answer. Yes, it is a precision problem and how Mathematica works inside. What I need is the evaluation of the function under any circumstance, because I will use fwC[ ] in other Mathematica functions such as NonlinearModelFit[ ]. I fixed the problem using: ` fwC[k1_, tau_, FE_, COH_, X_, t_] = 1 + (Exp[-k1 t] FE tau (-1 + Exp[k1 t] X (-1 + k1 tau) + Exp[t (k1 - 1/tau)] (1 + X - k1 X tau)))/(COH (-1 + k1 tau))//Expand; ` The issue I think is related on how Mathematica manages the internal calculus, the function is not so complex once simplified. $\endgroup$ Aug 15, 2020 at 11:36
  • $\begingroup$ @JavierJNavarro Welcome to numerical analysis. This isn't a Mathematica issue. When approximating, you must formulate the problem in a way that the automated arithmetic you're using can work in your domain. "Under any circumstance" is impossible with a computer. $\endgroup$
    – John Doty
    Aug 15, 2020 at 14:06
  • $\begingroup$ @bobhanlon thank you a lot for your suggestion. I this way I will have more control on the evaluation of my functions. $\endgroup$ Aug 15, 2020 at 17:39

Your Answer

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

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