function [x,DM]=fourdif(N,m); % % [x, DM] = fourdif(N,m) computes the m'th derivative Fourier spectral % differentiation matrix on grid with N equispaced points in [0,2pi) % % Input: % N: Size of differentiation matrix. % M: Derivative required (non-negative integer) % % Output: % x: Equispaced points 0, 2pi/N, 4pi/N, ... , (N-1)2pi/N % DM: m'th order diffrentiation matrix % % % Explicit formulas are used to compute the matrices for m=1 and 2. % A discrete Fouier approach is employed for m>2. The program % computes the first column and first row and then uses the % toeplitz command to create the matrix. % For m=1 and 2 the code implements a "flipping trick" to % improve accuracy suggested by W. Don and A. Solomonoff in % SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994). % The flipping trick is necesary since sin t can be computed to high % relative precision when t is small whereas sin (pi-t) cannot. % % S.C. Reddy, J.A.C. Weideman 1998. x=2*pi*(0:N-1)'/N; % gridpoints h=2*pi/N; % grid spacing zi=sqrt(-1); kk=(1:N-1)'; n1=floor((N-1)/2); n2=ceil((N-1)/2); if m==0, % compute first column col1=[1; zeros(N-1,1)]; % of zeroth derivative row1=col1; % matrix, which is identity elseif m==1, % compute first column if rem(N,2)==0 % of 1st derivative matrix topc=cot((1:n2)'*h/2); col1=[0; 0.5*((-1).^kk).*[topc; -flipud(topc(1:n1))]]; else topc=csc((1:n2)'*h/2); col1=[0; 0.5*((-1).^kk).*[topc; flipud(topc(1:n1))]]; end; row1=-col1; % first row elseif m==2, % compute first column if rem(N,2)==0 % of 2nd derivative matrix topc=csc((1:n2)'*h/2).^2; col1=[-pi^2/3/h^2-1/6; -0.5*((-1).^kk).*[topc; flipud(topc(1:n1))]]; else topc=csc((1:n2)'*h/2).*cot((1:n2)'*h/2); col1=[-pi^2/3/h^2+1/12; -0.5*((-1).^kk).*[topc; -flipud(topc(1:n1))]]; end; row1=col1; % first row else % employ FFT to compute N1=floor((N-1)/2); % 1st column of matrix for m>2 N2=(-N/2)*(rem(m,2)==0)*ones(rem(N,2)==0); mwave=zi*[(0:N1) N2 (-N1:-1)]; col1=real(ifft((mwave.^m).*fft([1 zeros(1,N-1)]))); if rem(m,2)==0, row1=col1; % first row even derivative else col1=[0 col1(2:N)]'; row1=-col1; % first row odd derivative end; end; DM=toeplitz(col1,row1);