📄 seven_point_centered_uni_d1.m
字号:
function [D]=seven_point_centered_uni_D1(z0,zL,n)
%...
%... A. Vande Wouwer, P. Saucez and W.E. Schiesser (2002)
%...
%... function seven_point_centered_uni_D1 returns the differentiation matrix for
%... computing the first derivative, xz, of a variable x over the spatial domain
%... z0 < z < zL from classical seven-point, sixth-order finite difference
%... approximations (this function replaces dss006)
%...
%...
%... argument list
%...
%... z0 lower boundary value of z (input)
%...
%... zL upper boundary value of z (input)
%...
%... n number of grid points in the z domain including the
%... boundary points (input)
%...
%... origin of the approximation
%...
%... the mathematical details of the following taylor series (or
%... polynomials) are given in functions cdm1d2 and cdm1d4.
%...
%... seven-point formulas
%...
%... (1) left end, point i = 1
%...
%... 2 3 4
%... a(x2 = x1 + x1 ( dz) + x1 ( dz) + x1 ( dz) + x1 ( dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 ( dz) + x1 ( dz) + x1 ( dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... 2 3 4
%... b(x3 = x1 + x1 (2dz) + x1 (2dz) + x1 (2dz) + x1 (2dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 (2dz) + x1 (2dz) + x1 (2dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... 2 3 4
%... c(x4 = x1 + x1 (3dz) + x1 (3dz) + x1 (3dz) + x1 (3dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 (3dz) + x1 (3dz) + x1 (3dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... 2 3 4
%... d(x5 = x1 + x1 (4dz) + x1 (4dz) + x1 (4dz) + x1 (4dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 (4dz) + x1 (4dz) + x1 (4dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... 2 3 4
%... e(x6 = x1 + x1 (5dz) + x1 (5dz) + x1 (5dz) + x1 (5dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 (5dz) + x1 (5dz) + x1 (5dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... 2 3 4
%... f(x7 = x1 + x1 (6dz) + x1 (6dz) + x1 (6dz) + x1 (6dz)
%... z 1f 2z 2f 3z 3f 4z 4f
%...
%... 5 6 7
%... + x1 (6dz) + x1 (6dz) + x1 (6dz) + ...)
%... 5z 5f 6z 6f 7z 7f
%...
%... constants a, b, c, d, e and f are selected so that the coeffi-
%... cients of the x1 terms sum to one and the coefficients of
%... z
%... the x1 , x1 , x1 , x1 and x1 terms sum to zero.
%... 2z 3z 4z 5z 6z
%...
%... 1 1 1 1 1
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 1
%...
%... 2 2 2 2 2
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 0
%...
%... 3 3 3 3 3
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 0
%...
%... 4 4 4 4 4
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 0
%...
%... 5 5 5 5 5
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 0
%...
%... 6 6 6 6 6
%... a + 2 b + 3 c + 4 d + 5 e + 6 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for x1 gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... x1 = (1/6f)(-1764*x1 + 4320x2 - 5400x3 + 4800x4
%... z 6 (1)
%... + 2700x5 + 864x6 - 120x7) + o(dz )
%...
%... (2) interior point, i = 2
%...
%... the preceding taylor series can be summarized for point i = 2
%... as
%...
%... a(x1 = x2 + ...)
%...
%... b(x3 = x2 + ...)
%...
%... c(x4 = x2 + ...)
%...
%... d(x5 = x2 + ...)
%...
%... e(x6 = x2 + ...)
%...
%... f(x7 = x2 + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 1
%...
%... 2 2 2 2 2 2
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 0
%...
%... 3 3 3 3 3 3
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 0
%...
%... 4 4 4 4 4 4
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 0
%...
%... 5 5 5 5 5 5
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 0
%...
%... 6 6 6 6 6 6
%... -1 a + 1 b + 2 c + 3 d + 4 e + 5 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for x2 gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... x2 = (1/6f)(-120x1 - 924x2 + 1800x3 - 1200x4
%... z 6 (2)
%... + 600x5 - 180x6 + 24x7) + o(dz )
%...
%... (3) interior point, i = 3
%...
%... the preceding taylor series can be summarized for point i = 3
%... as
%...
%... a(x1 = x3 + ...)
%...
%... b(x2 = x3 + ...)
%...
%... c(x4 = x3 + ...)
%...
%... d(x5 = x3 + ...)
%...
%... e(x6 = x3 + ...)
%...
%... f(x7 = x3 + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 1
%...
%... 2 2 2 2 2 2
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 0
%...
%... 3 3 3 3 3 3
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 0
%...
%... 4 4 4 4 4 4
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 0
%...
%... 5 5 5 5 5 5
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 0
%...
%... 6 6 6 6 6 6
%... -2 a + -1 b + 1 c + 2 d + 3 e + 4 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for x3 gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... x3 = (1/6f)(24x1 - 288x2 - 420x3 + 960x4
%... z 6 (3)
%... - 360x5 + 96x6 - 12x7) + o(dz )
%...
%... (4) interior point, i ne 2, 3, n-2, n-1
%...
%... the preceding taylor series can be summarized for point i = i
%... as
%...
%... a(xi-3 = xi + ...)
%...
%... b(xi-2 = xi + ...)
%...
%... c(xi-1 = xi + ...)
%...
%... d(xi+1 = xi + ...)
%...
%... e(xi+2 = xi + ...)
%...
%... f(xi+3 = xi + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 1
%...
%... 2 2 2 2 2 2
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 0
%...
%... 3 3 3 3 3 3
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 0
%...
%... 4 4 4 4 4 4
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 0
%...
%... 5 5 5 5 5 5
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 0
%...
%... 6 6 6 6 6 6
%... -3 a + -2 b + -1 c + 1 d + 2 e + 3 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for xi gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... xi = (1/6f)(-12xi-3 + 108xi-2 - 540xi-1 + 0xi
%... z 6 (4)
%... + 540xi+1 - 108xi+2 + 12xi+3) + o(dz )
%...
%... (5) interior point, i = n-2
%...
%... the preceding taylor series can be summarized for point i = n-2
%... as
%...
%... a(xn-6 = xn-2 + ...)
%...
%... b(xn-5 = xn-2 + ...)
%...
%... c(xn-4 = xn-2 + ...)
%...
%... d(xn-3 = xn-2 + ...)
%...
%... e(xn-1 = xn-2 + ...)
%...
%... f(xn = xn-2 + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 1
%...
%... 2 2 2 2 2 2
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 0
%...
%... 3 3 3 3 3 3
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 0
%...
%... 4 4 4 4 4 4
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 0
%...
%... 5 5 5 5 5 5
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 0
%...
%... 6 6 6 6 6 6
%... -4 a + -3 b + -2 c + -1 d + 1 e + 2 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for xn-2 gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... xn-2 = (1/6f)(12xn-6 - 96xn-5 + 360xn-4 - 960xn-3
%... z 6 (5)
%... + 420xn-2 + 288xn-1 - 24xn) + o(dz )
%...
%... (6) interior point, i = n-1
%...
%... the preceding taylor series can be summarized for point i = n-1
%... as
%...
%... a(xn-6 = xn-1 + ...)
%...
%... b(xn-5 = xn-1 + ...)
%...
%... c(xn-4 = xn-1 + ...)
%...
%... d(xn-3 = xn-1 + ...)
%...
%... e(xn-2 = xn-1 + ...)
%...
%... f(xn = xn-1 + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 1
%...
%... 2 2 2 2 2 2
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 0
%...
%... 3 3 3 3 3 3
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 0
%...
%... 4 4 4 4 4 4
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 0
%...
%... 5 5 5 5 5 5
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 0
%...
%... 6 6 6 6 6 6
%... -5 a + -4 b + -3 c + -2 d + -1 e + 1 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for xn-1 gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... xn-1 = (1/6f)(-24xn-6 + 180xn-5 - 600xn-4 + 1200xn-3
%... z 6 (6)
%... - 1800xn-2 + 924xn-1 + 120xn) + o(dz )
%...
%... (7) right end, point i = n
%...
%... the preceding taylor series can be summarized for point i = n
%... as
%...
%... a(xn-6 = xn + ...)
%...
%... b(xn-5 = xn + ...)
%...
%... c(xn-4 = xn + ...)
%...
%... d(xn-3 = xn + ...)
%...
%... e(xn-2 = xn + ...)
%...
%... f(xn-1 = xn + ...)
%...
%... the corresponding algebraic equations are
%...
%... 1 1 1 1 1 1
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 1
%...
%... 2 2 2 2 2 2
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 0
%...
%... 3 3 3 3 3 3
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 0
%...
%... 4 4 4 4 4 4
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 0
%...
%... 5 5 5 5 5 5
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 0
%...
%... 6 6 6 6 6 6
%... -6 a + -5 b + -4 c + -3 d + -2 e + -1 f = 0
%...
%... simultaneous solution for a, b, c, d, e and f followed by the
%... solution of the preceding taylor series, truncated after the
%... x terms, for xn gives the following seven-point approxi-
%... 6z z
%... tion
%...
%... xn = (1/6f)(120xn-6 - 864xn-5 + 2700xn-4 - 4800xn-3
%... z 6 (7)
%... + 5400xn-2 - 4320xn-1 + 1764xn) + o(dz )
%...
%... the weighting coefficients for equations (1) to (7) can be
%... summarized as
%...
%... -1764 4320 -5400 4800 -2700 864 -120
%...
%... -120 -924 1800 -1200 600 -180 24
%...
%... 24 -288 -420 960 -360 96 -12
%...
%... 1/6f -12 108 -540 0 540 -108 12
%...
%... 12 -96 360 -960 420 288 -24
%...
%... -24 180 -600 1200 -1800 924 120
%...
%... 120 -864 2700 -4800 5400 -4320 1764
%...
%... which are the coefficients reported by bickley for n = 6, m =
%... 1, p = 0 to 6 (bickley, w. g., formulae for numerical differ-
%... entiation, math. gaz., vol. 25, 1941).
%...
%... equations (1) to (7) can now be programmed to generate the
%... derivative x (z) of function x(z)
%... z
%...
%... compute the spatial increment
dz=(zL-z0)/(n-1);
r6fdz=1./(720.*dz);
%...
%... sparse discretization matrix
%...
%... interior points
D=diag(-12*ones(n-3,1),-3)+diag(+108*ones(n-2,1),-2)+diag(-540*ones(n-1,1),-1)+...
+diag(+540*ones(n-1,1),+1)+diag(-108*ones(n-2,1),+2)+diag(+12*ones(n-3,1),+3);
%...
%... boundary points
D([1 2 3],1:7) = [-1764 +4320 -5400 +4800 -2700 +864 -120; -120 -924 +1800 -1200 +600 -180 +24; +24 -288 -420 +960 -360 +96 -12];
%...
D([n-2 n-1 n],(n-6):n) = [+12 -96 +360 -960 +420 +288 -24; -24 +180 -600 +1200 -1800 +924 +120; +120 -864 +2700 -4800 +5400 -4320 +1764];
%...
D=r6fdz*D;
D=sparse(D);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -