square matrix containing Chebyshev polynomials becomes singular when its size becomes greater than a certain value

Hi, I've been trying to approximate functions with Chebyshev polynomials. In this case, I'm looking at , and I've approximated as where are coefficients and is the k-th Chebyshev polynomial. I've got N+1 points inside the interval [0,4] such that . My goal is to find what the coefficients are, and I did this by creating a matrix equation of the form and solving for x, where x is a vector of length N+1 containing the coefficients . So I applied the previous Chebyshev approximation on all N+1 points and got a matrix equation where A = , x = and B= , and I solved for x. Now, in order to simplify things, I chose my points to be evenly spaced within [0,4]. When trying out my code, I would get a certain vector for x and compare it via chebcoeffs with the actual coefficients. My estimated coefficients from x got increasingly closer to those obtained from chebcoeffs as I increased N starting from 1, as one would expect. However, when I hit N = 56 and onwards, my matrix A began to not be full rank, and hence be singular. I'm not sure why this happens, even though my definition of all matrices and vectors appears correct from a theoretical standpoint. Here's the full code:

t0 = 0;
tf = 4;
N = 56;
t = chebfun('t',[t0 tf]);

A = zeros(N+1,N+1);
d = zeros(N+1,1);

for c=0:N
    d(c+1,1)=tf*c./(N);
end

h = chebpoly(0:N,[t0 tf]);

for r=0:N
    A(r+1,:)= h(d(r+1,1));
end
A

B = zeros(N+1,1)
B(:,1) = (d).^3-(d)+(d).^2;
B

X= A\B

b=chebcoeffs((t).^3-(t)+(t).^2) %for comparison with actual chebyshev coefficients
rank(A)

[Torsten]

Why do you use the Chebyshev polynomials for approximation as if they were the usual polynomial basis 1,x,x^2,...,x^n ? Chebyshev polynomials are orthogonal with respect to a special scalar product, and this property is used to determine the coefficients a_n in your above expansion.

[Chang-Yu]

this is for a project of mine in which I need to find the coefficients a_n, which are approximates to the actual chebyshev coefficients used. If you read the Chebfun handbook available online, then you'd see that this sort of approximation is reasonable. I just dont get why it doesnt work for values of N greater than 56

[Torsten]

If you use the MATLAB function "chebyshevT" to create the polynomials and convert them to numerical functions using "matlabFunction", your code also works for n > 56.

[Walter Roberson]

Are you using https://www.mathworks.com/matlabcentral/fileexchange/47023-chebfun-current-version ?

[Chang-Yu]

@Torsten the function chebyshevT only creates chebyshev functions on the interval [-1,1], but the chebyshev functions that I'm using are on the interval [0,4] (this can be done via chebpoly, which scales them appropiately). Do you know how to make chebyshev functions scaled onto the interval [0 4] through chebyshevT?

[Chang-Yu]

@Walter Roberson yeah I'm using the current version and it still doesn't work for N > 56.

Answer 1

[Torsten]

t0 = 0;
tf = 4;
N = 75;
syms x
h = chebyshevT(0:N, x);

A = zeros(N+1,N+1);
d = zeros(N+1,1);
for c=1:N+1
    d(c,1) = 0.5*(t0+tf)+0.5*(tf-t0)*cos((2*c-1)/(2*(N+1))*pi);
end

for c=1:N+1
    A(c,:)= vpa(subs(h,x,cos((2*c-1)/(2*(N+1))*pi)));
end
A

A = 76×76

