Mercurial > matrix-functions
comparison toolbox/chebspec.m @ 2:c124219d7bfa draft
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Thu, 07 May 2015 18:36:24 +0200 |
| parents | 8f23314345f4 |
| children |
comparison
equal
deleted
inserted
replaced
| 1:e471a92d17be | 2:c124219d7bfa |
|---|---|
| 1 function C = chebspec(n, k) | |
| 2 %CHEBSPEC Chebyshev spectral differentiation matrix. | |
| 3 % C = CHEBSPEC(N, K) is a Chebyshev spectral differentiation | |
| 4 % matrix of order N. K = 0 (the default) or 1. | |
| 5 % For K = 0 (`no boundary conditions'), C is nilpotent, with | |
| 6 % C^N = 0 and it has the null vector ONES(N,1). | |
| 7 % C is similar to a Jordan block of size N with eigenvalue zero. | |
| 8 % For K = 1, C is nonsingular and well-conditioned, and its eigenvalues | |
| 9 % have negative real parts. | |
| 10 % For both K, the computed eigenvector matrix X from EIG is | |
| 11 % ill-conditioned (MESH(REAL(X)) is interesting). | |
| 12 | |
| 13 % References: | |
| 14 % C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral | |
| 15 % Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1988; p. 69. | |
| 16 % L.N. Trefethen and M.R. Trummer, An instability phenomenon in | |
| 17 % spectral methods, SIAM J. Numer. Anal., 24 (1987), pp. 1008-1023. | |
| 18 % D. Funaro, Computing the inverse of the Chebyshev collocation | |
| 19 % derivative, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 1050-1057. | |
| 20 | |
| 21 if nargin == 1, k = 0; end | |
| 22 | |
| 23 % k = 1 case obtained from k = 0 case with one bigger n. | |
| 24 if k == 1, n = n + 1; end | |
| 25 | |
| 26 n = n-1; | |
| 27 C = zeros(n+1); | |
| 28 | |
| 29 one = ones(n+1,1); | |
| 30 x = cos( (0:n)' * (pi/n) ); | |
| 31 d = ones(n+1,1); d(1) = 2; d(n+1) = 2; | |
| 32 | |
| 33 % eye(size(C)) on next line avoids div by zero. | |
| 34 C = (d * (one./d)') ./ (x*one'-one*x' + eye(size(C))); | |
| 35 | |
| 36 % Now fix diagonal and signs. | |
| 37 | |
| 38 C(1,1) = (2*n^2+1)/6; | |
| 39 for i=2:n+1 | |
| 40 if rem(i,2) == 0 | |
| 41 C(:,i) = -C(:,i); | |
| 42 C(i,:) = -C(i,:); | |
| 43 end | |
| 44 if i < n+1 | |
| 45 C(i,i) = -x(i)/(2*(1-x(i)^2)); | |
| 46 else | |
| 47 C(n+1,n+1) = -C(1,1); | |
| 48 end | |
| 49 end | |
| 50 | |
| 51 if k == 1 | |
| 52 C = C(2:n+1,2:n+1); | |
| 53 end |