Mercurial > repos > public > sbplib
changeset 267:f7ac3cd6eeaa operator_remake
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
line wrap: on
line diff
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_10.m --- a/+sbp/+implementations/d1_noneq_10.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_10.m Fri Sep 09 14:53:41 2016 +0200 @@ -30,7 +30,7 @@ x9 = 8.4309689056870e+00; x10 = 9.4309689056870e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -45,7 +45,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.0000000000000e-01; P1 = 5.8980851260667e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_12.m --- a/+sbp/+implementations/d1_noneq_12.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_12.m Fri Sep 09 14:53:41 2016 +0200 @@ -32,7 +32,7 @@ x11 = 1.1000000000000e+01; x12 = 1.2000000000000e+01; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -47,7 +47,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.0000000000011e-01; P1 = 5.9616216757547e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_4.m --- a/+sbp/+implementations/d1_noneq_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -24,7 +24,7 @@ x3 = 2.8022115125776e+00; x4 = 3.8022115125776e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -39,7 +39,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 2.1259737557798e-01; P1 = 1.0260290400758e+00;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_6.m --- a/+sbp/+implementations/d1_noneq_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -26,7 +26,7 @@ x5 = 4.2638953951239e+00; x6 = 5.2638953951239e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -41,7 +41,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.3030223027124e-01; P1 = 6.8851501587715e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_8.m --- a/+sbp/+implementations/d1_noneq_8.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_8.m Fri Sep 09 14:53:41 2016 +0200 @@ -28,7 +28,7 @@ x7 = 6.3192851303204e+00; x8 = 7.3192851303204e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -43,7 +43,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.0758368078310e-01; P1 = 6.1909685107891e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_minimal_10.m --- a/+sbp/+implementations/d1_noneq_minimal_10.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_10.m Fri Sep 09 14:53:41 2016 +0200 @@ -28,7 +28,7 @@ x7 = 7.0000000000000e+00; x8 = 8.0000000000000e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -43,7 +43,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.6717213975289e-01; P1 = 9.3675739171278e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_minimal_12.m --- a/+sbp/+implementations/d1_noneq_minimal_12.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_12.m Fri Sep 09 14:53:41 2016 +0200 @@ -30,7 +30,7 @@ x9 = 9.0000000000000e+00; x10 = 1.0000000000000e+01; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -45,7 +45,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.3013597111750e-01; P1 = 7.6146045079020e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_minimal_4.m --- a/+sbp/+implementations/d1_noneq_minimal_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -23,7 +23,7 @@ x2 = 1.7712298784256e+00; x3 = 2.7712298784256e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -38,7 +38,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 2.6864248295847e-01; P1 = 1.0094667153500e+00;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_minimal_6.m --- a/+sbp/+implementations/d1_noneq_minimal_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -25,7 +25,7 @@ x4 = 3.1968523189207e+00; x5 = 4.1968523189207e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -40,7 +40,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.2740260779883e-01; P1 = 6.1820981002054e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_noneq_minimal_8.m --- a/+sbp/+implementations/d1_noneq_minimal_8.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_noneq_minimal_8.m Fri Sep 09 14:53:41 2016 +0200 @@ -26,7 +26,7 @@ x5 = 4.4051531374839e+00; x6 = 5.4051531374839e+00; -xb = zeros(m+1,1); +xb = sparse(m+1,1); for i = 0:m xb(i+1) = eval(['x' num2str(i)]); end @@ -41,7 +41,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% -P = zeros(BP,1); +P = sparse(BP,1); %#ok<*NASGU> P0 = 1.4523997892351e-01; P1 = 7.6864793350174e-01;
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_2.m --- a/+sbp/+implementations/d1_upwind_2.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_2.m Fri Sep 09 14:53:41 2016 +0200 @@ -25,11 +25,11 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_3.m --- a/+sbp/+implementations/d1_upwind_3.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_3.m Fri Sep 09 14:53:41 2016 +0200 @@ -30,8 +30,8 @@ e_m=sparse(m,1); e_m(m)=1; - Dp=HI*(Qp-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_4.m --- a/+sbp/+implementations/d1_upwind_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -30,7 +30,7 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_5.m --- a/+sbp/+implementations/d1_upwind_5.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_5.m Fri Sep 09 14:53:41 2016 +0200 @@ -23,14 +23,14 @@ ]; 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(1:4,1:4) ,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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_6.m --- a/+sbp/+implementations/d1_upwind_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -32,8 +32,8 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_7.m --- a/+sbp/+implementations/d1_upwind_7.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_7.m Fri Sep 09 14:53:41 2016 +0200 @@ -31,10 +31,10 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_8.m --- a/+sbp/+implementations/d1_upwind_8.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_8.m Fri Sep 09 14:53:41 2016 +0200 @@ -37,7 +37,7 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d1_upwind_9.m --- a/+sbp/+implementations/d1_upwind_9.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d1_upwind_9.m Fri Sep 09 14:53:41 2016 +0200 @@ -34,10 +34,10 @@ 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') ; + Dp=HI*(Qp-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; - Dm=HI*(Qm-1/2*e_1*e_1'+1/2*e_m*e_m') ; + Dm=HI*(Qm-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; end
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_10.m --- a/+sbp/+implementations/d2_10.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_10.m Fri Sep 09 14:53:41 2016 +0200 @@ -7,26 +7,32 @@ H_U = [0.5261271563e10 / 0.18289152000e11 0 0 0 0 0 0 0 0 0 0; 0 0.2881040311e10 / 0.1828915200e10 0 0 0 0 0 0 0 0 0; 0 0 0.52175551e8 / 0.406425600e9 0 0 0 0 0 0 0 0; 0 0 0 0.11662993e8 / 0.6096384e7 0 0 0 0 0 0 0; 0 0 0 0 0.50124587e8 / 0.87091200e8 0 0 0 0 0 0; 0 0 0 0 0 0.50124587e8 / 0.72576000e8 0 0 0 0 0; 0 0 0 0 0 0 0.148333439e9 / 0.87091200e8 0 0 0 0; 0 0 0 0 0 0 0 0.63867949e8 / 0.152409600e9 0 0 0; 0 0 0 0 0 0 0 0 0.20608675e8 / 0.16257024e8 0 0; 0 0 0 0 0 0 0 0 0 0.1704508063e10 / 0.1828915200e10 0; 0 0 0 0 0 0 0 0 0 0 0.18425967263e11 / 0.18289152000e11;]; - H=eye(m); + H=speye(m); H(1:11,1:11)=H_U; - H(m-10:m,m-10:m)=flipud( fliplr(H_U(1:11,1:11) ) ); + H(m-10:m,m-10:m)=rot90( H_U(1:11,1:11) ,2 ); H=H*h; HI=inv(H); - Q=-(-1/1260*diag(ones(m-5,1),5)+5/504*diag(ones(m-4,1),4)-5/84*diag(ones(m-3,1),3)+5/21*diag(ones(m-2,1),2)-5/6*diag(ones(m-1,1),1)+5/6*diag(ones(m-1,1),-1)-5/21*diag(ones(m-2,1),-2)+5/84*diag(ones(m-3,1),-3)-5/504*diag(ones(m-4,1),-4)+1/1260*diag(ones(m-5,1),-5)); - +% Q=-(-1/1260*diag(ones(m-5,1),5)+5/504*diag(ones(m-4,1),4)-5/84*diag(ones(m-3,1),3)+5/21*diag(ones(m-2,1),2)-5/6*diag(ones(m-1,1),1)+5/6*diag(ones(m-1,1),-1)-5/21*diag(ones(m-2,1),-2)+5/84*diag(ones(m-3,1),-3)-5/504*diag(ones(m-4,1),-4)+1/1260*diag(ones(m-5,1),-5)); + + diags = -5:5; + stencil = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; + Q = stripeMatrix(stencil, diags, m); + + + Q_U = [-0.1e1 / 0.2e1 0.2300876759589119e16 / 0.3395198177280000e16 -0.99808615498093e14 / 0.2263465451520000e16 -0.34957747037683e14 / 0.212199886080000e15 -0.709586095717e12 / 0.13473008640000e14 0.325330433051e12 / 0.6218311680000e13 0.27953548723573e14 / 0.485028311040000e15 0.2690678501e10 / 0.412439040000e12 -0.2397491025029e13 / 0.70733295360000e14 -0.9959492094287e13 / 0.1131732725760000e16 0.5242772857661e13 / 0.522338181120000e15; -0.2300876759589119e16 / 0.3395198177280000e16 0 0.3103439505511e13 / 0.16643128320000e14 0.2700334568377e13 / 0.5052378240000e13 0.50587599589937e14 / 0.242514155520000e15 -0.5570893587157e13 / 0.40419025920000e14 -0.1496329934863e13 / 0.8083805184000e13 -0.322512443237e12 / 0.12482346240000e14 0.2275340833763e13 / 0.23096586240000e14 0.22922115021893e14 / 0.848799544320000e15 -0.143e3 / 0.5000e4; 0.99808615498093e14 / 0.2263465451520000e16 -0.3103439505511e13 / 0.16643128320000e14 0 0.15053664233879e14 / 0.40419025920000e14 -0.9306441440587e13 / 0.32335220736000e14 -0.945459729233e12 / 0.13473008640000e14 0.956829267413e12 / 0.5774146560000e13 0.446866085903e12 / 0.56586636288000e14 -0.41109372242993e14 / 0.754488483840000e15 0.1e1 / 0.500e3 0.17e2 / 0.2500e4; 0.34957747037683e14 / 0.212199886080000e15 -0.2700334568377e13 / 0.5052378240000e13 -0.15053664233879e14 / 0.40419025920000e14 0 0.3899174751943e13 / 0.10104756480000e14 0.4691717443831e13 / 0.10104756480000e14 -0.58571891887e11 / 0.396264960000e12 0.100791910589e12 / 0.1040195520000e13 -0.425149181e9 / 0.29719872000e11 -0.2376515922259e13 / 0.30314269440000e14 0.36894656431e11 / 0.1036385280000e13; 0.709586095717e12 / 0.13473008640000e14 -0.50587599589937e14 / 0.242514155520000e15 0.9306441440587e13 / 0.32335220736000e14 -0.3899174751943e13 / 0.10104756480000e14 0 -0.4552305973e10 / 0.444165120000e12 0.4984940784247e13 / 0.11548293120000e14 -0.19410791e8 / 0.146764800e9 -0.2912773695913e13 / 0.40419025920000e14 0.127067639161e12 / 0.3233522073600e13 -0.89277540287e11 / 0.37309870080000e14; -0.325330433051e12 / 0.6218311680000e13 0.5570893587157e13 / 0.40419025920000e14 0.945459729233e12 / 0.13473008640000e14 -0.4691717443831e13 / 0.10104756480000e14 0.4552305973e10 / 0.444165120000e12 0 0.31722122083e11 / 0.84913920000e11 -0.887187251021e12 / 0.10104756480000e14 -0.1661755478749e13 / 0.26946017280000e14 0.1505713246249e13 / 0.13473008640000e14 -0.38859042469e11 / 0.1036385280000e13; -0.27953548723573e14 / 0.485028311040000e15 0.1496329934863e13 / 0.8083805184000e13 -0.956829267413e12 / 0.5774146560000e13 0.58571891887e11 / 0.396264960000e12 -0.4984940784247e13 / 0.11548293120000e14 -0.31722122083e11 / 0.84913920000e11 0 0.9357094407023e13 / 0.20209512960000e14 0.52602356173249e14 / 0.161676103680000e15 -0.1435252677707e13 / 0.17322439680000e14 -0.33048158431e11 / 0.3109155840000e13; -0.2690678501e10 / 0.412439040000e12 0.322512443237e12 / 0.12482346240000e14 -0.446866085903e12 / 0.56586636288000e14 -0.100791910589e12 / 0.1040195520000e13 0.19410791e8 / 0.146764800e9 0.887187251021e12 / 0.10104756480000e14 -0.9357094407023e13 / 0.20209512960000e14 0 0.70089734285659e14 / 0.141466590720000e15 -0.105938137621e12 / 0.471555302400e12 0.4358988450443e13 / 0.65292272640000e14; 0.2397491025029e13 / 0.70733295360000e14 -0.2275340833763e13 / 0.23096586240000e14 0.41109372242993e14 / 0.754488483840000e15 0.425149181e9 / 0.29719872000e11 0.2912773695913e13 / 0.40419025920000e14 0.1661755478749e13 / 0.26946017280000e14 -0.52602356173249e14 / 0.161676103680000e15 -0.70089734285659e14 / 0.141466590720000e15 0 0.314274398580227e15 / 0.377244241920000e15 -0.97822819709e11 / 0.487710720000e12; 0.9959492094287e13 / 0.1131732725760000e16 -0.22922115021893e14 / 0.848799544320000e15 -0.1e1 / 0.500e3 0.2376515922259e13 / 0.30314269440000e14 -0.127067639161e12 / 0.3233522073600e13 -0.1505713246249e13 / 0.13473008640000e14 0.1435252677707e13 / 0.17322439680000e14 0.105938137621e12 / 0.471555302400e12 -0.314274398580227e15 / 0.377244241920000e15 0 0.7519148725913e13 / 0.9327467520000e13; -0.5242772857661e13 / 0.522338181120000e15 0.143e3 / 0.5000e4 -0.17e2 / 0.2500e4 -0.36894656431e11 / 0.1036385280000e13 0.89277540287e11 / 0.37309870080000e14 0.38859042469e11 / 0.1036385280000e13 0.33048158431e11 / 0.3109155840000e13 -0.4358988450443e13 / 0.65292272640000e14 0.97822819709e11 / 0.487710720000e12 -0.7519148725913e13 / 0.9327467520000e13 0;]; Q(1:11,1:11)=Q_U; - Q(m-10:m,m-10:m)=flipud( fliplr(-Q_U ) ); + Q(m-10:m,m-10:m)=rot90( -Q_U ,2 ); - D1=HI*Q; + D1=H\Q; - s18=-1.000000000000000; s19=0.195000000000000; % alpha 0.0605 +% s18=-1.000000000000000; s19=0.195000000000000; % alpha 0.0605 s18= -0.475000000000000; s19=0.110000000000000; % alpha 0.0350 %s18=0;s19=0; - DS=zeros(m,m); + DS=sparse(m,m); DS(1,1:9)=[0.49e2 / 0.20e2 - s18 - (7 * s19) -0.6e1 + 0.7e1 * s18 + (48 * s19) 0.15e2 / 0.2e1 - 0.21e2 * s18 - (140 * s19) -0.20e2 / 0.3e1 + 0.35e2 * s18 + (224 * s19) 0.15e2 / 0.4e1 - 0.35e2 * s18 - (210 * s19) -0.6e1 / 0.5e1 + 0.21e2 * s18 + (112 * s19) 0.1e1 / 0.6e1 - 0.7e1 * s18 - (28 * s19) s18 s19;]; DS(m,m-8:m)=fliplr(DS(1,1:9)); DS=DS/h; @@ -34,22 +40,27 @@ M_U = [0.12056593789671863908e1 -0.13378814169347239658e1 0.36847309286546532061e-2 0.15698288365600946515e0 -0.37472461482539197952e-2 -0.62491712449361657064e-2 -0.29164045872729581661e-1 0.54848184117832929161e-3 0.13613461413384884448e-1 -0.25059220258337808220e-2 -0.94113457993630916498e-3; -0.13378814169347239658e1 0.21749807117105597139e1 -0.12369059547124894597e0 -0.83712574037924152603e0 0.50065127254670973258e-1 0.81045853127317536361e-2 0.97405846039248226536e-1 -0.68942461520402214720e-3 -0.41326971493379188475e-1 0.75778529605774119402e-2 0.25800256160095691057e-2; 0.36847309286546532061e-2 -0.12369059547124894597e0 0.18361596652499065332e0 0.48289690013342693109e-1 -0.19719621435164680412e0 0.11406859029505842791e0 -0.29646295985488126964e-1 -0.16038463172861201306e-2 0.32879841528337653050e-2 -0.93242311589807387463e-3 0.12241332668787820533e-3; 0.15698288365600946515e0 -0.83712574037924152603e0 0.48289690013342693109e-1 0.12886524606662484673e1 -0.14403037739488789185e0 -0.44846291607489015475e0 -0.10598334599408054277e0 -0.15873275740355918053e-1 0.73988493386459608166e-1 -0.12508848749152899785e-1 -0.39290233894513005339e-2; -0.37472461482539197952e-2 0.50065127254670973258e-1 -0.19719621435164680412e0 -0.14403037739488789185e0 0.51482665719685186210e0 0.51199577887125103015e-1 -0.36233561810883077365e0 0.91356850268746392169e-1 0.24195916108052419451e-2 -0.18564214413731389338e-2 -0.70192677320704413827e-3; -0.62491712449361657064e-2 0.81045853127317536361e-2 0.11406859029505842791e0 -0.44846291607489015475e0 0.51199577887125103015e-1 0.68636003380365860083e0 -0.28358848290867614908e0 -0.13836006478253396528e0 0.76158070663111995297e-2 0.11447010307180005164e-1 -0.21349696610286552676e-2; -0.29164045872729581661e-1 0.97405846039248226536e-1 -0.29646295985488126964e-1 -0.10598334599408054277e0 -0.36233561810883077365e0 -0.28358848290867614908e0 0.15216081480839085990e1 -0.42653865162216293237e0 -0.42047484981879143123e0 0.19813359263872926304e-1 0.19221397241190103344e-1; 0.54848184117832929161e-3 -0.68942461520402214720e-3 -0.16038463172861201306e-2 -0.15873275740355918053e-1 0.91356850268746392169e-1 -0.13836006478253396528e0 -0.42653865162216293237e0 0.10656733504627815335e1 -0.66921872668484232217e0 0.12022033144141336599e0 -0.30157881394591483631e-1; 0.13613461413384884448e-1 -0.41326971493379188475e-1 0.32879841528337653050e-2 0.73988493386459608166e-1 0.24195916108052419451e-2 0.76158070663111995297e-2 -0.42047484981879143123e0 -0.66921872668484232217e0 0.24064247712949611684e1 -0.15150200315922263367e1 0.17373015320416595052e0; -0.25059220258337808220e-2 0.75778529605774119402e-2 -0.93242311589807387463e-3 -0.12508848749152899785e-1 -0.18564214413731389338e-2 0.11447010307180005164e-1 0.19813359263872926304e-1 0.12022033144141336599e0 -0.15150200315922263367e1 0.27682502485427255096e1 -0.15975407111468405444e1; -0.94113457993630916498e-3 0.25800256160095691057e-2 0.12241332668787820533e-3 -0.39290233894513005339e-2 -0.70192677320704413827e-3 -0.21349696610286552676e-2 0.19221397241190103344e-1 -0.30157881394591483631e-1 0.17373015320416595052e0 -0.15975407111468405444e1 0.29033627686681129471e1;]; - M=-(1/3150)*diag(ones(m-5,1),5)+(5/1008)*diag(ones(m-4,1),4)-(5/126)*diag(ones(m-3,1),3)+(5/21)*diag(ones(m-2,1),2)-(5/3)*diag(ones(m-1,1),1)... - -(1/3150)*diag(ones(m-5,1),-5)+(5/1008)*diag(ones(m-4,1),-4)-(5/126)*diag(ones(m-3,1),-3)+(5/21)*diag(ones(m-2,1),-2)-(5/3)*diag(ones(m-1,1),-1)... - +(5269/1800)*diag(ones(m,1),0); +% M=-(1/3150)*diag(ones(m-5,1),5)+(5/1008)*diag(ones(m-4,1),4)-(5/126)*diag(ones(m-3,1),3)+(5/21)*diag(ones(m-2,1),2)-(5/3)*diag(ones(m-1,1),1)... +% -(1/3150)*diag(ones(m-5,1),-5)+(5/1008)*diag(ones(m-4,1),-4)-(5/126)*diag(ones(m-3,1),-3)+(5/21)*diag(ones(m-2,1),-2)-(5/3)*diag(ones(m-1,1),-1)... +% +(5269/1800)*diag(ones(m,1),0); + + diags = -5:5; + left_stencil = [-1/3150,5/1008,-5/126,5/21,-5/3]; + stencil = [left_stencil,5269/1800,fliplr(left_stencil)]; + M = stripeMatrix(stencil, diags, m); M(1:11,1:11)=M_U; - M(m-10:m,m-10:m)=flipud( fliplr(M_U(1:11,1:11) ) ); + M(m-10:m,m-10:m)=rot90( M_U(1:11,1:11) ,2 ); - D2=HI*(-M/h+DS); + D2=H\(-M/h+DS); - e_1 = zeros(m,1); + e_1 = sparse(m,1); e_1(1)= 1; - e_m = zeros(m,1); + e_m = sparse(m,1); e_m(end)= 1; S_1 = -DS(1,:)'; S_m = DS(end,:)'; - Q = H*D1-(-e_1*e_1' + e_m*e_m'); + Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_12.m --- a/+sbp/+implementations/d2_12.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_12.m Fri Sep 09 14:53:41 2016 +0200 @@ -7,7 +7,7 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end -H=diag(ones(m,1),0); +H=speye(m,m); H(1:15,1:15)=diag([2.880607858916397e-01,... 1.560376162339675e+00,... 2.403139445479289e-01,... @@ -29,14 +29,9 @@ H=H*h; HI = inv(H); -a1 = 6/7; a2 = -15/56; a3 = 5/63; a4 = -1/56; a5 = 1/385; a6 = -1/5544; - -D1=a6*(diag(ones(m-6,1),6)-diag(ones(m-6,1),-6))+... - a5*(diag(ones(m-5,1),5)-diag(ones(m-5,1),-5))+... - a4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+... - a3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+... - a2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+... - a1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +diags = -6:6; +stencil = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; +D1 = stripeMatrix(stencil, diags, m); D1(1, 1)= -1.735744761135539e+00; @@ -281,14 +276,14 @@ D2=D1*D1; -e_1 = zeros(m,1); +e_1 = sparse(m,1); e_1(1)= 1; -e_m = zeros(m,1); +e_m = sparse(m,1); e_m(end)= 1; S_1 = (e_1'*D1)'; S_m = (e_m'*D1)'; -Q = H*D1-(-e_1*e_1' + e_m*e_m'); +Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_2.m --- a/+sbp/+implementations/d2_2.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_2.m Fri Sep 09 14:53:41 2016 +0200 @@ -5,29 +5,35 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - 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; - H=(eye(m,m));H(1,1)=0.5;H(m,m)=0.5; + H=(speye(m,m));H(1,1)=0.5;H(m,m)=0.5; H=h*H; HI=inv(H); - - D1=((.5*diag(ones(m-1,1),1)-.5*diag(ones(m-1,1),-1))); + + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrix(stencil, diags, m); + D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; D1(m,m-1)=-1;D1(m,m)=1; D1=D1/h; - Q=H*D1 + 1/2*e_1*e_1' - 1/2*e_m*e_m'; + Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); - D2=((diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)-2*diag(ones(m,1),0))); + diags = -1:1; + stencil = [1 -2 1]; + D2 = stripeMatrix(stencil, diags, m); + D2(1,1)=1;D2(1,2)=-2;D2(1,3)=1; D2(m,m-2)=1;D2(m,m-1)=-2;D2(m,m)=1; D2=D2/h^2; S_U=[-3/2, 2, -1/2]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:3)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-2:m)=fliplr(-S_U);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_4.m --- a/+sbp/+implementations/d2_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -5,12 +5,12 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - 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; e=ones(m,1); H=spdiags(e,0,m,m); - %H=diag(ones(m,1),0); + %H=speye(m,m); H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]); H(m-3:m,m-3:m)=rot90(diag([17/48 59/48 43/48 49/48]),2); @@ -24,7 +24,7 @@ Q(1:4,1:4)=Q_U; Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; @@ -37,13 +37,13 @@ M=M/h; S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; - S_1=sparse(zeros(1,m)); + S_1=sparse(sparse(1,m)); S_1(1:4)=S_U; - S_m=sparse(zeros(1,m)); + S_m=sparse(sparse(1,m)); S_m(m-3:m)=fliplr(-S_U); - D2=HI*(-M-e_1*S_1+e_m*S_m); + D2=H\(-M-e_1*S_1+e_m*S_m); S_1 = S_1'; S_m = S_m'; end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_6.m --- a/+sbp/+implementations/d2_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -5,10 +5,10 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - 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; - H=diag(ones(m,1),0); + H=speye(m,m); H(1:6,1:6)=diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, ... 43801/43200]); H(m-5:m,m-5:m)=rot90(diag([13649/43200,12013/8640, ... @@ -22,11 +22,9 @@ % Ett optimerat varde ger x1=0.70127127127127 = 331/472 x1=0.70127127127127; - - - D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); - - + diags = -3:3; + stencil = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; + D1 = stripeMatrix(stencil, diags, m); D1(1:6,1:9)=[-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ... 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ... @@ -37,19 +35,22 @@ D1(m-5:m,m-8:m)=rot90( -D1(1:6,1:9),2); D1=D1/h; - Q=H*D1 + 1/2*e_1*e_1' - 1/2*e_m*e_m'; + Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); - D2=(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180; - + %D2=(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180; + diags = -3:3; + stencil = 1/180*[2,-27,270,-490,270,-27,2]; + D2 = stripeMatrix(stencil, diags, m); + D2(1:6,1:9)=[114170/40947, -438107/54596, 336409/40947, -276997/81894, 3747/13649, 21035/163788, 0, 0, 0;6173/5860, -2066/879, 3283/1758, -303/293, 2111/3516, -601/4395, 0, 0, 0;-52391/81330, 134603/32532, -21982/2711, 112915/16266, -46969/16266, 30409/54220, 0, 0, 0;68603/321540, -12423/10718, 112915/32154, -75934/16077, 53369/21436, -54899/160770, 48/5359, 0, 0;-7053/39385, 86551/94524, -46969/23631, 53369/15754, -87904/23631, 820271/472620, -1296/7877, 96/7877, 0;21035/525612, -24641/131403, 30409/87602, -54899/131403, 820271/525612, -117600/43801, 64800/43801, -6480/43801, 480/43801]; D2(m-5:m,m-8:m)=rot90( D2(1:6,1:9) ,2 ); D2=D2/h^2; S_U=[-25/12, 4, -3, 4/3, -1/4]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:5)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-4:m)=fliplr(-S_U);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_8.m --- a/+sbp/+implementations/d2_8.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_8.m Fri Sep 09 14:53:41 2016 +0200 @@ -6,13 +6,15 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H(1:8,1:8)=diag([1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]); - H(m-7:m,m-7:m)=fliplr(flipud(diag([1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]))); + H(m-7:m,m-7:m)=rot90(diag([1498139/5080320, 1107307/725760, 20761/80640, 1304999/725760, 299527/725760, 103097/80640, 670091/725760, 5127739/5080320]),2); - D1=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4)); +% D1=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4)); - + diags = -4:4; + stencil = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; + D1 = stripeMatrix(stencil, diags, m); %r68 = -1022551/30481920; @@ -62,7 +64,7 @@ D1(1:8,1:12)=[-2540160/1498139, -142642467/5992556+50803200/1498139*r78+5080320/1498139*r67+25401600/1498139*r68, 705710031/5992556-228614400/1498139*r78-25401600/1498139*r67-121927680/1498139*r68, -3577778591/17977668+381024000/1498139*r78+50803200/1498139*r67+228614400/1498139*r68, 203718909/1498139-254016000/1498139*r78-50803200/1498139*r67-203212800/1498139*r68, -32111205/5992556+25401600/1498139*r67+76204800/1498139*r68, -652789417/17977668+76204800/1498139*r78-5080320/1498139*r67, 74517981/5992556-25401600/1498139*r78-5080320/1498139*r68, 0, 0, 0, 0;142642467/31004596-7257600/1107307*r78-725760/1107307*r67-3628800/1107307*r68, 0, -141502371/2214614+91445760/1107307*r78+10886400/1107307*r67+50803200/1107307*r68, 159673719/1107307-203212800/1107307*r78-29030400/1107307*r67-127008000/1107307*r68, -1477714693/13287684+152409600/1107307*r78+32659200/1107307*r67+127008000/1107307*r68, 11652351/2214614-17418240/1107307*r67-50803200/1107307*r68, 36069450/1107307-50803200/1107307*r78+3628800/1107307*r67, -536324953/46506894+17418240/1107307*r78+3628800/1107307*r68, 0, 0, 0, 0;-18095129/134148+3628800/20761*r78+403200/20761*r67+1935360/20761*r68, 47167457/124566-10160640/20761*r78-1209600/20761*r67-5644800/20761*r68, 0, -120219461/124566+25401600/20761*r78+4032000/20761*r67+16934400/20761*r68, 249289259/249132-25401600/20761*r78-6048000/20761*r67-22579200/20761*r68, -2611503/41522+3628800/20761*r67+10160640/20761*r68, -7149666/20761+10160640/20761*r78-806400/20761*r67, 37199165/290654-3628800/20761*r78-806400/20761*r68, 0, 0, 0, 0;3577778591/109619916-54432000/1304999*r78-7257600/1304999*r67-32659200/1304999*r68, -159673719/1304999+203212800/1304999*r78+29030400/1304999*r67+127008000/1304999*r68, 360658383/2609998-228614400/1304999*r78-36288000/1304999*r67-152409600/1304999*r68, 0, -424854441/5219996+127008000/1304999*r78+36288000/1304999*r67+127008000/1304999*r68, 22885113/2609998-29030400/1304999*r67-76204800/1304999*r68, 158096578/3914997-76204800/1304999*r78+7257600/1304999*r67, -296462325/18269986+29030400/1304999*r78+7257600/1304999*r68, 0, 0, 0, 0;-203718909/2096689+36288000/299527*r78+7257600/299527*r67+29030400/299527*r68, 1477714693/3594324-152409600/299527*r78-32659200/299527*r67-127008000/299527*r68, -747867777/1198108+228614400/299527*r78+54432000/299527*r67+203212800/299527*r68, 424854441/1198108-127008000/299527*r78-36288000/299527*r67-127008000/299527*r68, 0, -17380335/1198108+10886400/299527*r67+25401600/299527*r68, -67080435/1198108+25401600/299527*r78-3628800/299527*r67, 657798011/25160268-10886400/299527*r78-3628800/299527*r68, -2592/299527, 0, 0, 0;1529105/1237164-403200/103097*r67-1209600/103097*r68, -3884117/618582+1935360/103097*r67+5644800/103097*r68, 2611503/206194-3628800/103097*r67-10160640/103097*r68, -7628371/618582+3225600/103097*r67+8467200/103097*r68, 5793445/1237164-1209600/103097*r67-2822400/103097*r68, 0, 80640/103097*r67, 80640/103097*r68, 3072/103097, -288/103097, 0, 0;93255631/8041092-10886400/670091*r78+725760/670091*r67, -36069450/670091+50803200/670091*r78-3628800/670091*r67, 64346994/670091-91445760/670091*r78+7257600/670091*r67, -158096578/2010273+76204800/670091*r78-7257600/670091*r67, 67080435/2680364-25401600/670091*r78+3628800/670091*r67, -725760/670091*r67, 0, 725760/670091*r78, -145152/670091, 27648/670091, -2592/670091, 0;-3921999/1079524+25401600/5127739*r78+5080320/5127739*r68, 536324953/30766434-121927680/5127739*r78-25401600/5127739*r68, -334792485/10255478+228614400/5127739*r78+50803200/5127739*r68, 296462325/10255478-203212800/5127739*r78-50803200/5127739*r68, -657798011/61532868+76204800/5127739*r78+25401600/5127739*r68, -5080320/5127739*r68, -5080320/5127739*r78, 0, 4064256/5127739, -1016064/5127739, 193536/5127739, -18144/5127739]; - D1(m-7:m,m-11:m)=flipud( fliplr(-D1(1:8,1:12))); + D1(m-7:m,m-11:m)=rot90( -D1(1:8,1:12),2); D1=D1/h; @@ -72,20 +74,26 @@ % -35,-446,980,-980,2450/3,-490,196,-140/3,5, % 5,-80,-266,560,-350,560/3,-70,16,-5/3, % -5/3,20,-140,-126,350,-140,140/3,-10,1]; - D(m-7:m,m-11:m)=flipud( fliplr(-D1(1:8,1:12))); +% D(m-7:m,m-11:m)=rot90( -D1(1:8,1:12),2); + +% D2=(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0)); - D2=(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0)); + diags = -4:4; + left_stencil = [-1/560,8/315,-1/5,8/5]; + stencil = [left_stencil,-205/72,fliplr(left_stencil)]; + D2 = stripeMatrix(stencil, diags, m); + D2(1:8,1:12)=[4870382994799/1358976868290, -893640087518/75498714905,926594825119/60398971924, -1315109406200/135897686829,39126983272/15099742981, 12344491342/75498714905, -451560522577/2717953736580, 0, 0, 0, 0, 0;333806012194/390619153855, -154646272029/111605472530, 1168338040/33481641759, 82699112501/133926567036, -171562838/11160547253, -28244698346/167408208795, 11904122576/167408208795, -2598164715/312495323084, 0, 0, 0, 0;7838984095/52731029988, 1168338040/5649753213, -88747895/144865467, 423587231/627750357, -43205598281/22599012852, 4876378562/1883251071, -5124426509/3766502142, 10496900965/39548272491, 0, 0, 0, 0;-94978241528/828644350023, 82699112501/157837019052, 1270761693/13153084921, -167389605005/118377764289, 48242560214/39459254763, -31673996013/52612339684, 43556319241/118377764289, -44430275135/552429566682, 0, 0, 0, 0;1455067816/21132528431, -171562838/3018932633, -43205598281/36227191596, 48242560214/9056797899, -52276055645/6037865266, 57521587238/9056797899, -80321706377/36227191596, 8078087158/21132528431, -1296/299527, 0, 0, 0;10881504334/327321118845, -28244698346/140280479505, 4876378562/9352031967, -10557998671/12469375956, 57521587238/28056095901, -278531401019/93520319670, 73790130002/46760159835, -137529995233/785570685228, 2048/103097, -144/103097, 0, 0;-135555328849/8509847458140, 11904122576/101307707835, -5124426509/13507694378, 43556319241/60784624701, -80321706377/81046166268, 73790130002/33769235945, -950494905688/303923123505, 239073018673/141830790969, -145152/670091, 18432/670091, -1296/670091, 0;0, -2598164715/206729925524, 10496900965/155047444143, -44430275135/310094888286, 425162482/2720130599, -137529995233/620189776572, 239073018673/155047444143, -144648000000/51682481381, 8128512/5127739, -1016064/5127739, 129024/5127739, -9072/5127739]; - D2(m-7:m,m-11:m)=flipud( fliplr(D2(1:8,1:12) ) ); + D2(m-7:m,m-11:m)=rot90( D2(1:8,1:12) ,2 ); D2=D2/h^2; - DS=zeros(m,m); + DS=sparse(m,m); DS(1,1:7)=-[-4723/2100, 839/175, -157/35, 278/105, -103/140, -1/175, 6/175]; DS(m,m-6:m)=fliplr(-[-4723/2100, 839/175, -157/35, 278/105, -103/140, -1/175, 6/175]); @@ -104,13 +112,13 @@ %xlabel('Real part'); %ylabel('Imaginary part'); %title('Spectrum, minimal spectral radius'); - e_1 = zeros(m,1); + e_1 = sparse(m,1); e_1(1)= 1; - e_m = zeros(m,1); + e_m = sparse(m,1); e_m(end)= 1; S_1 = -DS(1,:)'; S_m = DS(end,:)'; - Q = H*D1-(-e_1*e_1' + e_m*e_m'); + Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_blocknorm_10.m --- a/+sbp/+implementations/d2_blocknorm_10.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_blocknorm_10.m Fri Sep 09 14:53:41 2016 +0200 @@ -8,56 +8,64 @@ H_U=[0.428081020217e12 / 0.2633637888000e13 0.779032713983e12 / 0.2633637888000e13 -0.1187642619571e13 / 0.2633637888000e13 0.1642279196603e13 / 0.2633637888000e13 -0.339289243121e12 / 0.526727577600e12 0.1261055176253e13 / 0.2633637888000e13 -0.658216413073e12 / 0.2633637888000e13 0.33968779823e11 / 0.376233984000e12 -0.764998e6 / 0.40186125e8; 0.779032713983e12 / 0.2633637888000e13 0.317907052061e12 / 0.164602368000e12 -0.1082918052397e13 / 0.658409472000e12 0.176473501369e12 / 0.82301184000e11 -0.521625191587e12 / 0.263363788800e12 0.200523313337e12 / 0.164602368000e12 -0.279496000009e12 / 0.658409472000e12 -0.832e3 / 0.40186125e8 0.129794760887e12 / 0.2633637888000e13; -0.1187642619571e13 / 0.2633637888000e13 -0.1082918052397e13 / 0.658409472000e12 0.1304108863849e13 / 0.329204736000e12 -0.533093695961e12 / 0.131681894400e12 0.213949854133e12 / 0.52672757760e11 -0.1883976009151e13 / 0.658409472000e12 0.51084128e8 / 0.40186125e8 -0.92763684343e11 / 0.658409472000e12 -0.59058717923e11 / 0.526727577600e12; 0.1642279196603e13 / 0.2633637888000e13 0.176473501369e12 / 0.82301184000e11 -0.533093695961e12 / 0.131681894400e12 0.217677310051e12 / 0.32920473600e11 -0.1509120465127e13 / 0.263363788800e12 0.170839232e9 / 0.40186125e8 -0.1404096707137e13 / 0.658409472000e12 0.11789520859e11 / 0.32920473600e11 0.85652315431e11 / 0.526727577600e12; -0.339289243121e12 / 0.526727577600e12 -0.521625191587e12 / 0.263363788800e12 0.213949854133e12 / 0.52672757760e11 -0.1509120465127e13 / 0.263363788800e12 0.21849109e8 / 0.3214890e7 -0.1134422468377e13 / 0.263363788800e12 0.602448430967e12 / 0.263363788800e12 -0.4910542309e10 / 0.10534551552e11 -0.83039945231e11 / 0.526727577600e12; 0.1261055176253e13 / 0.2633637888000e13 0.200523313337e12 / 0.164602368000e12 -0.1883976009151e13 / 0.658409472000e12 0.170839232e9 / 0.40186125e8 -0.1134422468377e13 / 0.263363788800e12 0.681437038097e12 / 0.164602368000e12 -0.1108257453763e13 / 0.658409472000e12 0.31631872327e11 / 0.82301184000e11 0.37820115539e11 / 0.376233984000e12; -0.658216413073e12 / 0.2633637888000e13 -0.279496000009e12 / 0.658409472000e12 0.51084128e8 / 0.40186125e8 -0.1404096707137e13 / 0.658409472000e12 0.602448430967e12 / 0.263363788800e12 -0.1108257453763e13 / 0.658409472000e12 0.623491124887e12 / 0.329204736000e12 -0.146643738067e12 / 0.658409472000e12 -0.98874149197e11 / 0.2633637888000e13; 0.33968779823e11 / 0.376233984000e12 -0.832e3 / 0.40186125e8 -0.92763684343e11 / 0.658409472000e12 0.11789520859e11 / 0.32920473600e11 -0.4910542309e10 / 0.10534551552e11 0.31631872327e11 / 0.82301184000e11 -0.146643738067e12 / 0.658409472000e12 0.174599973347e12 / 0.164602368000e12 0.4625165773e10 / 0.526727577600e12; -0.764998e6 / 0.40186125e8 0.129794760887e12 / 0.2633637888000e13 -0.59058717923e11 / 0.526727577600e12 0.85652315431e11 / 0.526727577600e12 -0.83039945231e11 / 0.526727577600e12 0.37820115539e11 / 0.376233984000e12 -0.98874149197e11 / 0.2633637888000e13 0.4625165773e10 / 0.526727577600e12 0.525286231387e12 / 0.526727577600e12;]; - H=eye(m); + H=speye(m); H(1:9,1:9)=H_U; - H(m-8:m,m-8:m)=flipud( fliplr(H_U(1:9,1:9) ) ); + H(m-8:m,m-8:m)=rot90( H_U(1:9,1:9) ,2 ); H=H*h; HI=inv(H); - Q=1/1260*diag(ones(m-5,1),5)-5/504*diag(ones(m-4,1),4)+5/84*diag(ones(m-3,1),3)-5/21*diag(ones(m-2,1),2)+5/6*diag(ones(m-1,1),1)-5/6*diag(ones(m-1,1),-1)+5/21*diag(ones(m-2,1),-2)-5/84*diag(ones(m-3,1),-3)+5/504*diag(ones(m-4,1),-4)-1/1260*diag(ones(m-5,1),-5) ; +% Q=1/1260*diag(ones(m-5,1),5)-5/504*diag(ones(m-4,1),4)+5/84*diag(ones(m-3,1),3)-5/21*diag(ones(m-2,1),2)+5/6*diag(ones(m-1,1),1)-5/6*diag(ones(m-1,1),-1)+5/21*diag(ones(m-2,1),-2)-5/84*diag(ones(m-3,1),-3)+5/504*diag(ones(m-4,1),-4)-1/1260*diag(ones(m-5,1),-5) ; + diags = -5:5; + stencil = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; + Q = stripeMatrix(stencil, diags, m); Q_U = [-0.1e1 / 0.2e1 0.78249683e8 / 0.71442000e8 -0.28290472447e11 / 0.18289152000e11 0.4285528063e10 / 0.2032128000e10 -0.3924872557e10 / 0.1828915200e10 0.2856344621e10 / 0.1828915200e10 -0.1598284927e10 / 0.2032128000e10 0.4652402687e10 / 0.18289152000e11 -0.731623e6 / 0.17860500e8; -0.78249683e8 / 0.71442000e8 0 0.5351945471e10 / 0.2032128000e10 -0.4078428731e10 / 0.1143072000e10 0.4879509877e10 / 0.1219276800e10 -0.1273709579e10 / 0.381024000e9 0.7433128649e10 / 0.3657830400e10 -0.1239008e7 / 0.1488375e7 0.21251483e8 / 0.124416000e9; 0.28290472447e11 / 0.18289152000e11 -0.5351945471e10 / 0.2032128000e10 0 0.6701592799e10 / 0.2612736000e10 -0.57535927e8 / 0.17418240e8 0.3058732543e10 / 0.870912000e9 -0.28821953e8 / 0.10206000e8 0.1858901437e10 / 0.1219276800e10 -0.2363118211e10 / 0.6096384000e10; -0.4285528063e10 / 0.2032128000e10 0.4078428731e10 / 0.1143072000e10 -0.6701592799e10 / 0.2612736000e10 0 0.52016695e8 / 0.20901888e8 -0.3777923e7 / 0.1275750e7 0.123333949e9 / 0.41472000e8 -0.2256321727e10 / 0.1143072000e10 0.1058261459e10 / 0.1828915200e10; 0.3924872557e10 / 0.1828915200e10 -0.4879509877e10 / 0.1219276800e10 0.57535927e8 / 0.17418240e8 -0.52016695e8 / 0.20901888e8 0 0.125390297e9 / 0.58060800e8 -0.609569351e9 / 0.261273600e9 0.438488399e9 / 0.243855360e9 -0.14244569e8 / 0.24385536e8; -0.2856344621e10 / 0.1828915200e10 0.1273709579e10 / 0.381024000e9 -0.3058732543e10 / 0.870912000e9 0.3777923e7 / 0.1275750e7 -0.125390297e9 / 0.58060800e8 0 0.4624381729e10 / 0.2612736000e10 -0.484096919e9 / 0.381024000e9 0.2676438019e10 / 0.6096384000e10; 0.1598284927e10 / 0.2032128000e10 -0.7433128649e10 / 0.3657830400e10 0.28821953e8 / 0.10206000e8 -0.123333949e9 / 0.41472000e8 0.609569351e9 / 0.261273600e9 -0.4624381729e10 / 0.2612736000e10 0 0.21500967689e11 / 0.18289152000e11 -0.7199454721e10 / 0.18289152000e11; -0.4652402687e10 / 0.18289152000e11 0.1239008e7 / 0.1488375e7 -0.1858901437e10 / 0.1219276800e10 0.2256321727e10 / 0.1143072000e10 -0.438488399e9 / 0.243855360e9 0.484096919e9 / 0.381024000e9 -0.21500967689e11 / 0.18289152000e11 0 0.761653e6 / 0.882000e6; 0.731623e6 / 0.17860500e8 -0.21251483e8 / 0.124416000e9 0.2363118211e10 / 0.6096384000e10 -0.1058261459e10 / 0.1828915200e10 0.14244569e8 / 0.24385536e8 -0.2676438019e10 / 0.6096384000e10 0.7199454721e10 / 0.18289152000e11 -0.761653e6 / 0.882000e6 0;]; Q(1:9,1:9)=Q_U; - Q(m-8:m,m-8:m)=flipud( fliplr(-Q_U(1:9,1:9) ) ); + Q(m-8:m,m-8:m)=rot90( -Q_U(1:9,1:9) ,2 ); - D1=HI*Q; + D1=H\Q; M_U =[0.3812926003e10 / 0.2438553600e10 -0.5433856529e10 / 0.1741824000e10 0.4187462879e10 / 0.1045094400e10 -0.65635105447e11 / 0.12192768000e11 0.1457682577e10 / 0.270950400e9 -0.27884016067e11 / 0.7315660800e10 0.22304839493e11 / 0.12192768000e11 -0.188132543e9 / 0.348364800e9 0.42711619e8 / 0.571536000e9; -0.5433885329e10 / 0.1741824000e10 0.23985229969e11 / 0.2286144000e10 -0.10208460799e11 / 0.571536000e9 0.8828370001e10 / 0.381024000e9 -0.12306735263e11 / 0.522547200e9 0.39313626089e11 / 0.2286144000e10 -0.4380200287e10 / 0.508032000e9 0.192498023e9 / 0.71442000e8 -0.14565232681e11 / 0.36578304000e11; 0.29313778073e11 / 0.7315660800e10 -0.10209211399e11 / 0.571536000e9 0.24157533391e11 / 0.653184000e9 -0.6561725111e10 / 0.130636800e9 0.26490755639e11 / 0.522547200e9 -0.352801289e9 / 0.9331200e7 0.401374423e9 / 0.20412000e8 -0.29541854057e11 / 0.4572288000e10 0.7396364989e10 / 0.7315660800e10; -0.65653163047e11 / 0.12192768000e11 0.8832999601e10 / 0.381024000e9 -0.6566554871e10 / 0.130636800e9 0.1575758731e10 / 0.21772800e8 -0.2571648133e10 / 0.34836480e8 0.556969019e9 / 0.10206000e8 -0.3139265911e10 / 0.108864000e9 0.165424529e9 / 0.16934400e8 -0.11623567549e11 / 0.7315660800e10; 0.1458632977e10 / 0.270950400e9 -0.86260417241e11 / 0.3657830400e10 0.3792749777e10 / 0.74649600e8 -0.2576699653e10 / 0.34836480e8 0.157840723e9 / 0.2041200e7 -0.29838889141e11 / 0.522547200e9 0.5142742211e10 / 0.174182400e9 -0.36792522023e11 / 0.3657830400e10 0.2438136689e10 / 0.1463132160e10; -0.27909676867e11 / 0.7315660800e10 0.39389053289e11 / 0.2286144000e10 -0.2478149663e10 / 0.65318400e8 0.559214069e9 / 0.10206000e8 -0.29914661941e11 / 0.522547200e9 0.2843819551e10 / 0.65318400e8 -0.14816015149e11 / 0.653184000e9 0.1677001673e10 / 0.228614400e9 -0.44441740171e11 / 0.36578304000e11; 0.3188985299e10 / 0.1741824000e10 -0.626864041e9 / 0.72576000e8 0.57539389e8 / 0.2916000e7 -0.3153500311e10 / 0.108864000e9 0.5162239811e10 / 0.174182400e9 -0.14840163949e11 / 0.653184000e9 0.419025709e9 / 0.31104000e8 -0.138945749e9 / 0.27216000e8 0.591880819e9 / 0.746496000e9; -0.1317440441e10 / 0.2438553600e10 0.192696473e9 / 0.71442000e8 -0.4230355151e10 / 0.653184000e9 0.165983249e9 / 0.16934400e8 -0.36905792423e11 / 0.3657830400e10 0.1679779433e10 / 0.228614400e9 -0.972870443e9 / 0.190512000e9 0.9129544111e10 / 0.2286144000e10 -0.13387742111e11 / 0.7315660800e10; 0.42721069e8 / 0.571536000e9 -0.14572922281e11 / 0.36578304000e11 0.7407199549e10 / 0.7315660800e10 -0.11649228349e11 / 0.7315660800e10 0.349038407e9 / 0.209018880e9 -0.44495912971e11 / 0.36578304000e11 0.29009849731e11 / 0.36578304000e11 -0.13387863071e11 / 0.7315660800e10 0.21585797479e11 / 0.7315660800e10;]; - T=[-0.1e1 / 0.3150e4 0.5e1 / 0.1008e4 -0.5e1 / 0.126e3 0.5e1 / 0.21e2 -0.5e1 / 0.3e1 0.5269e4 / 0.1800e4 -0.5e1 / 0.3e1 0.5e1 / 0.21e2 -0.5e1 / 0.126e3 0.5e1 / 0.1008e4 -0.1e1 / 0.3150e4;]; - M=(T(1)*diag(ones(m-5,1),5)+T(2)*diag(ones(m-4,1),4)+T(3)*diag(ones(m-3,1),3)+T(4)*diag(ones(m-2,1),2)+T(5)*diag(ones(m-1,1),1)+T(7)*diag(ones(m-1,1),-1)+T(8)*diag(ones(m-2,1),-2)+T(9)*diag(ones(m-3,1),-3)+T(10)*diag(ones(m-4,1),-4)+T(11)*diag(ones(m-5,1),-5)+T(6)*diag(ones(m,1),0)); +% T=[-0.1e1 / 0.3150e4 0.5e1 / 0.1008e4 -0.5e1 / 0.126e3 0.5e1 / 0.21e2 -0.5e1 / 0.3e1 0.5269e4 / 0.1800e4 -0.5e1 / 0.3e1 0.5e1 / 0.21e2 -0.5e1 / 0.126e3 0.5e1 / 0.1008e4 -0.1e1 / 0.3150e4;]; +% M=(T(1)*diag(ones(m-5,1),5)+T(2)*diag(ones(m-4,1),4)+T(3)*diag(ones(m-3,1),3)+T(4)*diag(ones(m-2,1),2)+T(5)*diag(ones(m-1,1),1)+T(7)*diag(ones(m-1,1),-1)+T(8)*diag(ones(m-2,1),-2)+T(9)*diag(ones(m-3,1),-3)+T(10)*diag(ones(m-4,1),-4)+T(11)*diag(ones(m-5,1),-5)+T(6)*diag(ones(m,1),0)); + + diags = -5:5; + left_stencil = [-1/3150,5/1008,-5/126,5/21,-5/3]; + stencil = [left_stencil,5269/1800,fliplr(left_stencil)]; + M = stripeMatrix(stencil, diags, m); M(1:9,1:9)=M_U; - M(m-8:m,m-8:m)=flipud( fliplr( M_U ) ); + M(m-8:m,m-8:m)=rot90( M_U ,2 ); M=M/h; DS_U=[0.761e3 / 0.280e3 -8 14 -0.56e2 / 0.3e1 0.35e2 / 0.2e1 -0.56e2 / 0.5e1 0.14e2 / 0.3e1 -0.8e1 / 0.7e1 0.1e1 / 0.8e1;]; - DS=zeros(m,m); + DS=sparse(m,m); DS(1,1:9)=DS_U; DS(m,m-8:m)=fliplr(DS_U); DS=DS/h; - D2=HI*(-M+DS); + D2=H\(-M+DS); % Try adding AD to boundary D1=HI*(Q-D9'*D9) - DD_9=zeros(m); - d9=[-1 9 -36 84 -126 126 -84 36 -9 1];t9=sum(abs(d9));%d9=d9/t9; - DD_9(1:1,1:10)=[d9]; - DD_9(m:m,m-9:m)=[d9]; - +% DD_9=sparse(m); +% d9=[-1 9 -36 84 -126 126 -84 36 -9 1];t9=sum(abs(d9));%d9=d9/t9; +% DD_9(1:1,1:10)=[d9]; +% DD_9(m:m,m-9:m)=[d9]; +% +% +% ADD=30*h/(t9)*DD_9'*DD_9; - ADD=30*h/(t9)*DD_9'*DD_9; - - e_1 = zeros(m,1); + e_1 = sparse(m,1); e_1(1)= 1; - e_m = zeros(m,1); + e_m = sparse(m,1); e_m(end)= 1; S_1 = -DS(1,:)'; S_m = DS(end,:)'; - Q = H*D1-(-e_1*e_1' + e_m*e_m'); + Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_blocknorm_4.m --- a/+sbp/+implementations/d2_blocknorm_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_blocknorm_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -8,56 +8,59 @@ H_U=[0.751e3 / 0.3456e4 0.661e3 / 0.3456e4 -0.515e3 / 0.3456e4 0.5e1 / 0.128e3; 0.661e3 / 0.3456e4 0.1405e4 / 0.1152e4 -0.3e1 / 0.128e3 0.29e2 / 0.3456e4; -0.515e3 / 0.3456e4 -0.3e1 / 0.128e3 0.989e3 / 0.1152e4 0.149e3 / 0.3456e4; 0.5e1 / 0.128e3 0.29e2 / 0.3456e4 0.149e3 / 0.3456e4 0.3407e4 / 0.3456e4;]; - H=eye(m); + H=speye(m); H(1:4,1:4)=H_U; - H(m-3:m,m-3:m)=flipud( fliplr(H_U(1:4,1:4) ) ); + H(m-3:m,m-3:m)=rot90( H_U(1:4,1:4) ,2 ); H=H*h; HI=inv(H); - Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); + e=ones(m,1); + Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; +% Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); Q_U = [-0.1e1 / 0.2e1 0.55e2 / 0.72e2 -0.47e2 / 0.144e3 0.1e1 / 0.16e2; -0.55e2 / 0.72e2 0 0.43e2 / 0.48e2 -0.19e2 / 0.144e3; 0.47e2 / 0.144e3 -0.43e2 / 0.48e2 0 0.47e2 / 0.72e2; -0.1e1 / 0.16e2 0.19e2 / 0.144e3 -0.47e2 / 0.72e2 0;]; Q(1:4,1:4)=Q_U; - Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) ); + Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); - D1=HI*Q; + D1=H\Q; M_U=[0.359e3 / 0.288e3 -0.443e3 / 0.288e3 0.97e2 / 0.288e3 -0.13e2 / 0.288e3; -0.51e2 / 0.32e2 0.325e3 / 0.96e2 -0.191e3 / 0.96e2 0.19e2 / 0.96e2; 0.43e2 / 0.96e2 -0.69e2 / 0.32e2 0.293e3 / 0.96e2 -0.137e3 / 0.96e2; -0.29e2 / 0.288e3 0.89e2 / 0.288e3 -0.427e3 / 0.288e3 0.727e3 / 0.288e3;]; - M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); +% M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); + M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; M(1:4,1:4)=M_U; - M(m-3:m,m-3:m)=flipud( fliplr( M_U ) ); + M(m-3:m,m-3:m)=rot90( M_U ,2 ); M=M/h; DS_U=[0.25e2 / 0.12e2 -4 3 -0.4e1 / 0.3e1 0.1e1 / 0.4e1;]; - DS=zeros(m,m); + DS=sparse(m,m); DS(1,1:5)=DS_U; DS(m,m-4:m)=fliplr(DS_U); DS=DS/h; - D2=HI*(-M+DS); + D2=H\(-M+DS); - d3=[-1 3 -3 1]; - t3=sum(abs(d3)); - DD_3(1:1,1:4)=[d3]; - DD_3(m:m,m-3:m)=[d3]; +% d3=[-1 3 -3 1]; +% t3=sum(abs(d3)); +% DD_3(1:1,1:4)=[d3]; +% DD_3(m:m,m-3:m)=[d3]; % This works for wave eq. % For studs interface in 1D no AD is needed. - ADD=1*h/(t3)*DD_3'*DD_3; +% ADD=1*h/(t3)*DD_3'*DD_3; - e_1 = zeros(m,1); + e_1 = sparse(m,1); e_1(1)= 1; - e_m = zeros(m,1); + e_m = sparse(m,1); e_m(end)= 1; S_1 = -DS(1,:)'; S_m = DS(end,:)'; - Q = H*D1-(-e_1*e_1' + e_m*e_m'); + Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_blocknorm_6.m --- a/+sbp/+implementations/d2_blocknorm_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_blocknorm_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -8,56 +8,63 @@ H_U=[0.8489084265971e13 / 0.45952647390720e14 0.24636450459943e14 / 0.98469958694400e14 -0.2796787072531e13 / 0.12308744836800e14 0.2793599068823e13 / 0.14360202309600e14 -0.66344569931569e14 / 0.689289710860800e15 0.3784697867191e13 / 0.137857942172160e15; 0.24636450459943e14 / 0.98469958694400e14 0.27815394775103e14 / 0.19693991738880e14 -0.445601472229e12 / 0.861612138576e12 0.3896159037731e13 / 0.17232242771520e14 -0.866505556741e12 / 0.27571588434432e14 -0.25625418493681e14 / 0.689289710860800e15; -0.2796787072531e13 / 0.12308744836800e14 -0.445601472229e12 / 0.861612138576e12 0.31409405327129e14 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.2651608170899e13 / 0.17232242771520e14 0.1434714163381e13 / 0.43080606928800e14; 0.2793599068823e13 / 0.14360202309600e14 0.3896159037731e13 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.6984350202787e13 / 0.5744080923840e13 -0.62662743973e11 / 0.861612138576e12 -0.435331581619e12 / 0.12308744836800e14; -0.66344569931569e14 / 0.689289710860800e15 -0.866505556741e12 / 0.27571588434432e14 0.2651608170899e13 / 0.17232242771520e14 -0.62662743973e11 / 0.861612138576e12 0.20320736807807e14 / 0.19693991738880e14 0.1368363924007e13 / 0.98469958694400e14; 0.3784697867191e13 / 0.137857942172160e15 -0.25625418493681e14 / 0.689289710860800e15 0.1434714163381e13 / 0.43080606928800e14 -0.435331581619e12 / 0.12308744836800e14 0.1368363924007e13 / 0.98469958694400e14 0.27414523542149e14 / 0.27571588434432e14;]; - H=eye(m); + H=speye(m); H(1:6,1:6)=H_U; - H(m-5:m,m-5:m)=flipud( fliplr(H_U(1:6,1:6) ) ); + H(m-5:m,m-5:m)=rot90( H_U(1:6,1:6) ,2 ); H=H*h; HI=inv(H); - Q=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); + diags = -3:3; + stencil = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; + Q = stripeMatrix(stencil, diags, m); + +% Q=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); Q_U = [-0.1e1 / 0.2e1 0.151864337282617e15 / 0.172322427715200e15 -0.251539972254817e15 / 0.344644855430400e15 0.61230525943549e14 / 0.114881618476800e15 -0.80987306509439e14 / 0.344644855430400e15 0.697178163343e12 / 0.13785794217216e14; -0.151864337282617e15 / 0.172322427715200e15 0 0.12350422095979e14 / 0.7658774565120e13 -0.78802251164141e14 / 0.68928971086080e14 0.4229407848431e13 / 0.7658774565120e13 -0.5372490790279e13 / 0.38293872825600e14; 0.251539972254817e15 / 0.344644855430400e15 -0.12350422095979e14 / 0.7658774565120e13 0 0.2217674201683e13 / 0.1723224277152e13 -0.13219134462287e14 / 0.22976323695360e14 0.19660399553981e14 / 0.114881618476800e15; -0.61230525943549e14 / 0.114881618476800e15 0.78802251164141e14 / 0.68928971086080e14 -0.2217674201683e13 / 0.1723224277152e13 0 0.62307836637379e14 / 0.68928971086080e14 -0.84068101764193e14 / 0.344644855430400e15; 0.80987306509439e14 / 0.344644855430400e15 -0.4229407848431e13 / 0.7658774565120e13 0.13219134462287e14 / 0.22976323695360e14 -0.62307836637379e14 / 0.68928971086080e14 0 0.44756810052211e14 / 0.57440809238400e14; -0.697178163343e12 / 0.13785794217216e14 0.5372490790279e13 / 0.38293872825600e14 -0.19660399553981e14 / 0.114881618476800e15 0.84068101764193e14 / 0.344644855430400e15 -0.44756810052211e14 / 0.57440809238400e14 0;]; Q(1:6,1:6)=Q_U; - Q(m-5:m,m-5:m)=flipud( fliplr(-Q_U(1:6,1:6) ) ); + Q(m-5:m,m-5:m)=rot90( -Q_U(1:6,1:6) ,2 ); - D1=HI*Q; + D1=H\Q; M_U=[0.960901171090739e15 / 0.689289710860800e15 -0.502032138770899e15 / 0.229763236953600e15 0.493085196645929e15 / 0.344644855430400e15 -0.329491854944251e15 / 0.344644855430400e15 0.89541920186441e14 / 0.229763236953600e15 -0.50617198740721e14 / 0.689289710860800e15; -0.100483015499831e15 / 0.45952647390720e14 0.807564929223191e15 / 0.137857942172160e15 -0.415779274818991e15 / 0.68928971086080e14 0.80693719872887e14 / 0.22976323695360e14 -0.196663473955997e15 / 0.137857942172160e15 0.37943821632959e14 / 0.137857942172160e15; 0.99938177941669e14 / 0.68928971086080e14 -0.84419552767043e14 / 0.13785794217216e14 0.106922123424097e15 / 0.11488161847680e14 -0.223356054245897e15 / 0.34464485543040e14 0.157526160982357e15 / 0.68928971086080e14 -0.10062402380533e14 / 0.22976323695360e14; -0.68310884976863e14 / 0.68928971086080e14 0.17038649985979e14 / 0.4595264739072e13 -0.231397767539273e15 / 0.34464485543040e14 0.232669188399619e15 / 0.34464485543040e14 -0.1657930371065e13 / 0.510584971008e12 0.34774771016773e14 / 0.68928971086080e14; 0.18789143112277e14 / 0.45952647390720e14 -0.213895716727517e15 / 0.137857942172160e15 0.171024751153381e15 / 0.68928971086080e14 -0.8523669967037e13 / 0.2552924855040e13 0.485768751245399e15 / 0.137857942172160e15 -0.229158724354277e15 / 0.137857942172160e15; -0.51766014925489e14 / 0.689289710860800e15 0.202930494289627e15 / 0.689289710860800e15 -0.54332868549353e14 / 0.114881618476800e15 0.180479548146281e15 / 0.344644855430400e15 -0.1146942437956153e16 / 0.689289710860800e15 0.211001773091419e15 / 0.76587745651200e14;]; - M=-(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180; +% M=-(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180; + diags = -3:3; + stencil = -1/180*[2,-27,270,-490,270,-27,2]; + M = stripeMatrix(stencil, diags, m); M(1:6,1:6)=M_U; - M(m-5:m,m-5:m)=flipud( fliplr( M_U ) ); + M(m-5:m,m-5:m)=rot90( M_U ,2 ); M=M/h; DS_U=[0.137e3 / 0.60e2 -5 5 -0.10e2 / 0.3e1 0.5e1 / 0.4e1 -0.1e1 / 0.5e1;]; - DS=zeros(m,m); + DS=sparse(m,m); DS(1,1:6)=DS_U; DS(m,m-5:m)=fliplr(DS_U); DS=DS/h; - D2=HI*(-M+DS); - - d5=[-1 5 -10 10 -5 1]; - t5=sum(abs(d5)); - DD_5(1:1,1:6)=[d5]; - DD_5(m:m,m-5:m)=[d5]; + D2=H\(-M+DS); - % This works for wave eq. - % For studs interface in 1D no AD is needed. - ADD=7*h/(t5)*DD_5'*DD_5; +% d5=[-1 5 -10 10 -5 1]; +% t5=sum(abs(d5)); +% DD_5(1:1,1:6)=[d5]; +% DD_5(m:m,m-5:m)=[d5]; +% +% % This works for wave eq. +% % For studs interface in 1D no AD is needed. +% ADD=7*h/(t5)*DD_5'*DD_5; - e_1 = zeros(m,1); + e_1 = sparse(m,1); e_1(1)= 1; - e_m = zeros(m,1); + e_m = sparse(m,1); e_m(end)= 1; S_1 = -DS(1,:)'; S_m = DS(end,:)'; - Q = H*D1-(-e_1*e_1' + e_m*e_m'); + Q = H*D1-(-(e_1*e_1') + (e_m*e_m')); M = -(H*D2-(-e_1*S_1' + e_m*S_m')); end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_blocknorm_8.m --- a/+sbp/+implementations/d2_blocknorm_8.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_blocknorm_8.m Fri Sep 09 14:53:41 2016 +0200 @@ -6,8 +6,8 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - 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; H_U=[ 0.704266523e9 / 0.4180377600e10 0.4579586639e10 / 0.16257024000e11 -0.3623870581e10 / 0.9754214400e10 0.12753127559e11 / 0.29262643200e11 -0.3687413731e10 / 0.9754214400e10 0.2169892891e10 / 0.9754214400e10 -0.13224544841e11 / 0.146313216000e12 0.4142047e7 / 0.199065600e9; 0.4579586639e10 / 0.16257024000e11 0.36543258551e11 / 0.20901888000e11 -0.8235820121e10 / 0.6967296000e10 0.1800520829e10 / 0.1393459200e10 -0.3725834681e10 / 0.4180377600e10 0.2588501879e10 / 0.6967296000e10 -0.10477621e8 / 0.995328000e9 -0.6589395529e10 / 0.146313216000e12; @@ -20,33 +20,40 @@ ]; - H=eye(m); + H=speye(m); H(1:8,1:8)=H_U; H(m-7:m,m-7:m)=rot90( H_U(1:8,1:8) ,2 ); H=H*h; HI=inv(H); - Q=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4)); - +% Q=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4)); + diags = -4:4; + stencil = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; + Q = stripeMatrix(stencil, diags, m); Q_U = [-0.1e1 / 0.2e1 0.16262381e8 / 0.15876000e8 -0.3770744693e10 / 0.3048192000e10 0.290431859e9 / 0.203212800e9 -0.363704879e9 / 0.304819200e9 0.206906927e9 / 0.304819200e9 -0.248883679e9 / 0.1016064000e10 0.2665637e7 / 0.62208000e8; -0.16262381e8 / 0.15876000e8 0 0.138589901e9 / 0.62208000e8 -0.31764881e8 / 0.12441600e8 0.461249e6 / 0.193536e6 -0.347232997e9 / 0.217728000e9 0.1788157e7 / 0.2488320e7 -0.5926349e7 / 0.37632000e8; 0.3770744693e10 / 0.3048192000e10 -0.138589901e9 / 0.62208000e8 0 0.11741773e8 / 0.5443200e7 -0.39109817e8 / 0.17418240e8 0.10216441e8 / 0.5376000e7 -0.245131109e9 / 0.217728000e9 0.92809903e8 / 0.304819200e9; -0.290431859e9 / 0.203212800e9 0.31764881e8 / 0.12441600e8 -0.11741773e8 / 0.5443200e7 0 0.4634999e7 / 0.2488320e7 -0.144219869e9 / 0.87091200e8 0.17445643e8 / 0.14515200e8 -0.38142949e8 / 0.101606400e9; 0.363704879e9 / 0.304819200e9 -0.461249e6 / 0.193536e6 0.39109817e8 / 0.17418240e8 -0.4634999e7 / 0.2488320e7 0 0.7992221e7 / 0.5443200e7 -0.817951e6 / 0.829440e6 0.4455517e7 / 0.13547520e8; -0.206906927e9 / 0.304819200e9 0.347232997e9 / 0.217728000e9 -0.10216441e8 / 0.5376000e7 0.144219869e9 / 0.87091200e8 -0.7992221e7 / 0.5443200e7 0 0.68487373e8 / 0.62208000e8 -0.1032638773e10 / 0.3048192000e10; 0.248883679e9 / 0.1016064000e10 -0.1788157e7 / 0.2488320e7 0.245131109e9 / 0.217728000e9 -0.17445643e8 / 0.14515200e8 0.817951e6 / 0.829440e6 -0.68487373e8 / 0.62208000e8 0 0.39529771e8 / 0.47628000e8; -0.2665637e7 / 0.62208000e8 0.5926349e7 / 0.37632000e8 -0.92809903e8 / 0.304819200e9 0.38142949e8 / 0.101606400e9 -0.4455517e7 / 0.13547520e8 0.1032638773e10 / 0.3048192000e10 -0.39529771e8 / 0.47628000e8 0;]; Q(1:8,1:8)=Q_U; Q(m-7:m,m-7:m)=rot90( -Q_U(1:8,1:8) ,2 ); - D1=HI*Q; + D1=H\Q; M_U =[0.27667249117e11 / 0.18289152000e11 -0.17100791927e11 / 0.6096384000e10 0.6123596021e10 / 0.2032128000e10 -0.12420079921e11 / 0.3657830400e10 0.3352522937e10 / 0.1219276800e10 -0.3030351383e10 / 0.2032128000e10 0.8955233071e10 / 0.18289152000e11 -0.448917533e9 / 0.6096384000e10; -0.2443029521e10 / 0.870912000e9 0.3279926909e10 / 0.373248000e9 -0.150833107e9 / 0.11612160e8 0.2418903029e10 / 0.174182400e9 -0.1195687489e10 / 0.104509440e9 0.1864443097e10 / 0.290304000e9 -0.275412413e9 / 0.124416000e9 0.26267539e8 / 0.74649600e8; 0.875033123e9 / 0.290304000e9 -0.754432799e9 / 0.58060800e8 0.262316881e9 / 0.10752000e8 -0.1615952663e10 / 0.58060800e8 0.1304948581e10 / 0.58060800e8 -0.59605951e8 / 0.4608000e7 0.270029509e9 / 0.58060800e8 -0.225377137e9 / 0.290304000e9; -0.71086111e8 / 0.20901888e8 0.2425994741e10 / 0.174182400e9 -0.1622486807e10 / 0.58060800e8 0.735382895e9 / 0.20901888e8 -0.1016121419e10 / 0.34836480e8 0.190014817e9 / 0.11612160e8 -0.447155539e9 / 0.74649600e8 0.179406911e9 / 0.174182400e9; 0.480578879e9 / 0.174182400e9 -0.6018333509e10 / 0.522547200e9 0.1319413093e10 / 0.58060800e8 -0.205032463e9 / 0.6967296e7 0.551889007e9 / 0.20901888e8 -0.887809303e9 / 0.58060800e8 0.914606453e9 / 0.174182400e9 -0.67482881e8 / 0.74649600e8; -0.434493809e9 / 0.290304000e9 0.1878773977e10 / 0.290304000e9 -0.423185977e9 / 0.32256000e8 0.964538597e9 / 0.58060800e8 -0.894343447e9 / 0.58060800e8 0.345461491e9 / 0.32256000e8 -0.1288081307e10 / 0.290304000e9 0.199200163e9 / 0.290304000e9; 0.183060319e9 / 0.373248000e9 -0.276656573e9 / 0.124416000e9 0.54579137e8 / 0.11612160e8 -0.3169984837e10 / 0.522547200e9 0.184339633e9 / 0.34836480e8 -0.1289417627e10 / 0.290304000e9 0.1431981949e10 / 0.373248000e9 -0.307164061e9 / 0.174182400e9; -0.449332253e9 / 0.6096384000e10 0.1290053923e10 / 0.3657830400e10 -0.1588745239e10 / 0.2032128000e10 0.1267377593e10 / 0.1219276800e10 -0.3326650673e10 / 0.3657830400e10 0.1396036981e10 / 0.2032128000e10 -0.2150231371e10 / 0.1219276800e10 0.52544801501e11 / 0.18289152000e11;]; - M=-(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0)); +% M=-(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0)); + diags = -4:4; + left_stencil = [-1/560,8/315,-1/5,8/5]; + stencil = [left_stencil,-205/72,fliplr(left_stencil)]; + M = stripeMatrix(stencil, diags, m); + M(1:8,1:8)=M_U; M(m-7:m,m-7:m)=rot90( M_U ,2 ); M=M/h; % DS_U=[0.363e3 / 0.140e3 -7 0.21e2 / 0.2e1 -0.35e2 / 0.3e1 0.35e2 / 0.4e1 -0.21e2 / 0.5e1 0.7e1 / 0.6e1 -0.1e1 / 0.7e1;]; - % DS=zeros(m,m); + % DS=sparse(m,m); % DS(1,1:8)=DS_U; % DS(m,m-7:m)=fliplr(DS_U); % DS=DS/h; @@ -54,13 +61,13 @@ %D2=HI*(-M+DS); S_U=-[0.363e3 / 0.140e3 -7 0.21e2 / 0.2e1 -0.35e2 / 0.3e1 0.35e2 / 0.4e1 -0.21e2 / 0.5e1 0.7e1 / 0.6e1 -0.1e1 / 0.7e1;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:8)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-7:m)=fliplr(-S_U); - D2=HI*(-M - e_1*S_1+e_m*S_m); + D2=H\(-M - e_1*S_1+e_m*S_m); S_1 = S_1'; S_m = S_m'; end \ No newline at end of file
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d2_variable_4.m --- a/+sbp/+implementations/d2_variable_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d2_variable_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -7,7 +7,7 @@ N = m; - H = speye(N); + H = spspeye(N); H(1,1) = 17/48; H(2,2) = 59/48; H(3,3) = 43/48; H(4,4) = 49/48; H(N,N) = 17/48; H(N-1,N-1) = 59/48; H(N-2,N-2) = 43/48; H(N-3,N-3) = 49/48; H = h*H; @@ -38,7 +38,7 @@ Q(1:4,1:4)=Q_U; Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); - D1=HI*(Q-1/2*e_l*e_l'+1/2*e_r*e_r') ; + D1=HI*(Q-1/2*(e_l*e_l')+1/2*(e_r*e_r')) ; function D2 = D2_fun(c) M = 78+(N-12)*5; @@ -62,7 +62,7 @@ - R = zeros(M,1); + R = sparse(M,1); R(1:24) = reshape(U(1:4,:)',24,1); R(25:30) = U(5,:); R(31) = -c(5+1) / 0.6e1 + c(5) / 0.8e1 + c(5+2) / 0.8e1; @@ -92,10 +92,10 @@ - BS = sparse(N,N); - BS(1,1:4) = -c(1)*1/h*[(-11/6);3;(-3/2);1/3]; - BS(N,N-3:N) = c(N)*1/h*[(-1/3);3/2;(-3);11/6]; - BS = sparse(BS); +% BS = sparse(N,N); +% BS(1,1:4) = -c(1)*1/h*[(-11/6);3;(-3/2);1/3]; +% BS(N,N-3:N) = c(N)*1/h*[(-1/3);3/2;(-3);11/6]; +% BS = sparse(BS); % %%Row and column indices%% M = 78+(N-12)*5; @@ -107,7 +107,7 @@ (N-4)*ones(7,1);... kron([N-3;N-2;N-1;N],ones(6,1))]; - cols = zeros(M,1); + cols = sparse(M,1); cols(1:24) = kron(ones(4,1),[1;2;3;4;5;6]); cols(25:39) = [(1:7)';(1:8)']; cols(M-23:M) = kron(ones(4,1),[N-5;N-4;N-3;N-2;N-1;N]);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_4.m --- a/+sbp/+implementations/d4_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -33,55 +33,58 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H_U=[0.35809e5 / 0.100800e6 0 0 0 0 0; 0 0.13297e5 / 0.11200e5 0 0 0 0; 0 0 0.5701e4 / 0.5600e4 0 0 0; 0 0 0 0.45109e5 / 0.50400e5 0 0; 0 0 0 0 0.35191e5 / 0.33600e5 0; 0 0 0 0 0 0.33503e5 / 0.33600e5;]; H(1:6,1:6)=H_U; - H(m-5:m,m-5:m)=fliplr(flipud(H_U)); + H(m-5:m,m-5:m)=rot90(H_U,2); H=H*h; HI=inv(H); % First derivative SBP operator, 1st order accurate at first 6 boundary points - q2=-1/12;q1=8/12; - Q=q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% q2=-1/12;q1=8/12; +% Q=q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + e=ones(m,1); + Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); Q_U = [0 0.526249e6 / 0.907200e6 -0.10819e5 / 0.777600e6 -0.50767e5 / 0.907200e6 -0.631e3 / 0.28800e5 0.91e2 / 0.7776e4; -0.526249e6 / 0.907200e6 0 0.1421209e7 / 0.2721600e7 0.16657e5 / 0.201600e6 -0.8467e4 / 0.453600e6 -0.33059e5 / 0.5443200e7; 0.10819e5 / 0.777600e6 -0.1421209e7 / 0.2721600e7 0 0.631187e6 / 0.1360800e7 0.400139e6 / 0.5443200e7 -0.8789e4 / 0.302400e6; 0.50767e5 / 0.907200e6 -0.16657e5 / 0.201600e6 -0.631187e6 / 0.1360800e7 0 0.496403e6 / 0.907200e6 -0.308533e6 / 0.5443200e7; 0.631e3 / 0.28800e5 0.8467e4 / 0.453600e6 -0.400139e6 / 0.5443200e7 -0.496403e6 / 0.907200e6 0 0.1805647e7 / 0.2721600e7; -0.91e2 / 0.7776e4 0.33059e5 / 0.5443200e7 0.8789e4 / 0.302400e6 0.308533e6 / 0.5443200e7 -0.1805647e7 / 0.2721600e7 0;]; Q(1:6,1:6)=Q_U; - Q(m-5:m,m-5:m)=flipud( fliplr( -Q_U ) ); + Q(m-5:m,m-5:m)=rot90( -Q_U ,2 ); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Second derivative, 1st order accurate at first 6 boundary points - m2=1/12;m1=-16/12;m0=30/12; - M=m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% m2=1/12;m1=-16/12;m0=30/12; +% M=m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); %M=(1/12*diag(ones(m-2,1),2)-16/12*diag(ones(m-1,1),1)-16/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)+30/12*diag(ones(m,1),0)); + M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; M_U=[0.2386127e7 / 0.2177280e7 -0.515449e6 / 0.453600e6 -0.10781e5 / 0.777600e6 0.61567e5 / 0.1360800e7 0.6817e4 / 0.403200e6 -0.1069e4 / 0.136080e6; -0.515449e6 / 0.453600e6 0.4756039e7 / 0.2177280e7 -0.1270009e7 / 0.1360800e7 -0.3751e4 / 0.28800e5 0.3067e4 / 0.680400e6 0.119459e6 / 0.10886400e8; -0.10781e5 / 0.777600e6 -0.1270009e7 / 0.1360800e7 0.111623e6 / 0.60480e5 -0.555587e6 / 0.680400e6 -0.551339e6 / 0.5443200e7 0.8789e4 / 0.453600e6; 0.61567e5 / 0.1360800e7 -0.3751e4 / 0.28800e5 -0.555587e6 / 0.680400e6 0.1025327e7 / 0.544320e6 -0.464003e6 / 0.453600e6 0.222133e6 / 0.5443200e7; 0.6817e4 / 0.403200e6 0.3067e4 / 0.680400e6 -0.551339e6 / 0.5443200e7 -0.464003e6 / 0.453600e6 0.5074159e7 / 0.2177280e7 -0.1784047e7 / 0.1360800e7; -0.1069e4 / 0.136080e6 0.119459e6 / 0.10886400e8 0.8789e4 / 0.453600e6 0.222133e6 / 0.5443200e7 -0.1784047e7 / 0.1360800e7 0.1812749e7 / 0.725760e6;]; M(1:6,1:6)=M_U; - M(m-5:m,m-5:m)=flipud( fliplr( M_U ) ); + M(m-5:m,m-5:m)=rot90( M_U ,2 ); M=M/h; S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:4)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-3:m)=fliplr(-S_U); - D2=HI*(-M-e_1*S_1+e_m*S_m); + D2=H\(-M-e_1*S_1+e_m*S_m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -91,7 +94,10 @@ % Third derivative, 1st order accurate at first 6 boundary points q3=-1/8;q2=1;q1=-13/8; - Q3=q3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q3=q3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + diags = -3:3; + stencil = [-q3,-q2,-q1,0,q1,q2,q3]; + Q3 = stripeMatrix(stencil, diags, m); %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3)); @@ -99,26 +105,30 @@ Q3_U = [0 -0.88471e5 / 0.67200e5 0.58139e5 / 0.33600e5 -0.1167e4 / 0.2800e4 -0.89e2 / 0.11200e5 0.7e1 / 0.640e3; 0.88471e5 / 0.67200e5 0 -0.43723e5 / 0.16800e5 0.46783e5 / 0.33600e5 -0.191e3 / 0.3200e4 -0.1567e4 / 0.33600e5; -0.58139e5 / 0.33600e5 0.43723e5 / 0.16800e5 0 -0.4049e4 / 0.2400e4 0.29083e5 / 0.33600e5 -0.71e2 / 0.1400e4; 0.1167e4 / 0.2800e4 -0.46783e5 / 0.33600e5 0.4049e4 / 0.2400e4 0 -0.8591e4 / 0.5600e4 0.10613e5 / 0.11200e5; 0.89e2 / 0.11200e5 0.191e3 / 0.3200e4 -0.29083e5 / 0.33600e5 0.8591e4 / 0.5600e4 0 -0.108271e6 / 0.67200e5; -0.7e1 / 0.640e3 0.1567e4 / 0.33600e5 0.71e2 / 0.1400e4 -0.10613e5 / 0.11200e5 0.108271e6 / 0.67200e5 0;]; Q3(1:6,1:6)=Q3_U; - Q3(m-5:m,m-5:m)=flipud( fliplr( -Q3_U ) ); + Q3(m-5:m,m-5:m)=rot90( -Q3_U ,2 ); Q3=Q3/h^2; S2_U=[2 -5 4 -1;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:4)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-3:m)=fliplr(S2_U); - D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; + D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx m3=-1/6;m2=2;m1=-13/2;m0=28/3; - M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + diags = -3:3; + left_stencil = [m3,m2,m1]; + stencil = [left_stencil,m0,fliplr(left_stencil)]; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); @@ -126,16 +136,16 @@ M4(1:6,1:6)=M4_U; - M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); + M4(m-5:m,m-5:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-1 3 -3 1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:4)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-3:m)=fliplr(-S3_U); - D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); + D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); % L=h*(m-1);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_6.m --- a/+sbp/+implementations/d4_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -31,30 +31,33 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H_U=[0.318365e6 / 0.1016064e7 0 0 0 0 0 0 0; 0 0.145979e6 / 0.103680e6 0 0 0 0 0 0; 0 0 0.139177e6 / 0.241920e6 0 0 0 0 0; 0 0 0 0.964969e6 / 0.725760e6 0 0 0 0; 0 0 0 0 0.593477e6 / 0.725760e6 0 0 0; 0 0 0 0 0 0.52009e5 / 0.48384e5 0 0; 0 0 0 0 0 0 0.141893e6 / 0.145152e6 0; 0 0 0 0 0 0 0 0.1019713e7 / 0.1016064e7;]; H(1:8,1:8)=H_U; - H(m-7:m,m-7:m)=fliplr(flipud(H_U)); + H(m-7:m,m-7:m)=rot90(H_U,2); H=H*h; HI=inv(H); % First derivative SBP operator, 1st order accurate at first 6 boundary points - q3=1/60;q2=-3/20;q1=3/4; - Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% q3=1/60;q2=-3/20;q1=3/4; +% Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + diags = -3:3; + stencil = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; + Q = stripeMatrix(stencil, diags, m); Q_U = [0 0.1547358409e10 / 0.2421619200e10 -0.422423e6 / 0.11211200e8 -0.1002751721e10 / 0.8717829120e10 -0.15605263e8 / 0.484323840e9 0.1023865e7 / 0.24216192e8 0.291943739e9 / 0.21794572800e11 -0.24659e5 / 0.2534400e7; -0.1547358409e10 / 0.2421619200e10 0 0.23031829e8 / 0.62899200e8 0.10784027e8 / 0.34594560e8 0.2859215e7 / 0.31135104e8 -0.45982103e8 / 0.345945600e9 -0.26681e5 / 0.1182720e7 0.538846039e9 / 0.21794572800e11; 0.422423e6 / 0.11211200e8 -0.23031829e8 / 0.62899200e8 0 0.28368209e8 / 0.69189120e8 -0.9693137e7 / 0.69189120e8 0.1289363e7 / 0.17740800e8 -0.39181e5 / 0.5491200e7 -0.168647e6 / 0.24216192e8; 0.1002751721e10 / 0.8717829120e10 -0.10784027e8 / 0.34594560e8 -0.28368209e8 / 0.69189120e8 0 0.5833151e7 / 0.10644480e8 0.4353179e7 / 0.69189120e8 0.2462459e7 / 0.155675520e9 -0.215471e6 / 0.10762752e8; 0.15605263e8 / 0.484323840e9 -0.2859215e7 / 0.31135104e8 0.9693137e7 / 0.69189120e8 -0.5833151e7 / 0.10644480e8 0 0.7521509e7 / 0.13837824e8 -0.1013231e7 / 0.11531520e8 0.103152839e9 / 0.8717829120e10; -0.1023865e7 / 0.24216192e8 0.45982103e8 / 0.345945600e9 -0.1289363e7 / 0.17740800e8 -0.4353179e7 / 0.69189120e8 -0.7521509e7 / 0.13837824e8 0 0.67795697e8 / 0.98841600e8 -0.17263733e8 / 0.151351200e9; -0.291943739e9 / 0.21794572800e11 0.26681e5 / 0.1182720e7 0.39181e5 / 0.5491200e7 -0.2462459e7 / 0.155675520e9 0.1013231e7 / 0.11531520e8 -0.67795697e8 / 0.98841600e8 0 0.1769933569e10 / 0.2421619200e10; 0.24659e5 / 0.2534400e7 -0.538846039e9 / 0.21794572800e11 0.168647e6 / 0.24216192e8 0.215471e6 / 0.10762752e8 -0.103152839e9 / 0.8717829120e10 0.17263733e8 / 0.151351200e9 -0.1769933569e10 / 0.2421619200e10 0;]; Q(1:8,1:8)=Q_U; - Q(m-7:m,m-7:m)=flipud( fliplr( -Q_U ) ); + Q(m-7:m,m-7:m)=rot90( -Q_U ,2 ); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -63,22 +66,26 @@ % Second derivative, 1st order accurate at first 6 boundary points m3=-1/90;m2=3/20;m1=-3/2;m0=49/18; - M=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + diags = -3:3; + stencil = [m3,m2,m1,m0,m1,m2,m3]; + M = stripeMatrix(stencil, diags, m); + M_U=[0.4347276223e10 / 0.3736212480e10 -0.1534657609e10 / 0.1210809600e10 0.68879e5 / 0.3057600e7 0.1092927401e10 / 0.13076743680e11 0.18145423e8 / 0.968647680e9 -0.1143817e7 / 0.60540480e8 -0.355447739e9 / 0.65383718400e11 0.56081e5 / 0.16473600e8; -0.1534657609e10 / 0.1210809600e10 0.42416226217e11 / 0.18681062400e11 -0.228654119e9 / 0.345945600e9 -0.12245627e8 / 0.34594560e8 -0.2995295e7 / 0.46702656e8 0.52836503e8 / 0.691891200e9 0.119351e6 / 0.12812800e8 -0.634102039e9 / 0.65383718400e11; 0.68879e5 / 0.3057600e7 -0.228654119e9 / 0.345945600e9 0.5399287e7 / 0.4193280e7 -0.24739409e8 / 0.34594560e8 0.7878737e7 / 0.69189120e8 -0.1917829e7 / 0.31449600e8 0.39727e5 / 0.3660800e7 0.10259e5 / 0.4656960e7; 0.1092927401e10 / 0.13076743680e11 -0.12245627e8 / 0.34594560e8 -0.24739409e8 / 0.34594560e8 0.7780367599e10 / 0.3736212480e10 -0.70085363e8 / 0.69189120e8 -0.500209e6 / 0.6289920e7 -0.311543e6 / 0.17962560e8 0.278191e6 / 0.21525504e8; 0.18145423e8 / 0.968647680e9 -0.2995295e7 / 0.46702656e8 0.7878737e7 / 0.69189120e8 -0.70085363e8 / 0.69189120e8 0.7116321131e10 / 0.3736212480e10 -0.545081e6 / 0.532224e6 0.811631e6 / 0.11531520e8 -0.84101639e8 / 0.13076743680e11; -0.1143817e7 / 0.60540480e8 0.52836503e8 / 0.691891200e9 -0.1917829e7 / 0.31449600e8 -0.500209e6 / 0.6289920e7 -0.545081e6 / 0.532224e6 0.324760747e9 / 0.138378240e9 -0.65995697e8 / 0.49420800e8 0.1469203e7 / 0.13759200e8; -0.355447739e9 / 0.65383718400e11 0.119351e6 / 0.12812800e8 0.39727e5 / 0.3660800e7 -0.311543e6 / 0.17962560e8 0.811631e6 / 0.11531520e8 -0.65995697e8 / 0.49420800e8 0.48284442317e11 / 0.18681062400e11 -0.1762877569e10 / 0.1210809600e10; 0.56081e5 / 0.16473600e8 -0.634102039e9 / 0.65383718400e11 0.10259e5 / 0.4656960e7 0.278191e6 / 0.21525504e8 -0.84101639e8 / 0.13076743680e11 0.1469203e7 / 0.13759200e8 -0.1762877569e10 / 0.1210809600e10 0.10117212851e11 / 0.3736212480e10;]; M(1:8,1:8)=M_U; - M(m-7:m,m-7:m)=flipud( fliplr( M_U ) ); + M(m-7:m,m-7:m)=rot90( M_U ,2 ); M=M/h; S_U=[-0.25e2 / 0.12e2 4 -3 0.4e1 / 0.3e1 -0.1e1 / 0.4e1;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:5)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-4:m)=fliplr(-S_U); - D2=HI*(-M-e_1*S_1+e_m*S_m); + D2=H\(-M-e_1*S_1+e_m*S_m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -88,54 +95,59 @@ % Third derivative, 1st order accurate at first 6 boundary points q4=7/240;q3=-3/10;q2=169/120;q1=-61/30; - Q3=q4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+q3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); - +% Q3=q4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+q3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + diags = -4:4; + stencil = [-q4,-q3,-q2,-q1,0,q1,q2,q3,q4]; + Q3 = stripeMatrix(stencil, diags, m); %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3)); Q3_U = [0 -0.10882810591e11 / 0.5811886080e10 0.398713069e9 / 0.132088320e9 -0.1746657571e10 / 0.1162377216e10 0.56050639e8 / 0.145297152e9 -0.11473393e8 / 0.1162377216e10 -0.38062741e8 / 0.1452971520e10 0.30473e5 / 0.4392960e7; 0.10882810591e11 / 0.5811886080e10 0 -0.3720544343e10 / 0.830269440e9 0.767707019e9 / 0.207567360e9 -0.1047978301e10 / 0.830269440e9 0.1240729e7 / 0.14826240e8 0.6807397e7 / 0.55351296e8 -0.50022767e8 / 0.1452971520e10; -0.398713069e9 / 0.132088320e9 0.3720544343e10 / 0.830269440e9 0 -0.2870078009e10 / 0.830269440e9 0.74962049e8 / 0.29652480e8 -0.12944857e8 / 0.30750720e8 -0.17846623e8 / 0.103783680e9 0.68707591e8 / 0.1162377216e10; 0.1746657571e10 / 0.1162377216e10 -0.767707019e9 / 0.207567360e9 0.2870078009e10 / 0.830269440e9 0 -0.727867087e9 / 0.276756480e9 0.327603877e9 / 0.207567360e9 -0.175223717e9 / 0.830269440e9 0.1353613e7 / 0.726485760e9; -0.56050639e8 / 0.145297152e9 0.1047978301e10 / 0.830269440e9 -0.74962049e8 / 0.29652480e8 0.727867087e9 / 0.276756480e9 0 -0.1804641793e10 / 0.830269440e9 0.311038417e9 / 0.207567360e9 -0.1932566239e10 / 0.5811886080e10; 0.11473393e8 / 0.1162377216e10 -0.1240729e7 / 0.14826240e8 0.12944857e8 / 0.30750720e8 -0.327603877e9 / 0.207567360e9 0.1804641793e10 / 0.830269440e9 0 -0.1760949511e10 / 0.830269440e9 0.2105883973e10 / 0.1452971520e10; 0.38062741e8 / 0.1452971520e10 -0.6807397e7 / 0.55351296e8 0.17846623e8 / 0.103783680e9 0.175223717e9 / 0.830269440e9 -0.311038417e9 / 0.207567360e9 0.1760949511e10 / 0.830269440e9 0 -0.1081094773e10 / 0.528353280e9; -0.30473e5 / 0.4392960e7 0.50022767e8 / 0.1452971520e10 -0.68707591e8 / 0.1162377216e10 -0.1353613e7 / 0.726485760e9 0.1932566239e10 / 0.5811886080e10 -0.2105883973e10 / 0.1452971520e10 0.1081094773e10 / 0.528353280e9 0;]; Q3(1:8,1:8)=Q3_U; - Q3(m-7:m,m-7:m)=flipud( fliplr( -Q3_U ) ); + Q3(m-7:m,m-7:m)=rot90( -Q3_U ,2 ); Q3=Q3/h^2; S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:5)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-4:m)=fliplr(S2_U); - D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; + D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8; - M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); - +% M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + diags = -4:4; + left_stencil = [m4,m3,m2,m1]; + stencil = [left_stencil,m0,fliplr(left_stencil)]; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); M4_U=[0.40833734273e11 / 0.10761070320e11 -0.162181998421e12 / 0.16397821440e11 0.4696168417e10 / 0.521748864e9 -0.245714671483e12 / 0.68870850048e11 0.2185939219e10 / 0.2869618752e10 -0.15248255797e11 / 0.114784750080e12 0.345156907e9 / 0.12298366080e11 0.6388381e7 / 0.1093188096e10; -0.162181998421e12 / 0.16397821440e11 0.147281127041e12 / 0.5380535160e10 -0.3072614435609e13 / 0.114784750080e12 0.320122985851e12 / 0.28696187520e11 -0.768046031383e12 / 0.344354250240e12 0.7861605187e10 / 0.14348093760e11 -0.803762437e9 / 0.4251287040e10 0.167394281e9 / 0.86088562560e11; 0.4696168417e10 / 0.521748864e9 -0.3072614435609e13 / 0.114784750080e12 0.139712483333e12 / 0.4782697920e10 -0.1634124842747e13 / 0.114784750080e12 0.90855193447e11 / 0.28696187520e11 -0.26412188989e11 / 0.38261583360e11 0.668741173e9 / 0.1793511720e10 -0.132673781e9 / 0.2342545920e10; -0.245714671483e12 / 0.68870850048e11 0.320122985851e12 / 0.28696187520e11 -0.1634124842747e13 / 0.114784750080e12 0.437353997177e12 / 0.43044281280e11 -0.172873969321e12 / 0.38261583360e11 0.34759553483e11 / 0.28696187520e11 -0.98928859751e11 / 0.344354250240e12 0.295000207e9 / 0.3587023440e10; 0.2185939219e10 / 0.2869618752e10 -0.768046031383e12 / 0.344354250240e12 0.90855193447e11 / 0.28696187520e11 -0.172873969321e12 / 0.38261583360e11 0.126711914423e12 / 0.21522140640e11 -0.520477408939e12 / 0.114784750080e12 0.49581230003e11 / 0.28696187520e11 -0.99640101991e11 / 0.344354250240e12; -0.15248255797e11 / 0.114784750080e12 0.7861605187e10 / 0.14348093760e11 -0.26412188989e11 / 0.38261583360e11 0.34759553483e11 / 0.28696187520e11 -0.520477408939e12 / 0.114784750080e12 0.19422074929e11 / 0.2391348960e10 -0.772894368601e12 / 0.114784750080e12 0.10579712849e11 / 0.4099455360e10; 0.345156907e9 / 0.12298366080e11 -0.803762437e9 / 0.4251287040e10 0.668741173e9 / 0.1793511720e10 -0.98928859751e11 / 0.344354250240e12 0.49581230003e11 / 0.28696187520e11 -0.772894368601e12 / 0.114784750080e12 0.456715296239e12 / 0.43044281280e11 -0.915425403107e12 / 0.114784750080e12; 0.6388381e7 / 0.1093188096e10 0.167394281e9 / 0.86088562560e11 -0.132673781e9 / 0.2342545920e10 0.295000207e9 / 0.3587023440e10 -0.99640101991e11 / 0.344354250240e12 0.10579712849e11 / 0.4099455360e10 -0.915425403107e12 / 0.114784750080e12 0.488029542379e12 / 0.43044281280e11;]; M4(1:8,1:8)=M4_U; - M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) ); + M4(m-7:m,m-7:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:5)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-4:m)=fliplr(-S3_U); - D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); + D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); - L=h*(m-1); +% L=h*(m-1); % x1=linspace(0,L,m)'; % x2=x1.^2/fac(2);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_2.m --- a/+sbp/+implementations/d4_compatible_2.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_2.m Fri Sep 09 14:53:41 2016 +0200 @@ -26,12 +26,12 @@ % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm, % vilket ?r n?dv?ndigt f?r stabilitet - BP = 1; + BP = 4; if(m<2*BP) error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0);H(1,1)=1/2;H(m,m)=1/2; + H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2; H=H*h; @@ -41,34 +41,42 @@ % First derivative SBP operator, 1st order accurate at first 6 boundary points q1=1/2; - Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + stencil = [-q1,0,q1]; + d = (length(stencil)-1)/2; + diags = -d:d; + Q = stripeMatrix(stencil, diags, m); %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Second derivative, 1st order accurate at first 6 boundary points - m1=-1;m0=2; - M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1; - M=M/h; +% m1=-1;m0=2; +% % M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1; +% stencil = [m2,m1,m0,m1,m2]; +% d = (length(stencil)-1)/2; +% diags = -d:d; +% M = stripeMatrix(stencil, diags, m); +% M=M/h; S_U=[-1 1]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:2)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-1:m)=fliplr(-S_U); - D2=HI*(-M-e_1*S_1+e_m*S_m); +% D2=H\(-M-e_1*S_1+e_m*S_m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -77,28 +85,32 @@ % Third derivative, 1st order accurate at first 6 boundary points - q2=1/2;q1=-1; - Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% q2=1/2;q1=-1; +% % Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% stencil = [-q2,-q1,0,q1,q2]; +% d = (length(stencil)-1)/2; +% diags = -d:d; +% Q3 = stripeMatrix(stencil, diags, m); %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3)); - Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;]; - Q3(1:4,1:4)=Q3_U; - Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) ); - Q3=Q3/h^2; +% Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;]; +% Q3(1:4,1:4)=Q3_U; +% Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 ); +% Q3=Q3/h^2; S2_U=[1 -2 1;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:3)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-2:m)=fliplr(S2_U); - D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; +% D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx @@ -112,16 +124,16 @@ M4(1:4,1:4)=M4_U; - M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) ); + M4(m-3:m,m-3:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-1 3 -3 1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:4)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-3:m)=fliplr(-S3_U); - D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); + D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_4.m --- a/+sbp/+implementations/d4_compatible_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -31,10 +31,10 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H_U=[0.3e1 / 0.11e2 0 0 0 0 0; 0 0.2125516311e10 / 0.1311004640e10 0 0 0 0; 0 0 0.278735189e9 / 0.1966506960e10 0 0 0; 0 0 0 0.285925927e9 / 0.163875580e9 0 0; 0 0 0 0 0.1284335339e10 / 0.1966506960e10 0; 0 0 0 0 0 0.4194024163e10 / 0.3933013920e10;]; H(1:6,1:6)=H_U; - H(m-5:m,m-5:m)=fliplr(flipud(H_U)); + H(m-5:m,m-5:m)=rot90(H_U,2); H=H*h; HI=inv(H); @@ -42,19 +42,23 @@ % First derivative SBP operator, 1st order accurate at first 6 boundary points q2=-1/12;q1=8/12; - Q=q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q=q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + stencil = [-q2,-q1,0,q1,q2]; + d = (length(stencil)-1)/2; + diags = -d:d; + Q = stripeMatrix(stencil, diags, m); %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); Q_U = [0 0.9e1 / 0.11e2 -0.9e1 / 0.22e2 0.1e1 / 0.11e2 0 0; -0.9e1 / 0.11e2 0 0.2595224893e10 / 0.2622009280e10 -0.151435707e9 / 0.327751160e9 0.1112665611e10 / 0.2622009280e10 -0.1290899e7 / 0.9639740e7; 0.9e1 / 0.22e2 -0.2595224893e10 / 0.2622009280e10 0 0.1468436423e10 / 0.983253480e9 -0.1194603401e10 / 0.983253480e9 0.72033031e8 / 0.238364480e9; -0.1e1 / 0.11e2 0.151435707e9 / 0.327751160e9 -0.1468436423e10 / 0.983253480e9 0 0.439819541e9 / 0.327751160e9 -0.215942641e9 / 0.983253480e9; 0 -0.1112665611e10 / 0.2622009280e10 0.1194603401e10 / 0.983253480e9 -0.439819541e9 / 0.327751160e9 0 0.1664113643e10 / 0.2622009280e10; 0 0.1290899e7 / 0.9639740e7 -0.72033031e8 / 0.238364480e9 0.215942641e9 / 0.983253480e9 -0.1664113643e10 / 0.2622009280e10 0;]; Q(1:6,1:6)=Q_U; - Q(m-5:m,m-5:m)=flipud( fliplr( -Q_U ) ); + Q(m-5:m,m-5:m)=rot90( -Q_U ,2 ); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -72,13 +76,13 @@ % M=M/h; % S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:4)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-3:m)=fliplr(-S_U); - % D2=HI*(-M-e_1*S_1+e_m*S_m); + % D2=H\(-M-e_1*S_1+e_m*S_m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -102,20 +106,25 @@ S2_U=[2 -5 4 -1;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:4)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-3:m)=fliplr(S2_U); - %D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; + %D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx m3=-1/6;m2=2;m1=-13/2;m0=28/3; - M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + + stencil = [m3,m2,m1,m0,m1,m2,m3]; + d = (length(stencil)-1)/2; + diags = -d:d; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); @@ -123,16 +132,16 @@ M4(1:6,1:6)=M4_U; - M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); + M4(m-5:m,m-5:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-1 3 -3 1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:4)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-3:m)=fliplr(-S3_U); - D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); + D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); S_1 = S_1'; S_m = S_m';
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_6.m --- a/+sbp/+implementations/d4_compatible_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -31,11 +31,11 @@ error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H_U=[0.7493827e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5534051e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.104561e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.260503e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.43237e5 / 0.103680e6 0 0 0; 0 0 0 0 0 0.514081e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3356179e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25631027e8 / 0.25401600e8;]; H(1:8,1:8)=H_U; - H(m-7:m,m-7:m)=fliplr(flipud(H_U)); + H(m-7:m,m-7:m)=rot90(H_U,2); H=H*h; HI=inv(H); @@ -43,18 +43,22 @@ % First derivative SBP operator, 3rd order accurate at first 8 boundary points q3=1/60;q2=-3/20;q1=3/4; - Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + stencil = [-q3,-q2,-q1,0,q1,q2,q3]; + d = (length(stencil)-1)/2; + diags = -d:d; + Q = stripeMatrix(stencil, diags, m); Q_U = [0 0.26431903e8 / 0.43545600e8 0.39791489e8 / 0.152409600e9 -0.12747751e8 / 0.16934400e8 0.76099447e8 / 0.152409600e9 -0.1397443e7 / 0.12192768e8 0 0; -0.26431903e8 / 0.43545600e8 0 -0.13847476213e11 / 0.19559232000e11 0.35844843977e11 / 0.11735539200e11 -0.63413503537e11 / 0.23471078400e11 0.4764412871e10 / 0.3911846400e10 -0.1668252557e10 / 0.5867769600e10 0.842644697e9 / 0.29338848000e11; -0.39791489e8 / 0.152409600e9 0.13847476213e11 / 0.19559232000e11 0 -0.73834802771e11 / 0.23471078400e11 0.1802732209e10 / 0.325987200e9 -0.65514173e8 / 0.16299360e8 0.79341409141e11 / 0.58677696000e11 -0.1282384321e10 / 0.7823692800e10; 0.12747751e8 / 0.16934400e8 -0.35844843977e11 / 0.11735539200e11 0.73834802771e11 / 0.23471078400e11 0 -0.5274106087e10 / 0.1173553920e10 0.33743985841e11 / 0.5867769600e10 -0.6482602549e10 / 0.2607897600e10 0.1506017269e10 / 0.3911846400e10; -0.76099447e8 / 0.152409600e9 0.63413503537e11 / 0.23471078400e11 -0.1802732209e10 / 0.325987200e9 0.5274106087e10 / 0.1173553920e10 0 -0.7165829063e10 / 0.2607897600e10 0.23903110999e11 / 0.11735539200e11 -0.5346675911e10 / 0.11735539200e11; 0.1397443e7 / 0.12192768e8 -0.4764412871e10 / 0.3911846400e10 0.65514173e8 / 0.16299360e8 -0.33743985841e11 / 0.5867769600e10 0.7165829063e10 / 0.2607897600e10 0 -0.1060918223e10 / 0.11735539200e11 0.628353989e9 / 0.3911846400e10; 0 0.1668252557e10 / 0.5867769600e10 -0.79341409141e11 / 0.58677696000e11 0.6482602549e10 / 0.2607897600e10 -0.23903110999e11 / 0.11735539200e11 0.1060918223e10 / 0.11735539200e11 0 0.25889988599e11 / 0.39118464000e11; 0 -0.842644697e9 / 0.29338848000e11 0.1282384321e10 / 0.7823692800e10 -0.1506017269e10 / 0.3911846400e10 0.5346675911e10 / 0.11735539200e11 -0.628353989e9 / 0.3911846400e10 -0.25889988599e11 / 0.39118464000e11 0;]; Q(1:8,1:8)=Q_U; - Q(m-7:m,m-7:m)=flipud( fliplr( -Q_U ) ); + Q(m-7:m,m-7:m)=rot90( -Q_U ,2 ); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -72,13 +76,13 @@ % M=M/h; S_U=[-0.12700800e8 / 0.7493827e7 0.185023321e9 / 0.89925924e8 0.39791489e8 / 0.44962962e8 -0.38243253e8 / 0.14987654e8 0.76099447e8 / 0.44962962e8 -0.34936075e8 / 0.89925924e8;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:6)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-5:m)=fliplr(-S_U); - %D2=HI*(-M-e_1*S_1+e_m*S_m); + %D2=H\(-M-e_1*S_1+e_m*S_m); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -102,37 +106,40 @@ % % S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:5)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-4:m)=fliplr(S2_U); % % % - % D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; + % D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8; - M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); - +% M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + stencil = [m4,m3,m2,m1,m0,m1,m2,m3,m4]; + d = (length(stencil)-1)/2; + diags = -d:d; + M4 = stripeMatrix(stencil, diags, m); M4_U=[0.600485868980522851e18 / 0.274825314114120000e18 -0.1421010223681841e16 / 0.348984525859200e15 0.38908412970187e14 / 0.1586293299360000e16 0.10224077451922837e17 / 0.2243471951952000e16 -0.7577302712815639e16 / 0.1744922629296000e16 0.138091642084013e15 / 0.59351109840000e14 -0.3775041725375197e16 / 0.4486943903904000e16 0.9907210230881393e16 / 0.61072292025360000e17; -0.1421010223681841e16 / 0.348984525859200e15 0.3985852497808703e16 / 0.407903991264000e15 -0.90048788923861e14 / 0.15579666333000e14 -0.4312795866499e13 / 0.997098645312e12 0.4414634708891947e16 / 0.448694390390400e15 -0.886174803100459e15 / 0.99709864531200e14 0.4333e4 / 0.1000e4 -0.13800578064893047e17 / 0.15704303663664000e17; 0.38908412970187e14 / 0.1586293299360000e16 -0.90048788923861e14 / 0.15579666333000e14 0.2071682582321887e16 / 0.113306664240000e15 -0.769471337294003e15 / 0.41545776888000e14 0.112191585452033e15 / 0.166183107552000e15 0.7204491902193671e16 / 0.623186653320000e15 -0.24847093554379e14 / 0.3115933266600e13 0.943854037768721e15 / 0.545288321655000e15; 0.10224077451922837e17 / 0.2243471951952000e16 -0.4312795866499e13 / 0.997098645312e12 -0.769471337294003e15 / 0.41545776888000e14 0.3086874339649421e16 / 0.81580798252800e14 -0.396009005312111e15 / 0.16618310755200e14 0.348854811893087e15 / 0.249274661328000e15 0.895954627955053e15 / 0.224347195195200e15 -0.184881685054543e15 / 0.166183107552000e15; -0.7577302712815639e16 / 0.1744922629296000e16 0.4414634708891947e16 / 0.448694390390400e15 0.112191585452033e15 / 0.166183107552000e15 -0.396009005312111e15 / 0.16618310755200e14 0.3774861828677557e16 / 0.112173597597600e15 -0.5693689108983593e16 / 0.249274661328000e15 0.803944126167107e15 / 0.99709864531200e14 -0.19547569411550791e17 / 0.15704303663664000e17; 0.138091642084013e15 / 0.59351109840000e14 -0.886174803100459e15 / 0.99709864531200e14 0.7204491902193671e16 / 0.623186653320000e15 0.348854811893087e15 / 0.249274661328000e15 -0.5693689108983593e16 / 0.249274661328000e15 0.73965842628398389e17 / 0.2492746613280000e16 -0.2184472662036043e16 / 0.124637330664000e15 0.46667e5 / 0.10000e5; -0.3775041725375197e16 / 0.4486943903904000e16 0.4333e4 / 0.1000e4 -0.24847093554379e14 / 0.3115933266600e13 0.895954627955053e15 / 0.224347195195200e15 0.803944126167107e15 / 0.99709864531200e14 -0.2184472662036043e16 / 0.124637330664000e15 0.37593640125444199e17 / 0.2243471951952000e16 -0.37e2 / 0.4e1; 0.9907210230881393e16 / 0.61072292025360000e17 -0.13800578064893047e17 / 0.15704303663664000e17 0.943854037768721e15 / 0.545288321655000e15 -0.184881685054543e15 / 0.166183107552000e15 -0.19547569411550791e17 / 0.15704303663664000e17 0.46667e5 / 0.10000e5 -0.37e2 / 0.4e1 0.12766926490502478779e20 / 0.1099301256456480000e19;]; M4(1:8,1:8)=M4_U; - M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) ); + M4(m-7:m,m-7:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:5)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-4:m)=fliplr(-S3_U); - D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); + D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); S_1 = S_1'; S_m = S_m';
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_halfvariable_2.m --- a/+sbp/+implementations/d4_compatible_halfvariable_2.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_halfvariable_2.m Fri Sep 09 14:53:41 2016 +0200 @@ -27,12 +27,12 @@ % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm, % vilket ?r n?dv?ndigt f?r stabilitet - BP = 1; + BP = 4; if(m<2*BP) error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0);H(1,1)=1/2;H(m,m)=1/2; + H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2; H=H*h; @@ -42,16 +42,20 @@ % First derivative SBP operator, 1st order accurate at first 6 boundary points q1=1/2; - Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + stencil = [-q1,0,q1]; + d = (length(stencil)-1)/2; + diags = -d:d; + Q = stripeMatrix(stencil, diags, m); %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; + D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -69,9 +73,9 @@ % Require a vector c with the koeffients S_U=[-3/2 2 -1/2]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:3)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-2:m)=fliplr(-S_U); S_1 = S_1'; @@ -115,35 +119,43 @@ % Third derivative, 1st order accurate at first 6 boundary points q2=1/2;q1=-1; - Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); +% Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); + stencil = [-q2,-q1,0,q1,q2]; + d = (length(stencil)-1)/2; + diags = -d:d; + Q3 = stripeMatrix(stencil, diags, m); %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3)); Q3_U = [0 -0.13e2 / 0.16e2 0.7e1 / 0.8e1 -0.1e1 / 0.16e2; 0.13e2 / 0.16e2 0 -0.23e2 / 0.16e2 0.5e1 / 0.8e1; -0.7e1 / 0.8e1 0.23e2 / 0.16e2 0 -0.17e2 / 0.16e2; 0.1e1 / 0.16e2 -0.5e1 / 0.8e1 0.17e2 / 0.16e2 0;]; Q3(1:4,1:4)=Q3_U; - Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) ); + Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 ); Q3=Q3/h^2; S2_U=[1 -2 1;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:3)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-2:m)=fliplr(S2_U); S2_1 = S2_1'; S2_m = S2_m'; - D3=HI*(Q3 - e_1*S2_1' + e_m*S2_m' +1/2*S_1*S_1' -1/2*S_m*S_m' ) ; + D3=HI*(Q3 - e_1*S2_1' + e_m*S2_m' +1/2*(S_1*S_1') -1/2*(S_m*S_m') ) ; % Fourth derivative, 0th order accurate at first 6 boundary points (still % yield 4th order convergence if stable: for example u_tt=-u_xxxx m2=1;m1=-4;m0=6; - M4=m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M4=m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + stencil = [m2,m1,m0,m1,m2]; + d = (length(stencil)-1)/2; + diags = -d:d; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); @@ -152,13 +164,13 @@ M4(1:4,1:4)=M4_U; - M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) ); + M4(m-3:m,m-3:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-1 3 -3 1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:4)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-3:m)=fliplr(-S3_U); S3_1 = S3_1'; S3_m = S3_m';
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_halfvariable_4.m --- a/+sbp/+implementations/d4_compatible_halfvariable_4.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_halfvariable_4.m Fri Sep 09 14:53:41 2016 +0200 @@ -13,45 +13,48 @@ %h=1/(m-1); %h=1; - BP = 4; + BP = 6; if(m<2*BP) error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - c=ones(m,1); +% c=ones(m,1); - H=diag(ones(m,1),0); + H=speye(m,m); H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]); - H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48]))); + H(m-3:m,m-3:m)=rot90(diag([17/48 59/48 43/48 49/48]),2); H=H*h; HI=inv(H); HI = sparse(HI); - Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); - Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;]; - Q(1:4,1:4)=Q_U; - Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) ); +% Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); + e=ones(m,1); +% Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; +% Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;]; +% Q(1:4,1:4)=Q_U; +% Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); - 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; - D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; +% D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;]; - M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); +% M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); + M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; M(1:4,1:4)=M_U; - M(m-3:m,m-3:m)=flipud( fliplr( M_U ) ); + M(m-3:m,m-3:m)=rot90( M_U ,2 ); M=M/h; S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:4)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-3:m)=fliplr(-S_U); S_1 = S_1'; S_m = S_m'; @@ -124,15 +127,19 @@ S2_U=[2 -5 4 -1;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:4)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-3:m)=fliplr(S2_U); S2_1 = S2_1'; S2_m = S2_m'; m3=-1/6;m2=2;m1=-13/2;m0=28/3; - M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M4=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + stencil = [m3,m2,m1,m0,m1,m2,m3]; + d = (length(stencil)-1)/2; + diags = -d:d; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); @@ -140,13 +147,13 @@ M4(1:6,1:6)=M4_U; - M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); + M4(m-5:m,m-5:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-1 3 -3 1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:4)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-3:m)=fliplr(-S3_U); S3_1 = S3_1'; S3_m = S3_m';
diff -r bfa130b7abf6 -r f7ac3cd6eeaa +sbp/+implementations/d4_compatible_halfvariable_6.m --- a/+sbp/+implementations/d4_compatible_halfvariable_6.m Fri Sep 09 11:03:13 2016 +0200 +++ b/+sbp/+implementations/d4_compatible_halfvariable_6.m Fri Sep 09 14:53:41 2016 +0200 @@ -26,61 +26,61 @@ % x=1:h:m*h;x=x'; % c=x.^0; - BP = 6; + BP = 8; if(m<2*BP) error(['Operator requires at least ' num2str(2*BP) ' grid points']); end - H=diag(ones(m,1),0); + H=speye(m,m); H(1:6,1:6)=diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, ... 43801/43200]); - H(m-5:m,m-5:m)=fliplr(flipud(diag([13649/43200,12013/8640, ... - 2711/4320,5359/4320,7877/8640,43801/43200]))); + H(m-5:m,m-5:m)=rot90(diag([13649/43200,12013/8640, ... + 2711/4320,5359/4320,7877/8640,43801/43200]),2); - x1=0.70127127127127; - - - D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); - +% x1=0.70127127127127; - D1(1:6,1:9)=[-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ... - 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ... - x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; ... - -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, ... - -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, ... - -34560/12013*x1+240569/120130, 0, 0, 0; ... - 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, ... - 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, ... - 25920/2711*x1-182261/27110, 0, 0, 0; ... - -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, ... - -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, ... - -17280/5359*x1+710473/321540, 72/5359, 0, 0; ... - 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, ... - 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, ... - 0, 8640/7877*x1, -1296/7877, 144/7877, 0; ... - -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602,... - -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ... - -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801]; - D1(m-5:m,m-8:m)=flipud( fliplr(-D1(1:6,1:9))); - D1=D1/h; +% D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); +% +% +% +% D1(1:6,1:9)=[-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ... +% 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ... +% x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; ... +% -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, ... +% -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, ... +% -34560/12013*x1+240569/120130, 0, 0, 0; ... +% 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, ... +% 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, ... +% 25920/2711*x1-182261/27110, 0, 0, 0; ... +% -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, ... +% -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, ... +% -17280/5359*x1+710473/321540, 72/5359, 0, 0; ... +% 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, ... +% 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, ... +% 0, 8640/7877*x1, -1296/7877, 144/7877, 0; ... +% -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602,... +% -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ... +% -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801]; +% D1(m-5:m,m-8:m)=rot90( -D1(1:6,1:9),2); +% D1=D1/h; - 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; S_U=[-25/12, 4, -3, 4/3, -1/4]/h; - S_1=zeros(1,m); + S_1=sparse(1,m); S_1(1:5)=S_U; - S_m=zeros(1,m); + S_m=sparse(1,m); S_m(m-4:m)=fliplr(-S_U); S_1 = S_1'; S_m = S_m'; - %DS=zeros(m,m); + %DS=sparse(m,m); %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4]; %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]); %DS=diag(c)*DS/h; @@ -127,9 +127,9 @@ D2 = @D2_fun; S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2; - S2_1=zeros(1,m); + S2_1=sparse(1,m); S2_1(1:5)=S2_U; - S2_m=zeros(1,m); + S2_m=sparse(1,m); S2_m(m-4:m)=fliplr(S2_U); S2_1 = S2_1'; S2_m = S2_m'; @@ -142,7 +142,11 @@ % yield 5th order convergence if stable: for example u_tt=-u_xxxx m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8; - M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); +% M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0); + stencil = [m4,m3,m2,m1,m0,m1,m2,m3,m4]; + d = (length(stencil)-1)/2; + diags = -d:d; + M4 = stripeMatrix(stencil, diags, m); %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0)); @@ -150,13 +154,13 @@ M4(1:8,1:8)=M4_U; - M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) ); + M4(m-7:m,m-7:m)=rot90( M4_U ,2 ); M4=M4/h^3; S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3; - S3_1=zeros(1,m); + S3_1=sparse(1,m); S3_1(1:5)=S3_U; - S3_m=zeros(1,m); + S3_m=sparse(1,m); S3_m(m-4:m)=fliplr(-S3_U); S3_1 = S3_1'; S3_m = S3_m';