UPD1: I loaded original data set in .txt format. There are two columns in the file: first for x values and second for phi[x] values.
Preface
The question arised during the work on a computational mathematics home task in my university. Please note that here in the preface I explain only the background of my question. For the question itself see the following chapters. The task was simple: based on a set ob "observed" values of an unknown function phi[x] construct a polynomial approximation of this function:
f[x_] := Sum[c[[i + 1]] * P[i, x], {i, 0, M}]
Here c is a list of expansion coefficients, P[i, x] is a polynomial of i-th degree of x (from any orthogonal system, e.g. Legendre's polynomial) and M is some pre-defined (set by the instructor) number. We had to find coefficients from the list c.
If the unknown function phi[x] was measured in N points from a list x = {x1, x2, ..., xN} then this problem can be re-written in form of reduntant system of linear equations with N rows and M columns:
c[[1]] * P[0, x1] + c[[2]] * P[1, x1] + ... + c[[M + 1]] * P[M, x1] == phi[x1];
c[[1]] * P[0, x2] + c[[2]] * P[1, x2] + ... + c[[M + 1]] * P[M, x2] == phi[x2];
...
c[[1]] * P[0, xN] + c[[2]] * P[1, xN] + ... + c[[M + 1]] * P[M, xN] == phi[xN];
Or in matrix form:
Dot[A, c] == phi[x]
We had to solve this system with two methods: using QR-decomposition of the matrix A and simply multiplying both sides by Transpose[A] and then solving a new system with a square positive-defined matrix using, for example, Cholesky method:
Dot[Transpose[A], A].c == Dot[Transpose[A], phi[x]]
The problem
All the calculations were not difficult and the resulting approximation was quite good. As a part of the task report we had to compute condition number of both matrices A and Transpose[A].A in spectral norm:
condA = First[#] / Last[#]&@SingularValueList[A];
condATA = First[#] / Last[#]&@SingularValueList[Traspose[A].A]
And here quite surprising things happened. According to "obvious" (and naive actually) estimate condATA should grow with increase of the maximal polynomial degree as square of condA:
condATA ~ condA^2 (*as max degree M grows*)
And this rule was satisfied for small (1 - 7) values of M. But then condATA stopped growing and finally for M ~ 10-20 both condition numbers started oscillating around level of 10^9 - 10^13.
Since we were not restricted in the programming language that we use in this task, a lot of my groupmates chose Python as their working tool. And in Python the rule mentioned above worked even for big M: condATA grew as square of condA reaching values of 10^18 and higher. MATLAB computations gave similar results.
The question
So, the question is: which results are correct - those obtained in Mathematica or those computed in Python's numpy library and MATLAB? I worked in machine precision, accuracy of the condition number condATA varied from ~11 for small degrees to ~2.6 for M = 10 (I mean results of the Accuracy[] function). No warnings arised during its computation.
And if Mathematica is wrong how can I avoid such situations in future?
Some words in the end
This is my first post here in the Wolfram Community, so if you have suggestions on the style please write them in the comments. I also apologize for all possible mistakes in my English.
Attachments:
