Recently I am repeating a numerical discrete Fourier transform. Since this problem is related to a real experiment.

So the authors use some complex functions to fit the results.

Below formula is the Laplace transform of the waiting time PDF:

PHI30EQ[u_] := (u - 52/15*u^2 - 16/15*u^3)*Exp[u]*
    BesselK[1, u] + (7*u + 4*u^2 + 16/15*u^3)*Exp[u]*
    BesselK[0, u] - (2*u/Pi)^(1/2)*16/5;

Table[N[PHI30EQ[u], 30], {u, 0, 0.1, 0.001}]

The aim is to calculate the above formula numerically. The following information can be ignored.

The results given by MMA are:


Do you think the above numerical result is correct?

  • $\begingroup$ What do you mean by "Until now, I found two types of result."? Why do you believe it's Mathematica and Matlab that's wrong, rather than the formula itself? $\endgroup$
    – xzczd
    Jul 20 '20 at 8:43
  • $\begingroup$ @xzczd "Until now, I found two types of results." I mean I use two versions of MMA, they give me different answers. $\endgroup$
    – Blueka
    Jul 20 '20 at 8:50
  • $\begingroup$ @xzczd At the first stage, I think the authors make a mistake. I am sorry, now I changed my mind. $\endgroup$
    – Blueka
    Jul 20 '20 at 8:52
  • 1
    $\begingroup$ Why is Matlab mentioned in the title? It is not mentioned or referenced in the text of the Question. $\endgroup$ Jul 20 '20 at 22:47

The line Table[N[PHI30EQ[u], 30], {u, 0, 0.1, 0.001}] doesn't do what you think it does. You're asking for 30 digits of precision, but you supply u as a machine number. If you mix arbitrary precision and machine precision like that, you'll get machine precision answers. I suspect you instead want is:

Table[N[PHI30EQ[u], 30], {u, 0, Rationalize[0.1], Rationalize[0.001]}]

{Indeterminate, 0.966028446491522706752530103591, 0.969759610917190434881076174700,...}

It also looks like this formula goes to 1 in the limit u -> 0:

Limit[PHI30EQ[u], u -> 0]


You could fill in that hole by defining:

PHI30EQ[_?(EqualTo[0])] = 1;
  • $\begingroup$ Since the above formula is the Laplace transform of a normailzed PDF. Thus at $u=0$, the corresponding value is 1. $\endgroup$
    – Blueka
    Jul 20 '20 at 11:07
  • $\begingroup$ But the other values for different u are not correct. $\endgroup$
    – Blueka
    Jul 20 '20 at 11:09
  • $\begingroup$ There is little reason to doubt the numerical accuracy of Mathematica's result in this case. When you do arbitrary precision calculations, the rounding errors are tracked internally and the formula is really not that complicated. If something is wrong, it's probably the formula itself. $\endgroup$ Jul 20 '20 at 11:12
  • $\begingroup$ Thank you. I get it! $\endgroup$
    – Blueka
    Jul 20 '20 at 11:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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