1.0000 0.9998 0.9991 0.9981 0.9966 0.9947 0.9923 0.9896 0.9864 0.9827 0.9787 0.9743 0.9694 0.9641 0.9584 0.9523 0.9458 0.9389 0.9316 0.9239 0.9158 0.9073 0.8984 0.8891 0.8795 0.8694 0.8591 0.8483 0.8372 0.8257 1.0000 0.9981 0.9923 0.9827 0.9694 0.9523 0.9316 0.9073 0.8795 0.8483 0.8138 0.7763 0.7357 0.6923 0.6463 0.5978 0.5469 0.4940 0.4392 0.3827 0.3247 0.2655 0.2052 0.1442 0.0826 0.0207 -0.0413 -0.1032 -0.1646 -0.2254 1.0000 0.9947 0.9787 0.9523 0.9158 0.8694 0.8138 0.7496 0.6773 0.5978 0.5119 0.4205 0.3247 0.2254 0.1237 0.0207 -0.0826 -0.1849 -0.2853 -0.3827 -0.4759 -0.5641 -0.6463 -0.7216 -0.7891 -0.8483 -0.8984 -0.9389 -0.9694 -0.9896 1.0000 0.9896 0.9584 0.9073 0.8372 0.7496 0.6463 0.5295 0.4017 0.2655 0.1237 -0.0207 -0.1646 -0.3051 -0.4392 -0.5641 -0.6773 -0.7763 -0.8591 -0.9239 -0.9694 -0.9947 -0.9991 -0.9827 -0.9458 -0.8891 -0.8138 -0.7216 -0.6142 -0.4940 1.0000 0.9827 0.9316 0.8483 0.7357 0.5978 0.4392 0.2655 0.0826 -0.1032 -0.2853 -0.4577 -0.6142 -0.7496 -0.8591 -0.9389 -0.9864 -0.9998 -0.9787 -0.9239 -0.8372 -0.7216 -0.5811 -0.4205 -0.2455 -0.0620 0.1237 0.3051 0.4759 0.6304 1.0000 0.9743 0.8984 0.7763 0.6142 0.4205 0.2052 -0.0207 -0.2455 -0.4577 -0.6463 -0.8017 -0.9158 -0.9827 -0.9991 -0.9641 -0.8795 -0.7496 -0.5811 -0.3827 -0.1646 0.0620 0.2853 0.4940 0.6773 0.8257 0.9316 0.9896 0.9966 0.9523 1.0000 0.9641 0.8591 0.6923 0.4759 0.2254 -0.0413 -0.3051 -0.5469 -0.7496 -0.8984 -0.9827 -0.9966 -0.9389 -0.8138 -0.6304 -0.4017 -0.1442 0.1237 0.3827 0.6142 0.8017 0.9316 0.9947 0.9864 0.9073 0.7631 0.5641 0.3247 0.0620 1.0000 0.9523 0.8138 0.5978 0.3247 0.0207 -0.2853 -0.5641 -0.7891 -0.9389 -0.9991 -0.9641 -0.8372 -0.6304 -0.3635 -0.0620 0.2455 0.5295 0.7631 0.9239 0.9966 0.9743 0.8591 0.6619 0.4017 0.1032 -0.2052 -0.4940 -0.7357 -0.9073 1.0000 0.9389 0.7631 0.4940 0.1646 -0.1849 -0.5119 -0.7763 -0.9458 -0.9998 -0.9316 -0.7496 -0.4759 -0.1442 0.2052 0.5295 0.7891 0.9523 0.9991 0.9239 0.7357 0.4577 0.1237 -0.2254 -0.5469 -0.8017 -0.9584 -0.9981 -0.9158 -0.7216 1.0000 0.9239 0.7071 0.3827 -0.0000 -0.3827 -0.7071 -0.9239 -1.0000 -0.9239 -0.7071 -0.3827 0.0000 0.3827 0.7071 0.9239 1.0000 0.9239 0.7071 0.3827 -0.0000 -0.3827 -0.7071 -0.9239 -1.0000 -0.9239 -0.7071 -0.3827 0.0000 0.3827

B = zeros(N+1,1);
B(:,1) = (d).^3-(d)+(d).^2;

X= A\B

X = 76×1

24.0000 36.0000 14.0000 2.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000

hold on
plot(d,(d).^3-(d)+(d).^2)
fplot(subs(h,x,(x-0.5*(t0+tf))/(0.5*(tf-t0)))*X,[t0 tf])
hold off
grid on

[Bruno Luong]

That is exactky what I have suggested in my answer: use the Chebyschev nodes

t0 = 0;
tf = 4;
N = 75;
syms x
h = chebyshevT(0:N, x);
A = zeros(N+1,N+1);
for c=1:N+1
    A(c,:)= vpa(subs(h,x,cos((2*c-1)/(2*(N+1))*pi)));
end

and the matrix is absolutely well conditioned for any (large) N

cond(A)

ans = 1.4142

[Torsten]

In comparison to equidistant points:

t0 = 0;
tf = 4;
N = 75;
syms x
h = chebyshevT(0:N, x);
A = zeros(N+1,N+1);
d = zeros(N+1,1);
for c=0:N
    d(c+1,1) = 0.5*(t0+tf)+0.5*(tf-t0)*(-1+2*c/N);     
end
for c=0:N
    A(c+1,:)= vpa(subs(h,x,-1+2*c/N));
end
A

A = 76×76

