square matrix containing Chebyshev polynomials becomes singular when its size becomes greater than a certain value
9 views (last 30 days)
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:
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
Accepted Answer
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
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
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)
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
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
If you select discrete points d the Chebychev nodes see definition wiki, and not uniform, it will become well posed.
2 Comments
See Also
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!