3
$\begingroup$

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:

Indeterminate,0.966028,0.96976,0.977704,0.987383,0.997947,1.00901,1.02035,1.03185,1.04345,1.05507,1.0667,1.0783,1.08987,1.10138,1.11284,1.12423,1.13555,1.1468,1.15798,1.16908,1.18011,1.19106,1.20193,1.21274,1.22346,1.23412,1.2447,1.25521,1.26565,1.27602,1.28632,1.29656,1.30672,1.31683,1.32687,1.33684,1.34676,1.35661,1.3664,1.37614,1.38581,1.39543,1.405,1.41451,1.42396,1.43336,1.44271,1.45201,1.46126,1.47046,1.47961,1.48871,1.49776,1.50677,1.51573,1.52464,1.53351,1.54234,1.55112,1.55986,1.56856,1.57721,1.58583,1.59441,1.60294,1.61144,1.6199,1.62832,1.6367,1.64504,1.65335,1.66162,1.66986,1.67806,1.68623,1.69436,1.70246,1.71052,1.71856,1.72656,1.73452,1.74246,1.75036,1.75824,1.76608,1.77389,1.78167,1.78942,1.79715,1.80484,1.81251,1.82014,1.82775,1.83533,1.84289,1.85041,1.85791,1.86539,1.87283,1.88026

Do you think the above numerical result is correct?

$\endgroup$
  • $\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 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 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 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$ – Eric Towers Jul 20 at 22:47
9
$\begingroup$

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]

1

You could fill in that hole by defining:

PHI30EQ[_?(EqualTo[0])] = 1;
| improve this answer | |
$\endgroup$
  • $\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 at 11:07
  • $\begingroup$ But the other values for different u are not correct. $\endgroup$ – Blueka Jul 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$ – Sjoerd Smit Jul 20 at 11:12
  • $\begingroup$ Thank you. I get it! $\endgroup$ – Blueka Jul 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.