1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 1.0000 -0.9733 0.8948 -0.7685 0.6012 -0.4018 0.1811 0.0494 -0.2772 0.4902 -0.6771 0.8278 -0.9344 0.9912 -0.9951 0.9460 -0.8463 0.7016 -0.5194 0.3095 -0.0832 -0.1477 0.3706 -0.5738 0.7464 -0.8791 0.9650 -0.9994 0.9805 -0.9094 1.0000 -0.9467 0.7924 -0.5535 0.2557 0.0695 -0.3872 0.6636 -0.8693 0.9822 -0.9903 0.8929 -0.7001 0.4327 -0.1192 -0.2071 0.5113 -0.7609 0.9294 -0.9988 0.9616 -0.8218 0.5944 -0.3036 -0.0196 0.3408 -0.6255 0.8435 -0.9716 0.9960 1.0000 -0.9200 0.6928 -0.3548 -0.0401 0.4285 -0.7483 0.9484 -0.9968 0.8857 -0.6329 0.2788 0.1199 -0.4994 0.7990 -0.9708 0.9872 -0.8457 0.5688 -0.2010 -0.1990 0.5672 -0.8446 0.9869 -0.9712 0.8002 -0.5012 0.1219 0.2768 -0.6313 1.0000 -0.8933 0.5961 -0.1717 -0.2894 0.6887 -0.9411 0.9927 -0.8325 0.4948 -0.0515 -0.4028 0.7712 -0.9750 0.9709 -0.7596 0.3863 0.0695 -0.5104 0.8424 -0.9947 0.9348 -0.6755 0.2721 0.1894 -0.6104 0.9013 -0.9998 0.8851 -0.5815 1.0000 -0.8667 0.5022 -0.0039 -0.4955 0.8628 -1.0000 0.8705 -0.5089 0.0116 0.4888 -0.8589 0.9999 -0.8743 0.5155 -0.0193 -0.4821 0.8549 -0.9997 0.8780 -0.5221 0.0270 0.4753 -0.8509 0.9995 -0.8816 0.5286 -0.0347 -0.4685 0.8468 1.0000 -0.8400 0.4112 0.1492 -0.6618 0.9627 -0.9555 0.6425 -0.1240 -0.4343 0.8535 -0.9997 0.8259 -0.3879 -0.1743 0.6807 -0.9693 0.9477 -0.6228 0.0987 0.4571 -0.8665 0.9987 -0.8113 0.3643 0.1993 -0.6991 0.9752 -0.9392 0.6027 1.0000 -0.8133 0.3230 0.2879 -0.7913 0.9993 -0.8342 0.3577 0.2524 -0.7682 0.9973 -0.8540 0.3919 0.2165 -0.7441 0.9939 -0.8726 0.4256 0.1803 -0.7189 0.9891 -0.8901 0.4587 0.1439 -0.6928 0.9830 -0.9063 0.4912 0.1073 -0.6657 1.0000 -0.7867 0.2377 0.4127 -0.8870 0.9829 -0.6594 0.0545 0.5736 -0.9569 0.9320 -0.5094 -0.1305 0.7148 -0.9941 0.8492 -0.3420 -0.3111 0.8315 -0.9971 0.7373 -0.1629 -0.4810 0.9196 -0.9659 0.6001 0.0218 -0.6344 0.9763 -0.9017 1.0000 -0.7600 0.1552 0.5241 -0.9518 0.9227 -0.4506 -0.2377 0.8119 -0.9965 0.7027 -0.0716 -0.5938 0.9742 -0.8870 0.3740 0.3185 -0.8581 0.9859 -0.6404 -0.0125 0.6594 -0.9897 0.8450 -0.2947 -0.3971 0.8983 -0.9683 0.5735 0.0965

B = zeros(N+1,1);
B(:,1) = (d).^3-(d)+(d).^2;
X= A\B

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.466622e-19.

X = 76×1

14.4987 29.0637 -4.5242 -4.4614 -17.1230 -5.5462 -14.8989 -4.2590 -12.0098 -2.6990

hold on
plot(d,(d).^3-(d)+(d).^2)
fplot(subs(h,x,(x-0.5*(t0+tf))/(0.5*(tf-t0)))*X,[t0 tf])
hold off
grid on

Answer 2

[Bruno Luong]

If you select discrete points d the Chebychev nodes see definition wiki, and not uniform, it will become well posed.

[Chang-Yu]

if you read the Chebfun handbook available online, there it says you can use any points on the interval. Its just that I'm not sure why it becomes singular after a certain N, but for N less than 56, it gives good values.

[Bruno Luong]

Try it, if it works then I give you a mathematic explanation.