Example 6.10

Calculation of derivatives using discrete Chebychev Transform

Contents

Initialization
Using Chebyshev Polynomial :
Using Finite differences :

Initialization

close all;
clear example N u du j X1 U1 a b j Ud1
clear dX X2 U2 A Ud2 x

example = 'b';


N = 5;


switch example
    case 'a'
        u  = @(x) x.^4;
        du = @(x) 4*x.^3;

    case 'b'
        u  = @(x) 4*(x.^2).*(1 - x.^2).*exp(-x/2);
        du = @(x) (8*x - 2*x.^2 - 16*x.^3 + 2*x.^4).*exp(-x/2);
end;

Using Chebyshev Polynomial :

j  = [N:-1:0]';
X1 = cos(j*pi/N);
U1 = u(X1);

a  = dct1(U1);
b  = zeros(N+2,1);

for j = [N-1:-1:1]
    b(j+1) = b(j+3) + 2*(j+1)*a(j+2);
end;
b(1) = a(2) + 0.5*b(3);

b  = -b(1:N+1,1);
b(1) = b(1)*2;

Ud1 = dct1(b)*N;
Ud1(2:N,1) = Ud1(2:N,1)/2;

Using Finite differences :

dX = 2/N;
X2 = [-1:dX:1]';
U2 = u(X2);

A  = -diag(ones(N,1),-1)+diag(ones(N,1),1);
A(1,1:3) = [-3 4 -1];
A(N+1,N-1:N+1) = [1 -4 3];
A  = A/(2*dX);

Ud2 = A*U2;


x = [-1:0.01:1];
plot(X1,Ud1,'ks',X2,Ud2,'k*',x,du(x),'k-');
legend(['Chebyshev, N = ' int2str(N)],...
       ['Central FD, N = ' int2str(N)],...
        'Exact');
xlabel('x'); ylabel('u`');