Mercurial > repos > public > sbplib
changeset 52:4f5a65f49035
Fixed the upwind operators.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 11 Nov 2015 15:43:15 -0800 |
| parents | 0be702829bb9 |
| children | cca09cc5121d |
| files | +sbp/Upwind.m +sbp/upwind3_3bp.m +sbp/upwind4.m +sbp/upwind6.m +sbp/upwind8.m |
| diffstat | 5 files changed, 84 insertions(+), 51 deletions(-) [+] |
line wrap: on
line diff
--- a/+sbp/Upwind.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/Upwind.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,4 +1,4 @@ -classdef Ordinary < sbp.OpSet +classdef Upwind < sbp.OpSet properties norms % Struct containing norm matrices such as H,Q, M boundary % Struct contanging vectors for boundry point approximations @@ -9,10 +9,10 @@ end methods - function obj = Ordinary(m,h,order) + function obj = Upwind(m,h,order) if order == 3 - [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind3(m,h); + [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind3_3bp(m,h); elseif order == 4 [H, HI, Dp, Dm, e_1, e_m] = sbp.upwind4(m,h); elseif order == 6 @@ -28,7 +28,6 @@ obj.norms.H = H; obj.norms.HI = HI; - obj.norms.Q = Q; obj.boundary.e_1 = e_1; obj.boundary.e_m = e_m;
--- a/+sbp/upwind3_3bp.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/upwind3_3bp.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,22 +1,29 @@ function [H, HI, Dp, Dm, e_1, e_m] = upwind3(m,h) - H(1:3,1:3)=[0.3e1 / 0.8e1 0 0; 0 0.7e1 / 0.6e1 0; 0 0 0.23e2 / 0.24e2;]; - H(m-2:m,m-2:m)=fliplr(flipud(H(1:3,1:3))); - H=H*h; - H=spdiag(ones(m,1),0); - HI=inv(H); + Hv = ones(m,1); + Hv(1:3) = [3/8; 7/6; 23/24]; + Hv(m-2:m) = rot90(Hv(1:3),2); + Hv = Hv*h; + H = spdiag(Hv,0); + HI = spdiag(1./Hv,0); - Qp=(-1/3*spdiag(ones(m-1,1),-1)-1/2*spdiag(ones(m,1),0)+1*spdiag(ones(m-1,1),+1)-1/6*spdiag(ones(m-2,1),+2)); - %Q_U = [-0.49e2 / 0.2640e4 0.11167e5 / 0.15840e5 -0.1541e4 / 0.7920e4 0.43e2 / 0.5280e4; -0.9403e4 / 0.15840e5 -0.551e3 / 0.2640e4 0.1779e4 / 0.1760e4 -0.2311e4 / 0.7920e4; 0.659e3 / 0.7920e4 -0.751e3 / 0.1760e4 -0.1541e4 / 0.2640e4 0.21287e5 / 0.15840e5; 0.51e2 / 0.1760e4 -0.551e3 / 0.7920e4 -0.3683e4 / 0.15840e5 -0.713e3 / 0.880e3;]; - %Q_U = [-0.29e2 / 0.1680e4 0.7067e4 / 0.10080e5 -0.961e3 / 0.5040e4 0.23e2 / 0.3360e4; -0.6023e4 / 0.10080e5 -0.331e3 / 0.1680e4 0.1119e4 / 0.1120e4 -0.1451e4 / 0.5040e4; 0.439e3 / 0.5040e4 -0.491e3 / 0.1120e4 -0.961e3 / 0.1680e4 0.13507e5 / 0.10080e5; 0.31e2 / 0.1120e4 -0.331e3 / 0.5040e4 -0.2383e4 / 0.10080e5 -0.453e3 / 0.560e3;]; + q_diags = [-1, 0, 1, 2]; + q_stencil = [-1/3 -1/2 1 -1/6]; + Qp = stripeMatrix(q_stencil, q_diags,m); - Q_U = [-0.1e1 / 0.24e2 0.17e2 / 0.24e2 -0.1e1 / 0.6e1; -0.13e2 / 0.24e2 -0.1e1 / 0.4e1 0.23e2 / 0.24e2; 0.1e1 / 0.12e2 -0.11e2 / 0.24e2 -0.11e2 / 0.24e2;]; + Q_U = [ + -1/24 17/24 -1/6; + -13/24 -1/4 23/24; + 1/12 -11/24 -11/24; + ]; Qp(1:3,1:3)=Q_U; - Qp(m-2:m,m-2:m)=flipud( fliplr(Q_U(1:3,1:3) ) )'; %%% This is different from standard SBP + Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; - e_1=sparse(m,1);e_1(1)=1; - e_m=sparse(m,1);e_m(m)=1; + e_1=sparse(m,1); + e_1(1)=1; + e_m=sparse(m,1); + e_m(m)=1; Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
--- a/+sbp/upwind4.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/upwind4.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,25 +1,30 @@ function [H, HI, Dp, Dm, e_1, e_m] = upwind4(m,h) - H=spdiag(ones(m,1),0); - H(1:4,1:4)=spdiag([49/144 61/48 41/48 149/144]); - H(m-3:m,m-3:m)=fliplr(flipud(spdiag([49/144 61/48 41/48 149/144]))); - H=H*h; - HI=inv(H); + Hv=ones(m,1); + Hv(1:4)=[49/144 61/48 41/48 149/144]; + Hv(m-3:m)=rot90(Hv(1:4),2); + Hv = Hv*h; + H = spdiag(Hv,0); + HI = spdiag(1./Hv,0); + q_diags = [-1, 0, 1, 2 3]; + q_stencil = [-1/4 -5/6 3/2 -1/2 1/12]; + Qp = stripeMatrix(q_stencil, q_diags,m); - Qp=(-1/4*spdiag(ones(m-1,1),-1)-5/6*spdiag(ones(m,1),0)+3/2*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/12*spdiag(ones(m-3,1),+3)); - %Q_U = [-0.49e2 / 0.2640e4 0.11167e5 / 0.15840e5 -0.1541e4 / 0.7920e4 0.43e2 / 0.5280e4; -0.9403e4 / 0.15840e5 -0.551e3 / 0.2640e4 0.1779e4 / 0.1760e4 -0.2311e4 / 0.7920e4; 0.659e3 / 0.7920e4 -0.751e3 / 0.1760e4 -0.1541e4 / 0.2640e4 0.21287e5 / 0.15840e5; 0.51e2 / 0.1760e4 -0.551e3 / 0.7920e4 -0.3683e4 / 0.15840e5 -0.713e3 / 0.880e3;]; - %Q_U = [-0.29e2 / 0.1680e4 0.7067e4 / 0.10080e5 -0.961e3 / 0.5040e4 0.23e2 / 0.3360e4; -0.6023e4 / 0.10080e5 -0.331e3 / 0.1680e4 0.1119e4 / 0.1120e4 -0.1451e4 / 0.5040e4; 0.439e3 / 0.5040e4 -0.491e3 / 0.1120e4 -0.961e3 / 0.1680e4 0.13507e5 / 0.10080e5; 0.31e2 / 0.1120e4 -0.331e3 / 0.5040e4 -0.2383e4 / 0.10080e5 -0.453e3 / 0.560e3;]; - - Q_U = [-0.1e1 / 0.48e2 0.205e3 / 0.288e3 -0.29e2 / 0.144e3 0.1e1 / 0.96e2; -0.169e3 / 0.288e3 -0.11e2 / 0.48e2 0.33e2 / 0.32e2 -0.43e2 / 0.144e3; 0.11e2 / 0.144e3 -0.13e2 / 0.32e2 -0.29e2 / 0.48e2 0.389e3 / 0.288e3; 0.1e1 / 0.32e2 -0.11e2 / 0.144e3 -0.65e2 / 0.288e3 -0.13e2 / 0.16e2;]; + Q_U = [ + -0.1e1/0.48e2 0.205e3/0.288e3 -0.29e2/0.144e3 0.1e1/0.96e2; + -0.169e3/0.288e3 -0.11e2/0.48e2 0.33e2/0.32e2 -0.43e2/0.144e3; + 0.11e2/0.144e3 -0.13e2/0.32e2 -0.29e2/0.48e2 0.389e3/0.288e3; + 0.1e1/0.32e2 -0.11e2/0.144e3 -0.65e2/0.288e3 -0.13e2/0.16e2; + ]; Qp(1:4,1:4)=Q_U; - Qp(m-3:m,m-3:m)=flipud( fliplr(Q_U(1:4,1:4) ) )'; %%% This is different from standard SBP + Qp(m-3:m,m-3:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; - e_1=zeros(m,1);e_1(1)=1; - e_m=zeros(m,1);e_m(m)=1; + e_1=sparse(m,1);e_1(1)=1; + e_m=sparse(m,1);e_m(m)=1; Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
--- a/+sbp/upwind6.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/upwind6.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,22 +1,31 @@ function [H, HI, Dp, Dm, e_1, e_m] = upwind6(m,h) - H=spdiag(ones(m,1),0); - H(1:6,1:6)=[0.13613e5 / 0.43200e5 0 0 0 0 0; 0 0.12049e5 / 0.8640e4 0 0 0 0; 0 0 0.535e3 / 0.864e3 0 0 0; 0 0 0 0.1079e4 / 0.864e3 0 0; 0 0 0 0 0.7841e4 / 0.8640e4 0; 0 0 0 0 0 0.43837e5 / 0.43200e5;]; - H(m-5:m,m-5:m)=fliplr(flipud(H(1:6,1:6))); - H=H*h; - HI=inv(H); + Hv = ones(m,1); + Hv(1:6)=[0.13613e5/0.43200e5; 0.12049e5/0.8640e4 ; 0.535e3/0.864e3 ; 0.1079e4/0.864e3 ; 0.7841e4/0.8640e4 ; 0.43837e5/0.43200e5]; + Hv(m-5:m)=rot90(Hv(1:6),2); + Hv = Hv*h; + H = spdiag(Hv,0); + HI = spdiag(1./Hv,0); + q_diags = [-2 -1 0 1 2 3 4]; + q_stencil = [1/30 -2/5 -7/12 4/3 -1/2 2/15 -1/60]; + Qp = stripeMatrix(q_stencil, q_diags,m); - Qp=(1/30*spdiag(ones(m-2,1),-2)-2/5*spdiag(ones(m-1,1),-1)-7/12*spdiag(ones(m,1),0)+4/3*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+2/15*spdiag(ones(m-3,1),+3)-1/60*spdiag(ones(m-4,1),+4)); - %Q_U = [-0.375523e6 / 0.237369744e9 0.73393569793e11 / 0.111599208000e12 -0.1244182786433e13 / 0.14954293872000e14 -0.245555527093e12 / 0.2492382312000e13 0.131952309637e12 / 0.14954293872000e14 0.62907543023e11 / 0.3738573468000e13; -0.4818422888131e13 / 0.7477146936000e13 -0.10098911e8 / 0.307701520e9 0.796148045807e12 / 0.1495429387200e13 0.487507785187e12 / 0.2990858774400e13 0.8887918703e10 / 0.249238231200e12 -0.804641655323e12 / 0.14954293872000e14; 0.925771462433e12 / 0.14954293872000e14 -0.601394962127e12 / 0.1495429387200e13 -0.674340047e9 / 0.4153970520e10 0.560418114917e12 / 0.747714693600e12 -0.1043208352973e13 / 0.2990858774400e13 0.42207619661e11 / 0.356054616000e12; 0.285728443093e12 / 0.2492382312000e13 -0.873874460227e12 / 0.2990858774400e13 -0.34960363091e11 / 0.106816384800e12 -0.1528791703e10 / 0.4153970520e10 0.1741615892927e13 / 0.1495429387200e13 -0.6107723941553e13 / 0.14954293872000e14; -0.213784143637e12 / 0.14954293872000e14 0.1011411511e10 / 0.35605461600e11 0.141331778093e12 / 0.2990858774400e13 -0.632571739487e12 / 0.1495429387200e13 -0.1464298201e10 / 0.2769313680e10 0.9523869708211e13 / 0.7477146936000e13; -0.60750297023e11 / 0.3738573468000e13 0.614597810123e12 / 0.14954293872000e14 -0.16980260027e11 / 0.2492382312000e13 -0.343200536047e12 / 0.14954293872000e14 -0.389757912373e12 / 0.1068163848000e13 -0.964053659e9 / 0.1661588208e10;]; - Q_U =[-0.265e3 / 0.128688e6 0.1146190567e10 / 0.1737288000e10 -0.1596619e7 / 0.18384000e8 -0.55265831e8 / 0.579096000e9 0.26269819e8 / 0.3474576000e10 0.2464501e7 / 0.144774000e9; -0.1116490567e10 / 0.1737288000e10 -0.8839e4 / 0.214480e6 0.190538869e9 / 0.347457600e9 0.102705469e9 / 0.694915200e9 0.413741e6 / 0.9651600e7 -0.191689861e9 / 0.3474576000e10; 0.1096619e7 / 0.18384000e8 -0.135385429e9 / 0.347457600e9 -0.61067e5 / 0.321720e6 0.45137333e8 / 0.57909600e8 -0.253641811e9 / 0.694915200e9 0.70665929e8 / 0.579096000e9; 0.66965831e8 / 0.579096000e9 -0.208765789e9 / 0.694915200e9 -0.17623253e8 / 0.57909600e8 -0.18269e5 / 0.45960e5 0.410905829e9 / 0.347457600e9 -0.477953317e9 / 0.1158192000e10; -0.49219819e8 / 0.3474576000e10 0.293299e6 / 0.9651600e7 0.26422771e8 / 0.694915200e9 -0.141938309e9 / 0.347457600e9 -0.346583e6 / 0.643440e6 0.2217185207e10 / 0.1737288000e10; -0.2374501e7 / 0.144774000e9 0.142906261e9 / 0.3474576000e10 -0.3137129e7 / 0.579096000e9 -0.29884283e8 / 0.1158192000e10 -0.630168407e9 / 0.1737288000e10 -0.3559e4 / 0.6128e4;]; + Q_U =[ + -0.265e3/0.128688e6 0.1146190567e10/0.1737288000e10 -0.1596619e7/0.18384000e8 -0.55265831e8/0.579096000e9 0.26269819e8/0.3474576000e10 0.2464501e7/0.144774000e9; + -0.1116490567e10/0.1737288000e10 -0.8839e4/0.214480e6 0.190538869e9/0.347457600e9 0.102705469e9/0.694915200e9 0.413741e6/0.9651600e7 -0.191689861e9/0.3474576000e10; + 0.1096619e7/0.18384000e8 -0.135385429e9/0.347457600e9 -0.61067e5/0.321720e6 0.45137333e8/0.57909600e8 -0.253641811e9/0.694915200e9 0.70665929e8/0.579096000e9; + 0.66965831e8/0.579096000e9 -0.208765789e9/0.694915200e9 -0.17623253e8/0.57909600e8 -0.18269e5/0.45960e5 0.410905829e9/0.347457600e9 -0.477953317e9/0.1158192000e10; + -0.49219819e8/0.3474576000e10 0.293299e6/0.9651600e7 0.26422771e8/0.694915200e9 -0.141938309e9/0.347457600e9 -0.346583e6/0.643440e6 0.2217185207e10/0.1737288000e10; + -0.2374501e7/0.144774000e9 0.142906261e9/0.3474576000e10 -0.3137129e7/0.579096000e9 -0.29884283e8/0.1158192000e10 -0.630168407e9/0.1737288000e10 -0.3559e4/0.6128e4; + ]; Qp(1:6,1:6)=Q_U; - Qp(m-5:m,m-5:m)=flipud( fliplr(Q_U(1:6,1:6) ) )'; %%% This is different from standard SBP + Qp(m-5:m,m-5:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; - e_1=zeros(m,1);e_1(1)=1; - e_m=zeros(m,1);e_m(m)=1; + e_1=sparse(m,1);e_1(1)=1; + e_m=sparse(m,1);e_m(m)=1; Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;
--- a/+sbp/upwind8.m Wed Nov 11 15:15:26 2015 -0800 +++ b/+sbp/upwind8.m Wed Nov 11 15:43:15 2015 -0800 @@ -1,22 +1,35 @@ function [H, HI, Dp, Dm, e_1, e_m] = upwind8(m,h) - H=spdiag(ones(m,1),0); - H(1:8,1:8)=[0.7489399e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5537831e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.103373e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.261259e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.298231e6 / 0.725760e6 0 0 0; 0 0 0 0 0 0.515917e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3349159e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25639991e8 / 0.25401600e8;]; + Hv=ones(m,1); + Hv(1:8) = [0.7489399e7/0.25401600e8; 0.5537831e7/0.3628800e7; 0.103373e6/0.403200e6; 0.261259e6/0.145152e6; 0.298231e6/0.725760e6; 0.515917e6/0.403200e6; 0.3349159e7/0.3628800e7; 0.25639991e8/0.25401600e8]; + Hv(m-7:m)=rot90(Hv(1:8),2); + Hv=Hv*h; + + H = spdiag(Hv,0); + HI = spdiag(1./Hv,0); - H(m-7:m,m-7:m)=fliplr(flipud(H(1:8,1:8))); - H=H*h; - HI=inv(H); + q_diags = [-3 -2 -1 0 1 2 3 4 5]; + q_stencil = [-1/168 +1/14 -1/2 -9/20 +5/4 -1/2 +1/6 -1/28 +1/280]; + Qp = stripeMatrix(q_stencil, q_diags,m); + - Qp=(-1/168*spdiag(ones(m-3,1),-3)+1/14*spdiag(ones(m-2,1),-2)-1/2*spdiag(ones(m-1,1),-1)-9/20*spdiag(ones(m,1),0)+5/4*spdiag(ones(m-1,1),+1)-1/2*spdiag(ones(m-2,1),+2)+1/6*spdiag(ones(m-3,1),+3)-1/28*spdiag(ones(m-4,1),+4)+1/280*spdiag(ones(m-5,1),+5)); - %Q_U = [-0.375523e6 / 0.237369744e9 0.73393569793e11 / 0.111599208000e12 -0.1244182786433e13 / 0.14954293872000e14 -0.245555527093e12 / 0.2492382312000e13 0.131952309637e12 / 0.14954293872000e14 0.62907543023e11 / 0.3738573468000e13; -0.4818422888131e13 / 0.7477146936000e13 -0.10098911e8 / 0.307701520e9 0.796148045807e12 / 0.1495429387200e13 0.487507785187e12 / 0.2990858774400e13 0.8887918703e10 / 0.249238231200e12 -0.804641655323e12 / 0.14954293872000e14; 0.925771462433e12 / 0.14954293872000e14 -0.601394962127e12 / 0.1495429387200e13 -0.674340047e9 / 0.4153970520e10 0.560418114917e12 / 0.747714693600e12 -0.1043208352973e13 / 0.2990858774400e13 0.42207619661e11 / 0.356054616000e12; 0.285728443093e12 / 0.2492382312000e13 -0.873874460227e12 / 0.2990858774400e13 -0.34960363091e11 / 0.106816384800e12 -0.1528791703e10 / 0.4153970520e10 0.1741615892927e13 / 0.1495429387200e13 -0.6107723941553e13 / 0.14954293872000e14; -0.213784143637e12 / 0.14954293872000e14 0.1011411511e10 / 0.35605461600e11 0.141331778093e12 / 0.2990858774400e13 -0.632571739487e12 / 0.1495429387200e13 -0.1464298201e10 / 0.2769313680e10 0.9523869708211e13 / 0.7477146936000e13; -0.60750297023e11 / 0.3738573468000e13 0.614597810123e12 / 0.14954293872000e14 -0.16980260027e11 / 0.2492382312000e13 -0.343200536047e12 / 0.14954293872000e14 -0.389757912373e12 / 0.1068163848000e13 -0.964053659e9 / 0.1661588208e10;]; - Q_U =[-0.16683e5 / 0.63018560e8 0.29345822969987e14 / 0.43354248537600e14 -0.2734625426591e13 / 0.40644608004000e14 -0.113480208109603e15 / 0.780376473676800e15 -0.830250230261e12 / 0.26012549122560e14 0.2500519492033e13 / 0.32515686403200e14 0.2235718279643e13 / 0.390188236838400e15 -0.388481888477e12 / 0.26543417472000e14; -0.29227665839987e14 / 0.43354248537600e14 -0.493793e6 / 0.63018560e8 0.8302717120817e13 / 0.26543417472000e14 0.3739408501537e13 / 0.9290196115200e13 0.2684481534461e13 / 0.13935294172800e14 -0.4450185662513e13 / 0.18580392230400e14 -0.1221838279381e13 / 0.37160784460800e14 0.90595000956023e14 / 0.1950941184192000e16; 0.2505689537591e13 / 0.40644608004000e14 -0.7312922392817e13 / 0.26543417472000e14 -0.69332623e8 / 0.1323389760e10 0.10994933811709e14 / 0.18580392230400e14 -0.9270952411151e13 / 0.18580392230400e14 0.3191238635141e13 / 0.20644880256000e14 0.4442211176987e13 / 0.92901961152000e14 -0.940661365031e12 / 0.32515686403200e14; 0.118016946570403e15 / 0.780376473676800e15 -0.4173878828737e13 / 0.9290196115200e13 -0.7990503962509e13 / 0.18580392230400e14 -0.207799621e9 / 0.1323389760e10 0.2044021560341e13 / 0.2477385630720e13 0.511197701761e12 / 0.18580392230400e14 0.1237681717213e13 / 0.13935294172800e14 -0.7784834666617e13 / 0.130062745612800e15; 0.68609076271e11 / 0.2364777192960e13 -0.2235651762161e13 / 0.13935294172800e14 0.6527681584751e13 / 0.18580392230400e14 -0.1115980068821e13 / 0.2477385630720e13 -0.55386253e8 / 0.189055680e9 0.3208334350649e13 / 0.3716078446080e13 -0.407569013461e12 / 0.844563283200e12 0.136474842626653e15 / 0.780376473676800e15; -0.2487637785013e13 / 0.32515686403200e14 0.4244231077313e13 / 0.18580392230400e14 -0.1550378843141e13 / 0.20644880256000e14 -0.5726967564961e13 / 0.18580392230400e14 -0.1017898941929e13 / 0.3716078446080e13 -0.526653889e9 / 0.1323389760e10 0.45241297077547e14 / 0.37160784460800e14 -0.2279608411897e13 / 0.5080576000500e13; -0.2164019088443e13 / 0.390188236838400e15 0.1263196075861e13 / 0.37160784460800e14 -0.6600697610987e13 / 0.92901961152000e14 0.556610591687e12 / 0.13935294172800e14 0.926842346471e12 / 0.9290196115200e13 -0.18757693936747e14 / 0.37160784460800e14 -0.584765899e9 / 0.1323389760e10 0.204646287449e12 / 0.168431424000e12; 0.387091928477e12 / 0.26543417472000e14 -0.90231551688023e14 / 0.1950941184192000e16 0.1032404418251e13 / 0.32515686403200e14 0.3502353445417e13 / 0.130062745612800e15 -0.15385068876253e14 / 0.780376473676800e15 0.262499068919e12 / 0.10161152001000e14 -0.867004691939e12 / 0.1852745664000e13 -0.85017967e8 / 0.189055680e9;]; + Q_U =[ + -0.16683e5/0.63018560e8 0.29345822969987e14/0.43354248537600e14 -0.2734625426591e13/0.40644608004000e14 -0.113480208109603e15/0.780376473676800e15 -0.830250230261e12/0.26012549122560e14 0.2500519492033e13/0.32515686403200e14 0.2235718279643e13/0.390188236838400e15 -0.388481888477e12/0.26543417472000e14; + -0.29227665839987e14/0.43354248537600e14 -0.493793e6/0.63018560e8 0.8302717120817e13/0.26543417472000e14 0.3739408501537e13/0.9290196115200e13 0.2684481534461e13/0.13935294172800e14 -0.4450185662513e13/0.18580392230400e14 -0.1221838279381e13/0.37160784460800e14 0.90595000956023e14/0.1950941184192000e16; + 0.2505689537591e13/0.40644608004000e14 -0.7312922392817e13/0.26543417472000e14 -0.69332623e8/0.1323389760e10 0.10994933811709e14/0.18580392230400e14 -0.9270952411151e13/0.18580392230400e14 0.3191238635141e13/0.20644880256000e14 0.4442211176987e13/0.92901961152000e14 -0.940661365031e12/0.32515686403200e14; + 0.118016946570403e15/0.780376473676800e15 -0.4173878828737e13/0.9290196115200e13 -0.7990503962509e13/0.18580392230400e14 -0.207799621e9/0.1323389760e10 0.2044021560341e13/0.2477385630720e13 0.511197701761e12/0.18580392230400e14 0.1237681717213e13/0.13935294172800e14 -0.7784834666617e13/0.130062745612800e15; + 0.68609076271e11/0.2364777192960e13 -0.2235651762161e13/0.13935294172800e14 0.6527681584751e13/0.18580392230400e14 -0.1115980068821e13/0.2477385630720e13 -0.55386253e8/0.189055680e9 0.3208334350649e13/0.3716078446080e13 -0.407569013461e12/0.844563283200e12 0.136474842626653e15/0.780376473676800e15; + -0.2487637785013e13/0.32515686403200e14 0.4244231077313e13/0.18580392230400e14 -0.1550378843141e13/0.20644880256000e14 -0.5726967564961e13/0.18580392230400e14 -0.1017898941929e13/0.3716078446080e13 -0.526653889e9/0.1323389760e10 0.45241297077547e14/0.37160784460800e14 -0.2279608411897e13/0.5080576000500e13; + -0.2164019088443e13/0.390188236838400e15 0.1263196075861e13/0.37160784460800e14 -0.6600697610987e13/0.92901961152000e14 0.556610591687e12/0.13935294172800e14 0.926842346471e12/0.9290196115200e13 -0.18757693936747e14/0.37160784460800e14 -0.584765899e9/0.1323389760e10 0.204646287449e12/0.168431424000e12; + 0.387091928477e12/0.26543417472000e14 -0.90231551688023e14/0.1950941184192000e16 0.1032404418251e13/0.32515686403200e14 0.3502353445417e13/0.130062745612800e15 -0.15385068876253e14/0.780376473676800e15 0.262499068919e12/0.10161152001000e14 -0.867004691939e12/0.1852745664000e13 -0.85017967e8/0.189055680e9; + ]; Qp(1:8,1:8)=Q_U; - Qp(m-7:m,m-7:m)=flipud( fliplr(Q_U(1:8,1:8) ) )'; %%% This is different from standard SBP + Qp(m-7:m,m-7:m)=rot90(Q_U,2)'; %%% This is different from standard SBP Qm=-Qp'; - e_1=zeros(m,1);e_1(1)=1; - e_m=zeros(m,1);e_m(m)=1; + e_1=sparse(m,1);e_1(1)=1; + e_m=sparse(m,1);e_m(m)=1; Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ;