Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_halfvariable_4.m @ 267:f7ac3cd6eeaa operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
11 11
12 %m=20; %problemstorlek 12 %m=20; %problemstorlek
13 %h=1/(m-1); 13 %h=1/(m-1);
14 %h=1; 14 %h=1;
15 15
16 BP = 4; 16 BP = 6;
17 if(m<2*BP) 17 if(m<2*BP)
18 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 18 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
19 end 19 end
20 20
21 c=ones(m,1); 21 % c=ones(m,1);
22 22
23 23
24 H=diag(ones(m,1),0); 24 H=speye(m,m);
25 H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]); 25 H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]);
26 H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48]))); 26 H(m-3:m,m-3:m)=rot90(diag([17/48 59/48 43/48 49/48]),2);
27 H=H*h; 27 H=H*h;
28 HI=inv(H); 28 HI=inv(H);
29 HI = sparse(HI); 29 HI = sparse(HI);
30 30
31 31
32 32
33 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)); 33 % 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));
34 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;]; 34 e=ones(m,1);
35 Q(1:4,1:4)=Q_U; 35 % Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
36 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) ); 36 % 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;];
37 % Q(1:4,1:4)=Q_U;
38 % Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
37 39
38 e_1=zeros(m,1);e_1(1)=1; 40 e_1=sparse(m,1);e_1(1)=1;
39 e_m=zeros(m,1);e_m(m)=1; 41 e_m=sparse(m,1);e_m(m)=1;
40 42
41 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 43 % D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
42 44
43 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;]; 45 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;];
44 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)); 46 % 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));
47 M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12;
45 48
46 M(1:4,1:4)=M_U; 49 M(1:4,1:4)=M_U;
47 50
48 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) ); 51 M(m-3:m,m-3:m)=rot90( M_U ,2 );
49 M=M/h; 52 M=M/h;
50 53
51 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h; 54 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
52 S_1=zeros(1,m); 55 S_1=sparse(1,m);
53 S_1(1:4)=S_U; 56 S_1(1:4)=S_U;
54 S_m=zeros(1,m); 57 S_m=sparse(1,m);
55 S_m(m-3:m)=fliplr(-S_U); 58 S_m(m-3:m)=fliplr(-S_U);
56 S_1 = S_1'; 59 S_1 = S_1';
57 S_m = S_m'; 60 S_m = S_m';
58 61
59 62
122 end 125 end
123 D2 = @D2_fun; 126 D2 = @D2_fun;
124 127
125 128
126 S2_U=[2 -5 4 -1;]/h^2; 129 S2_U=[2 -5 4 -1;]/h^2;
127 S2_1=zeros(1,m); 130 S2_1=sparse(1,m);
128 S2_1(1:4)=S2_U; 131 S2_1(1:4)=S2_U;
129 S2_m=zeros(1,m); 132 S2_m=sparse(1,m);
130 S2_m(m-3:m)=fliplr(S2_U); 133 S2_m(m-3:m)=fliplr(S2_U);
131 S2_1 = S2_1'; 134 S2_1 = S2_1';
132 S2_m = S2_m'; 135 S2_m = S2_m';
133 136
134 m3=-1/6;m2=2;m1=-13/2;m0=28/3; 137 m3=-1/6;m2=2;m1=-13/2;m0=28/3;
135 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); 138 % 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);
139 stencil = [m3,m2,m1,m0,m1,m2,m3];
140 d = (length(stencil)-1)/2;
141 diags = -d:d;
142 M4 = stripeMatrix(stencil, diags, m);
136 143
137 %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)); 144 %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));
138 145
139 M4_U=[0.5762947e7 / 0.2316384e7 -0.6374287e7 / 0.1158192e7 0.573947e6 / 0.165456e6 -0.124637e6 / 0.289548e6 0.67979e5 / 0.2316384e7 -0.60257e5 / 0.1158192e7; -0.6374287e7 / 0.1158192e7 0.30392389e8 / 0.2316384e7 -0.2735053e7 / 0.289548e6 0.273109e6 / 0.165456e6 0.83767e5 / 0.1158192e7 0.245549e6 / 0.2316384e7; 0.573947e6 / 0.165456e6 -0.2735053e7 / 0.289548e6 0.5266855e7 / 0.579096e6 -0.1099715e7 / 0.289548e6 0.869293e6 / 0.1158192e7 -0.10195e5 / 0.144774e6; -0.124637e6 / 0.289548e6 0.273109e6 / 0.165456e6 -0.1099715e7 / 0.289548e6 0.3259225e7 / 0.579096e6 -0.324229e6 / 0.72387e5 0.1847891e7 / 0.1158192e7; 0.67979e5 / 0.2316384e7 0.83767e5 / 0.1158192e7 0.869293e6 / 0.1158192e7 -0.324229e6 / 0.72387e5 0.2626501e7 / 0.330912e6 -0.7115491e7 / 0.1158192e7; -0.60257e5 / 0.1158192e7 0.245549e6 / 0.2316384e7 -0.10195e5 / 0.144774e6 0.1847891e7 / 0.1158192e7 -0.7115491e7 / 0.1158192e7 0.21383077e8 / 0.2316384e7;]; 146 M4_U=[0.5762947e7 / 0.2316384e7 -0.6374287e7 / 0.1158192e7 0.573947e6 / 0.165456e6 -0.124637e6 / 0.289548e6 0.67979e5 / 0.2316384e7 -0.60257e5 / 0.1158192e7; -0.6374287e7 / 0.1158192e7 0.30392389e8 / 0.2316384e7 -0.2735053e7 / 0.289548e6 0.273109e6 / 0.165456e6 0.83767e5 / 0.1158192e7 0.245549e6 / 0.2316384e7; 0.573947e6 / 0.165456e6 -0.2735053e7 / 0.289548e6 0.5266855e7 / 0.579096e6 -0.1099715e7 / 0.289548e6 0.869293e6 / 0.1158192e7 -0.10195e5 / 0.144774e6; -0.124637e6 / 0.289548e6 0.273109e6 / 0.165456e6 -0.1099715e7 / 0.289548e6 0.3259225e7 / 0.579096e6 -0.324229e6 / 0.72387e5 0.1847891e7 / 0.1158192e7; 0.67979e5 / 0.2316384e7 0.83767e5 / 0.1158192e7 0.869293e6 / 0.1158192e7 -0.324229e6 / 0.72387e5 0.2626501e7 / 0.330912e6 -0.7115491e7 / 0.1158192e7; -0.60257e5 / 0.1158192e7 0.245549e6 / 0.2316384e7 -0.10195e5 / 0.144774e6 0.1847891e7 / 0.1158192e7 -0.7115491e7 / 0.1158192e7 0.21383077e8 / 0.2316384e7;];
140 147
141 M4(1:6,1:6)=M4_U; 148 M4(1:6,1:6)=M4_U;
142 149
143 M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); 150 M4(m-5:m,m-5:m)=rot90( M4_U ,2 );
144 M4=M4/h^3; 151 M4=M4/h^3;
145 152
146 S3_U=[-1 3 -3 1;]/h^3; 153 S3_U=[-1 3 -3 1;]/h^3;
147 S3_1=zeros(1,m); 154 S3_1=sparse(1,m);
148 S3_1(1:4)=S3_U; 155 S3_1(1:4)=S3_U;
149 S3_m=zeros(1,m); 156 S3_m=sparse(1,m);
150 S3_m(m-3:m)=fliplr(-S3_U); 157 S3_m(m-3:m)=fliplr(-S3_U);
151 S3_1 = S3_1'; 158 S3_1 = S3_1';
152 S3_m = S3_m'; 159 S3_m = S3_m';
153 160
154 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 161 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');