square matrix containing Chebyshev polynomials becomes singular when its size becomes greater than a certain value
Show older comments
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:
, 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
Torsten
on 27 Nov 2023
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.
Torsten
on 28 Nov 2023
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
on 28 Nov 2023
Chang-Yu
on 28 Nov 2023
Chang-Yu
on 28 Nov 2023
Accepted Answer
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
B = zeros(N+1,1);
B(:,1) = (d).^3-(d)+(d).^2;
X= A\B
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
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)
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
B = zeros(N+1,1);
B(:,1) = (d).^3-(d)+(d).^2;
X= A\B
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
More Answers (1)
Bruno Luong
on 27 Nov 2023
0 votes
If you select discrete points d the Chebychev nodes see definition wiki, and not uniform, it will become well posed.
2 Comments
Chang-Yu
on 27 Nov 2023
Bruno Luong
on 27 Nov 2023
Try it, if it works then I give you a mathematic explanation.
Categories
Find more on Polynomials in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
