# Some Strange results from MMA and Matlab

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?

• 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? Jul 20 '20 at 8:43
• @xzczd "Until now, I found two types of results." I mean I use two versions of MMA, they give me different answers. Jul 20 '20 at 8:50
• @xzczd At the first stage, I think the authors make a mistake. I am sorry, now I changed my mind. Jul 20 '20 at 8:52
• Why is Matlab mentioned in the title? It is not mentioned or referenced in the text of the Question. 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]


1

You could fill in that hole by defining:

PHI30EQ[_?(EqualTo[0])] = 1;

• Since the above formula is the Laplace transform of a normailzed PDF. Thus at $u=0$, the corresponding value is 1. Jul 20 '20 at 11:07
• But the other values for different u are not correct. Jul 20 '20 at 11:09
• 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. Jul 20 '20 at 11:12
• Thank you. I get it! Jul 20 '20 at 11:28