Mercurial > repos > public > sbplib
diff +sbp/+implementations/d4_compatible_halfvariable_6.m @ 267:f7ac3cd6eeaa operator_remake
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Fri, 09 Sep 2016 14:53:41 +0200 |
| parents | bfa130b7abf6 |
| children |
line wrap: on
line diff
--- 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';