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

4 views (last 30 days)
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)
  6 Comments
Chang-Yu
Chang-Yu on 28 Nov 2023
@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?

Sign in to comment.

Accepted Answer

Torsten
Torsten on 28 Nov 2023
Edited: Torsten on 28 Nov 2023
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
  2 Comments
Bruno Luong
Bruno Luong on 28 Nov 2023
Moved: Bruno Luong on 28 Nov 2023
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
Torsten on 28 Nov 2023
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

Sign in to comment.

More Answers (1)

Bruno Luong
Bruno Luong on 27 Nov 2023
If you select discrete points d the Chebychev nodes see definition wiki, and not uniform, it will become well posed.
  2 Comments
Chang-Yu
Chang-Yu on 27 Nov 2023
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.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!