forked from chebfun/chebfun
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cheb2jac.m
145 lines (113 loc) · 5.32 KB
/
cheb2jac.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
function c_jac = cheb2jac( c_cheb, alpha, beta )
%CHEB2JAC Convert Chebyshev coefficients to Jacobi coefficients.
% C_JAC = CHEB2LEG(C_CHEB, A, B) converts the vector C_CHEB of Chebyshev
% coefficients to a vector C_JAC of Jacobi coefficients such that
% C_CHEB(1)*T_0(x) + ... + C_CHEB(N)*T{N-1}(x) = ...
% C_LEG(1)*P_0^{(A,B)}(x) + ... + C_LEG(N)*P{N-1}^{(A,B)}(x),
% where P_k^{(A,B)} is the degree k Jacobi polynomial corresponding to the
% weight function w(x) = (1-X)^A * (1+X)^B.
%
% See also JAC2CHEB, JAC2JAC.
% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.
%%%%%%%%% For developers %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For more information on the algorithm for N>512, see Section 5.3 of
%
% A. Townsend, M. Webb, and S. Olver, "Fast polynomial transforms based on
% Toeplitz and Hankel matrices", submitted, 2016.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = size(c_cheb, 1); % Length of coefficients.
if ( alpha == 0 && beta == 0 )
% Use cheb2leg for alpha = beta = 0:
c_jac = cheb2leg( c_cheb );
elseif ( alpha == -.5 && beta == -.5 )
% Undo scaling:
% Convert T_n -> P_n^(-1/2,1/2):
scl = [1 cumprod((1/2:1/2+N-2)./(1:N-1))]'; % P_n^(-1/2,-1/2)(1)
c_jac = spdiags( scl, 0, N, N) \ c_cheb;
elseif ( N <= 512 )
c_jac = cheb2jac_direct(c_cheb, alpha, beta);
else
% Convert Chebyshev to Jacobi (-1/2,-1/2) and then call jac2jac:
% Convert T_n -> P_n^(-1/2,1/2):
scl = [1 cumprod((1/2:1/2+N-2)./(1:N-1))]'; % P_n^(-1/2,-1/2)(1)
c_jac = spdiags( scl, 0, N, N) \ c_cheb;
% Now convert P_n^(-1/2,-1/2) -> P_n^(alpha,beta):
c_jac = jac2jac( c_jac, -1/2, -1/2, alpha, beta );
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% DIRECT METHOD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c_jac = cheb2jac_direct(c_cheb, a, b)
%CHEB2JAC_DIRECT Convert Cheb to Jacobi (A,B) coeffs using direct method.
[N, m] = size(c_cheb); % Number of columns.
N = N - 1; % Degree of polynomial.
% Don't let N be too big:
if ( N <= 0 ), c_jac = c_cheb; return, end % Trivial case.
f = chebtech2.coeffs2vals([c_cheb ; zeros(N, m)]); % Values on 2*N+1 Cheb grid.
% 2*N+1 Chebyshev grid (reversed order) and Clenshaw-Curtis-Jacobi weights:
[w, x] = ccjQuadwts(2*N+1, a, b);
% Make the Jacobi-Chebyshev Vandemonde matrix:
apb = a + b; aa = a * a; bb = b * b;
P = zeros(2*N+1, N+1); P(:,1) = 1;
P(:,2) = 0.5*(2*(a + 1) + (apb + 2)*(x - 1));
for k = 2:N
k2 = 2*k;
k2apb = k2 + apb;
q1 = k2*(k + apb)*(k2apb - 2);
q2 = (k2apb - 1)*(aa - bb);
q3 = (k2apb - 2)*(k2apb - 1)*k2apb;
q4 = 2*(k + a - 1)*(k + b - 1)*k2apb;
P(:,k+1) = ((q2 + q3*x).*P(:,k) - q4*P(:,k-1)) / q1;
end
% Scaling:
% NN = (0:N)';
% scale = 2^(a+b+1)*gamma(NN+a+1).*gamma(NN+b+1) ./ ...
% ((2*NN+a+b+1).*gamma(NN+a+b+1).*factorial(NN))
scale = zeros(N+1, 1);
scale(1) = beta(a+1, b+1);
for n = 0:N-1
scale(n+2) = (2*n+a+b+1)*(n+a+1)*(n+b+1) / ...
((n+1)*(2*n+a+b+3)*(n+a+b+1))*scale(n+1);
end
scale = 2^(a+b+1)*scale;
% Jacobi coefficients:
c_jac = bsxfun(@times, P.'*(bsxfun(@times, f , w.')), 1./scale);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% DCT METHODS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [w, x] = ccjQuadwts(n, a, b)
%CCJQUADWTS Clenshaw-Curtis-Jacobi quadrature weights.
% [W, X] = CCJQUADWTS(N, A, B) returns the N-point Clenshaw-Curtis-Jacobi
% quadrature nodes, X = CHEBPTS(N), and weights, W, corresponding to the
% weight function w(t) = (1-t)^A * (1+t)^B on the interval [-1,1].
% TODO: Move this to somewhere more sensible / accessible.
if ( a == b && a == 0 ) % Clenshaw-Curtis
c = 2./[1, 1-(2:2:(n-1)).^2]; % Standard Chebyshev moments
c = [c, c(floor(n/2):-1:2)]; % Mirror for DCT via FFT
w = ifft(c); % Interior weights
w([1, n]) = w(1)/2; % Boundary weights
elseif ( a == b ) % Gegenbauer
l = a + .5; % Gegenbauer parameter
g0 = gamma(l+.5)*sqrt(pi)/gamma(l+1);
k = 1:floor((n-1)/2);
c = g0*[1, cumprod((k-l-1)./(k+l))]; % Chebyshev moments for (1-x)^a(1+x)^b
c = [c, c(floor(n/2):-1:2)]; % Mirror for DCT via FFT
w = ifft(c); % Interior weights
w([1, n]) = w(1)/2; % Boundary weights
else % Jacobi
c = [1, (a-b)/(a+b+2), zeros(1, n-2)]; % Initialise moments
for r = 1:n % Recurrence relation for 3F2([r, -r, b +1 ; .5, a+b+2] ; 1 ):
c(r+2) = - (2*(b-a)*c(r+1) + (a+b+2-r)*c(r)) / (a+b+2+r);
end
c = 2^(a+b+1)*gamma(a+1)*gamma(b+1)/gamma(a+b+2) * c; % Moments (with const)
v = ifft([c(1:n), c(n-1:-1:2)]); % Mirror for DCT via FFT
w = [v(1), 2*v(2:n-1), v(n)]; % Rescale interior weights
end
if ( nargout > 1 )
x = chebtech2.chebpts(n); % 2nd-kind Chebyshev points.
end
end