Mercurial > repos > public > sbplib
diff +sbp/+implementations/d2_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/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);