No BSD License
-
pade_firstder(F,pade,Isym)
Presentation
-
pade_init(nx,l0,lx,bd0,bd1)
Presentation
-
pade_secder(F,pade,Isym)
Presentation
-
pade_test(nperiod,b0,b1)
Presentation
-
View all files
from
6th order Pade Derivative
by Julien Reveillon
A set of functions to compute high order (6th) finite difference derivatives (first and second orde
|
| pade_secder(F,pade,Isym) |
%% 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 [D2F] = pade_secder(F,pade,Isym)
%% Head
%
% F : function to be derivated
% pade : structure initialized in pade_init.m
% Isym : symmetrical (=1) or anti-symmetrical (~=1)
% function (if necessary)
%
N = pade.nx;
DN2 = (1/pade.dx)^2;
NxM1 = N-1;
NxM2 = N-2;
NxM3 = N-3;
NxM4 = N-4;
NxM5 = N-5;
if(Isym ~= 1 )
Isym = 2;
end
ax(1:N) = pade.ax(1:N,Isym,2);
bx(1:N) = pade.bx(1:N,Isym,2);
cx(1:N) = pade.cx(1:N,Isym,2);
A1 = 5.5;
B1 = 4*(A1 -1)*DN2/3;
C1 = (10- A1)*DN2/3*0.25;
B2 = 12*DN2;
coef = ax(1)+ax(NxM2)+bx(1);
%% Derivative process : left boundary
D2F(3:NxM2) = ( (F(4:NxM1)-F(3:NxM2)) + ...
(F(2:NxM3)-F(3:NxM2)))* B1 + ...
( (F(5:N )-F(3:NxM2)) + ...
(F(1:NxM4)-F(3:NxM2)))* C1;
if ( pade.bd0 == 0 ) %periodic
D2F(1) = ( (F(2) -F(1)) +...
(F(NxM1) -F(1))) * B1 +...
( (F(3 ) -F(1)) +...
(F(NxM2) -F(1))) * C1;
D2F(2) = ( (F(3) -F(2)) +...
(F(1) -F(2))) * B1 +...
( (F(4) -F(2)) +...
(F(NxM1) -F(2))) * C1;
elseif ( pade.bd0 == 1 ) %bounded
D2F(1) = ( (F(1) -F(2))*27 + ...
(F(3) -F(1))*15 - ...
(F(4) -F(1)))*DN2;
D2F(2) = ( (F(3) -F(2)) + ...
(F(1) -F(2)))*B2;
elseif( pade.bd0 == 2 ) %symmetrical centered
if( Isym == 2 )
D2F(1) = ( (F(2) -F(1)) + ...
(-F(1) -F(1))) *B1 + ...
( (F(3) -F(1)) + ...
(-F(2) -F(1))) *C1;
D2F(2) = ( (F(3) -F(2)) + ...
(F(1) -F(2))) *B1 + ...
( (F(4) -F(2)) + ...
(-F(1) -F(2))) *C1;
else
D2F(1) = (F(2) -F(1)) *B1 + ...
( (F(3) -F(1)) + ...
(F(2) -F(1))) *C1;
D2F(2) = ( (F(3) -F(2)) + ...
(F(1) -F(2))) *B1 + ...
( (F(4) -F(2)) + ...
(F(1) -F(2))) *C1;
end
elseif( pade.bd0 == 3 ) %symmetrical non centered
if( Isym == 2 )
D2F(1) = -2*F(1)*B1 -2*F(1)*C1;
D2F(2) = ( (F(3) -F(2)) + ...
(F(1) -F(2))) *B1 + ...
(F(4) -3*F(2))*C1;
else
D2F(1) = (2*F(2) -2*F(1))*B1 + ...
(2*F(3) -2*F(1))*C1;
D2F(2) = ( (F(3) -F(2)) + ...
(F(1) -F(2))) *B1 + ...
( (F(4) -F(2))) *C1;
end
else
fprintf(1, 'error SecondDerivative \n');
pause;
end
%% Derivative process : right boundary
if ( pade.bd1 == 0 )
D2F(NxM1) = ( (F(1) -F(NxM1)) +...
(F(NxM2) -F(NxM1))) * B1 +...
( (F(2) -F(NxM1)) +...
(F(NxM3) -F(NxM1))) * C1;
elseif ( pade.bd1 == 1 );
D2F(N) = ( (F(N) -F(NxM1))*27 + ...
(F(NxM2) -F(N))*15 - ...
(F(NxM3) -F(N)))*DN2;
D2F(NxM1) = ( (F(N) -F(NxM1)) + ...
(F(NxM2) -F(NxM1)))*B2;
elseif( pade.bd1 == 2 )
if( Isym == 2 )
D2F(N) = ( (F(NxM1) -F(N)) + ...
(-F(N) -F(N))) *B1 + ...
( (F(NxM2) -F(N)) + ...
(-F(NxM1) -F(N))) *C1;
D2F(NxM1) = ( (F(NxM2) -F(NxM1)) + ...
(F(N) -F(NxM1))) *B1 +...
( (F(NxM3) -F(NxM1)) +...
(-F(N) -F(NxM1))) *C1;
else
D2F(N) = (F(NxM1) -F(N)) *B1 + ...
( (F(NxM2) -F(N)) + ...
(F(NxM1) -F(N))) *C1;
D2F(NxM1) = ( (F(NxM2) -F(NxM1)) + ...
(F(N) -F(NxM1))) *B1 +...
( (F(NxM3) -F(NxM1)) +...
(F(N) -F(NxM1))) *C1;
end
elseif( pade.bd1 == 3 );
if( Isym == 2 )
D2F(N) = -2*F(N)*B1-2*F(N)*C1;
D2F(NxM1) = ( (F(NxM2) -F(NxM1)) + ...
(F(N) -F(NxM1))) *B1 + ...
(F(NxM3) -3*F(NxM1)) *C1 ;
else
D2F(N) = (2*F(NxM1) -2*F(N))*B1 + ...
(2*F(NxM2) -2*F(N))*C1;
D2F(NxM1) = ( (F(NxM2) -F(NxM1)) + ...
(F(N) -F(NxM1))) *B1 + ...
( (F(NxM3) -F(NxM1))) *C1;
end
else
fprintf(1, 'error SecondDerivative \n');
pause;
end
%% Thomas Algorithm phase II
if ( pade.bd0 == 0 )
% periodic
D2F(1) = D2F(1)*cx(1);
for I = 2:NxM2
D2F(I) = (D2F(I)-D2F(I-1))*cx(I);
end
for I = NxM3:-1:1
D2F(I) = D2F(I)-D2F(I+1)*cx(I);
end
%
D2F(NxM1) = (D2F(NxM1)-D2F(1)-D2F(NxM2))/coef;
%
for I = 1:NxM2
D2F(I) = D2F(I)+D2F(NxM1)*ax(I);
end
%
D2F(N) = D2F(1);
else
% non-periodic
D2F(1) = D2F(1)/bx(1);
for I = 2:NxM1
D2F(I) = (D2F(I)-D2F(I-1))/bx(I);
end
%
D2F(N) = (D2F(N)-D2F(NxM1)*ax(N))/bx(N);
%
for I = NxM1:-1:1
D2F(I) = D2F(I)-D2F(I+1)*cx(I);
end
end
end
|
|
Contact us at files@mathworks.com
