%% Presentation
%
% Pade : 6th order finite difference derivative
%
% .m file in the package:
% pade_init.m : initialization of the Pade coefficients
% pade_firstder : first order derivative
% pade_secder : second order derivative
%
%% Author
% Julien Reveillon -- Julien.Reveillon(at)coria.fr
% CORIA, University of Rouen, France
%
%% History
% January 27, 2007 : Version 1.0 = 1D, regular grid
%
%% Remarks
%
% 1) A sixth order FD derivative can not be used with too strong gradients
% otherwise it leads to oscilloations of the derivatives
%
% 2) The grid (xp) is generated inside the pade_init function.
% A slight modification would authorize external grids
% but beware to the boundary conditions
function [ax,bx,cx,xp,dx] = pade_init(nx,l0,lx,bd0,bd1)
%% Head
% Initilization of the Pade coefficients
% pade structure:
% nx : number of grid points
% l0 : first grid point
% lx : length of the grid
% bd0 : left boundary condition
% bd1 : right boundary condition
% =0 : periodical
% =1 : bounded
% =2 : centered symmetry
% =3 : non-centered symmetry
%% Data Initialization
if ((bd0 == 0) && (bd1 ~= 0)) || ((bd1 == 0) && (bd0 ~= 0))
fprintf(1,'Periodic boundary error\n');
pause;
end
% Preallocation of the xp grid and
% the tridigonal matrix coefficients : ax, bx, cx
% ax(Number of node, 1 symetric/ 2 antisym, derivative 1 or 2)
ax(1:nx,1:2,1:2) = 0;
bx(1:nx,1:2,1:2) = 0;
cx(1:nx,1:2,1:2) = 0;
xp(1:nx) = 0;
%% Grid generation
if ((bd0 ~= 2) && (bd1 ~= 2))
dx = lx/(nx-1);
for I = 1:nx
xp(I) = l0+(I-1)*dx;
end
elseif ((bd0 == 2) && (bd1 ~= 2))
dx = lx/(nx-0.5);
for I = 1:nx
xp(I) = l0+(I-0.5)*dx;
end
elseif ((bd0 ~= 2) && (bd1 == 2))
dx = lx/(nx-0.5);
for I = 1:nx
xp(I) = l0+(I-1)*dx;
end
elseif ((bd0 == 2) && (bd1 == 2))
dx = lx/nx;
for I = 1:nx
xp(I) = l0+(I-0.5)*dx;
end
end
%% Matrix coefficient initialization
%
% First derivative generic initialization
%
A1 = 3;
A2 = 16*(2*A1 + 1)/(40 - A1);
ax(1:nx,1:2,1) = 1;
bx(1:nx,1:2,1) = A1;
cx(1:nx,1:2,1) = 1;
%
% Second derivative generic initialisation
%
SECA1 = 5.5;
ax(1:nx,1:2,2) = 1;
bx(1:nx,1:2,2) = SECA1;
cx(1:nx,1:2,2) = 1;
%
% Modifications according to the boundary conditions
%
if ( bd0 == 1)
%first derivative
bx(1 ,1:2,1) = 2;
cx(1 ,1:2,1) = 4;
bx(2 ,1:2,1) = 4;
bx(3 ,1:2,1) = A2;
%second derivative
cx(1 ,1:2,2) = 11;
bx(1 ,1:2,2) = 1;
bx(2 ,1:2,2) = 10;
elseif ( bd0 == 2 )
%first derivative
bx(1 ,1,1) = 2;
bx(1 ,2,1) = 4;
%second derivative
bx(1 ,1,2) = 6.5;
bx(1 ,2,2) = 4.5;
elseif ( bd0 == 3 )
%first derivative
cx(1 ,1,1) = 0;
cx(1 ,2,1) = 2;
%second derivative
cx(1 ,1,2) = 2;
cx(1 ,2,2) = 0;
end
if ( bd1 == 1)
%first derivative
ax(nx ,1:2,1) = 4;
bx(nx ,1:2,1) = 2;
bx(nx-1,1:2,1) = 4;
bx(nx-2,1:2,1) = A2;
%second derivative
ax(nx ,1:2,2) = 11;
bx(nx ,1:2,2) = 1;
bx(nx-1,1:2,2) = 10;
elseif ( bd1 == 2 )
%first derivative
bx(nx ,1,1) = 2;
bx(nx ,2,1) = 4;
%second derivative
bx(nx ,1,2) = 6.5;
bx(nx ,2,2) = 4.5;
elseif ( bd1 == 3 )
%first derivative
ax(nx ,1,1) = 0;
ax(nx ,2,1) = 2;
%second derivative
ax(nx ,1,2) = 2;
ax(nx ,2,2) = 0;
end
%% Thomas Algorithm phase I
if(bd0 == 0)
cx(nx-2,1:2,1:2) = 1;
cx(1,1:2,1:2) = 1/bx(1,1:2,1:2);
for I = 2:nx-2
bx(I,1:2,1:2) = bx(I,1:2,1:2)-cx(I-1,1:2,1:2);
cx(I,1:2,1:2) = 1/bx(I,1:2,1:2);
end
%
ax(1,1:2,1:2) = -1/bx(1,1:2,1:2);
for I = 2:nx-3
ax(I,1:2,1:2) = -ax(I-1,1:2,1:2)./bx(I,1:2,1:2);
end
ax(nx-2,1:2,1:2) = (-1-ax(nx-3,1:2,1:2))./bx(nx-2,1:2,1:2);
%
for I = nx-3:-1:1
ax(I,1:2,1:2) = ax(I,1:2,1:2)-ax(I+1,1:2,1:2).*cx(I,1:2,1:2);
end
else
cx(nx,1:2,1:2) = 1;
cx(1,1:2,1:2) = cx(1,1:2,1:2)./bx(1,1:2,1:2);
for I = 2:nx
bx(I,1:2,1:2) = bx(I,1:2,1:2)-cx(I-1,1:2,1:2).*ax(I,1:2,1:2);
cx(I,1:2,1:2) = cx(I,1:2,1:2)./bx(I,1:2,1:2);
end
end
%
%End of Subroutine
%
fprintf(1, 'Initialization OK \n');
end
