sandbox/easystab/david/dif1D.m
- Function dif1D
- Finite differences (optimized case)
- Finite differences (basic implementation, equivalent to the previous)
- Finite differences on a periodic domain
- Chebyshev collocation
- Fourier collocation
- Hermite collocation on infinite domain.
- Stretched Chebyshev ; infinite ; exponential
- Stretched Chebyshev ; infinite ; algebraic
- Stretched Chebyshev ; semi-infinite ; exponential
- Stretched Chebyshev ; semi-infinite ; algebraic
Function dif1D
function [d,dd,w,s]=dif1D(type,s0,L,N,opt);
This function creates a mesh s, differentiation matrices of first and second order d and dd, as well as the weight w.
On output: s is a column-vector containing the mesh points, starting from the point at the left. d and dd are the derivation matrices (quare matrices, in full or sparse storing format) w is a line-vector containing the weights.
Last parameter ‘opt’ is optional. For Finite-Difference it is interpreted as the number of points of the stencil (default :3) and for Stretched Chebyschev it is interpreted as the parameter b controlling the stretching.
switch type
Finite differences (optimized case)
For a domain [s0,s0+L]
Last parameter ‘pts’ is optional, if absent it is set to 3 (relevant value in almost all cases)
case 'fd'
if(nargin==4)
pts = 3;
else
pts = opt;
end
scale=L/2;
[s,d] = fddif(N,1,pts);
[s,dd] = fddif(N,2,pts);
s=s*scale;
d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
w=([diff(s'),0]+[0,diff(s')])/2;
Finite differences (basic implementation, equivalent to the previous)
case 'fds'
dx = L/(N-1); % grid spacing
s = [ s0 : L/(N-1) : s0+L ]'; % creates an equispaced mesh with N point on the interval
for i=2:N-1 % centered three-point formulas
d(i,i) = 0;
d(i,i+1) = 1/(2*dx);
d(i,i-1) = -1/(2*dx);
dd(i,i) = -2/dx^2;
dd(i,i+1) = 1/(dx^2);
dd(i,i-1) = 1/(dx^2);
end
% uncentred formulas for first and last grid points
d(1,1) = -3/(2*dx); d(1,2) = 4/(2*dx) ; d(1,3) = -1/(2*dx);
d(N,N-2) = 1/(2*dx); d(N,N-1) = -4/(2*dx); d(N,N) = 3/(2*dx);
dd(1,1) = 2/dx^2; dd(1,2) = -5/dx^2; dd(1,3) = 4/dx^2; dd(1,4) = -1/dx^2;
dd(N,N) = 2/dx^2; dd(N,N-1) = -5/dx^2; dd(N,N-2) = 4/dx^2; dd(N,N-3) = -1/dx^2;
% "weight" to compute integrals using the trapeze rule
for i=2:N-1
w(i) = dx;
end
w(1) = dx/2; w(N) = dx/2;
Finite differences on a periodic domain
Domain is defined as [s0,s0+L]
case 'fp' % periodic finite differences
scale=L/2;
[s,d] = fdper(N,1,opt);
[s,dd] = fdper(N,2,opt);
s=s*scale;
d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
w=ones(1,N)*(s(2)-s(1));
Chebyshev collocation
on domain [s0,s0+L]
case 'cheb' % chebychev
scale=-L/2;
[s,DM] = chebdif(N,2);
d=DM(:,:,1);
dd=DM(:,:,2);
s=s*scale;
d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
w=L*clencurt(N)/2;
Fourier collocation
(requires a periodic domain [s0,s0+L])
case 'fou' % Fourier
scale=L/(2*pi);
[s, d] = fourdif(N, 1);
[s, dd] = fourdif(N, 2);
s=s*scale;
w=ones(1,N)*(s(2)-s(1));
d=d/scale; dd=dd/scale^2; s=s-s(1)+s0;
Hermite collocation on infinite domain.
Domain is centred on s0, L controls the density
(increase L to cluster the points towards s0)
case 'her' % hermite
scale=1;
[s,DM] = herdif(N,2,L);
d=DM(:,:,1);
dd=DM(:,:,2);
s=s*scale;
d=d/scale; dd=dd/scale^2; s=s+s0;
w=zeros(1,N); % weight is not implemented yet !!!
Stretched Chebyshev ; infinite ; exponential
- s0 is the center of the domain.
- L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
- optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
case 'chebInfExp' % chebychev with coordinate stretching to treat large domains
if(nargin==4)
b = 0.999;
else
b = opt;
end
[dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
s = L*atanh(xi)+s0; % coordinate stretching (exponential)
G(:, 1) = - (xi.^2 - 1) /L; % G1 = ds / dxi
G(:, 2) = 2 *xi .* (xi.^2 - 1) /L^2; % G2 = d^2 s /d xi^2
G(:, 3) = - 2 * (3 * xi.^2 - 1) .* (xi.^2 - 1) / L^3;
G(:, 4) = 8 * xi .* (3 * xi.^2 - 2) .* (xi.^2 - 1) / L^4;
d = G(:,1).*dxi;
dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
w = wxi./G(:,1)';
Stretched Chebyshev ; infinite ; algebraic
- s0 is the center of the domain.
- L controls the clustering of points around the origin : about half points will be localised in [s0-L,s0+L]
- optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
case 'chebInfAlg' % chebychev with coordinate stretching to treat large domains
if(nargin==4)
b = 0.999;
else
b = opt;
end
[dxi,ddxi,wxi,xi]=dif1D('cheb',-b,2*b,N); % Standart Chebyshev for xi \in [-b,b]
s = L * xi ./ sqrt(1 - xi.^2);
K = sqrt(s.^2 + L^2);
G(:, 1) = L^2 ./ K.^3;
G(:, 2) = - 3 * s * L^2 ./ K.^5;
G(:, 3) = 3 * (4 * s.^2 - L^2) * L^2 ./ K.^7;
G(:, 4) = - 15 * s .* (4 * s.^2 - 3 * L^2) * L^2 ./ K.^9;
d = G(:,1).*dxi;
dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
w = wxi./G(:,1)';
Stretched Chebyshev ; semi-infinite ; exponential
- s0 is the left bound of the domain.
- L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
- optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
case 'chebSemiInfExp' % chebychev with coordinate stretching to treat large domains
if(nargin==4)
b = 0.999;
else
b = opt;
end
[dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,b]
s = L*atanh(xi)+s0; % coordinate stretching (exponential)
G(:, 1) = - (xi.^2 - 1) /L; % G1 = ds / dxi
G(:, 2) = 2 * xi .* (xi.^2 - 1) /L^2; % G2 = d^2 s /d xi^2
d = G(:,1).*dxi;
dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
w = wxi./G(:,1)';
Stretched Chebyshev ; semi-infinite ; algebraic
- s0 is the left bound of the domain.
- L controls the clustering of points around the origin : about half points will be localised in [s0,s0+L]
- optional parameter b (default value 0.999) controls the stretching (b should be smaller than 1).
case 'chebSemiInfAlg' % chebychev with coordinate stretching to treat large domains
if(nargin==4)
b = 0.999;
else
b = opt;
end
[dxi,ddxi,wxi,xi]=dif1D('cheb',0,b,N); % Standart Chebyshev for xi \in [0,1]
s = L * xi ./ sqrt(1 - xi.^2);
K = sqrt(s.^2 + L^2);
G(:, 1) = L^2 ./ K.^3;
G(:, 2) = - 3 * s .* L^2 ./ K.^5;
d = G(:,1).*dxi;
dd = G(:,1).^2.*ddxi+G(:,2).*dxi;
w = wxi./G(:,1)';
end