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 

If you see the above result carefully, you will find the result is wrong. So I use the Matlab to calculate it again. The results are the same.

We can see that both MMA and Matlab did not give a good answer.

Let us expand the first formula at small $ u$

Series[(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, {u, 0, 10}] 

Until now, I found two types of result.

Result A:

 1 - 16/5 Sqrt[2/\[Pi]] Sqrt[u] +   1/15 (-37 - 105 EulerGamma + 105 Log[2] - 105 Log[u]) u +   1/60 (-257 - 630 EulerGamma + 630 Log[2] - 630 Log[u]) u^2 +   1/60 (-16 - 693 EulerGamma + 693 Log[2] - 693 Log[u]) u^3 + ((7519 -      25740 EulerGamma + 25740 Log[2] -      25740 Log[u]) u^4)/2880 + ((3275 - 6006 EulerGamma +      6006 Log[2] - 6006 Log[u]) u^5)/1152 + ((52483 -      69615 EulerGamma + 69615 Log[2] -      69615 Log[u]) u^6)/28800 + ((2082133 - 2238390 EulerGamma +      2238390 Log[2] - 2238390 Log[u]) u^7)/2419200 + ((25195103 -      23279256 EulerGamma + 23279256 Log[2] -      23279256 Log[u]) u^8)/77414400 + ((95727847 -      78738660 EulerGamma + 78738660 Log[2] -      78738660 Log[u]) u^9)/928972800 + ((156692681 -      117417300 EulerGamma + 117417300 Log[2] -      117417300 Log[u]) u^10)/5573836800 

Resut B:

(1 - 16/5 Sqrt[2 \[Pi]] Sqrt[u] +    1/15 (-37 - 105 EulerGamma + 105 Log[2] - 105 Log[u]) u +    1/60 (-257 - 630 EulerGamma + 630 Log[2] - 630 Log[u]) u^2 +    1/60 (-16 - 693 EulerGamma + 693 Log[2] -       693 Log[u]) u^3 + ((7519 - 25740 EulerGamma + 25740 Log[2] -       25740 Log[u]) u^4)/   2880 + ((3275 - 6006 EulerGamma + 6006 Log[2] - 6006 Log[u]) u^5)/   1152 + ((52483 - 69615 EulerGamma + 69615 Log[2] -       69615 Log[u]) u^6)/   28800 + ((2082133 - 2238390 EulerGamma + 2238390 Log[2] -       2238390 Log[u]) u^7)/   2419200 + ((25195103 - 23279256 EulerGamma + 23279256 Log[2] -       23279256 Log[u]) u^8)/   77414400 + ((95727847 - 78738660 EulerGamma + 78738660 Log[2] -       78738660 Log[u]) u^9)/   928972800 + ((156692681 - 117417300 EulerGamma + 117417300 Log[2] -       117417300 Log[u]) u^10)/5573836800) 

You can check that one of them is wrong by checking the value at a specific $ u$ and MMA did not give a valid result for the first formula, i.e.,

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}] 

I believe the correct result is

Indeterminate,0.793116,0.725224,0.67821,0.641557,0.611302,0.585457,0.562865,0.542783,0.524707,0.508273,0.493211,0.479315,0.466422,0.454402,0.44315,0.432578,0.422613,0.413194,0.404268,0.39579,0.387721,0.380025,0.372674,0.36564,0.3589,0.352432,0.346218,0.340241,0.334485,0.328937,0.323584,0.318414,0.313416,0.308581,0.303901,0.299366,0.294969,0.290704,0.286564,0.282542,0.278633,0.274832,0.271134,0.267535,0.264029,0.260614,0.257284,0.254038,0.250871,0.24778,0.244762,0.241815,0.238936,0.236122,0.233371,0.23068,0.228048,0.225473,0.222952,0.220484,0.218066,0.215698,0.213378,0.211104,0.208874,0.206688,0.204544,0.202441,0.200377,0.198351,0.196362,0.19441,0.192493,0.19061,0.18876,0.186942,0.185156,0.1834,0.181675,0.179978,0.178309,0.176667,0.175053,0.173464,0.171901,0.170363,0.168849,0.167359,0.165892,0.164447,0.163024,0.161623,0.160242,0.158883,0.157543,0.156222,0.154921,0.153639,0.152375,0.151129