Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_halfvariable_2.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
25 % randsplutningar f?r D1 och D2 25 % randsplutningar f?r D1 och D2
26 26
27 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm, 27 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
28 % vilket ?r n?dv?ndigt f?r stabilitet 28 % vilket ?r n?dv?ndigt f?r stabilitet
29 29
30 BP = 1; 30 BP = 4;
31 if(m<2*BP) 31 if(m<2*BP)
32 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 32 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
33 end 33 end
34 34
35 H=diag(ones(m,1),0);H(1,1)=1/2;H(m,m)=1/2; 35 H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2;
36 36
37 37
38 H=H*h; 38 H=H*h;
39 HI=inv(H); 39 HI=inv(H);
40 40
41 41
42 % First derivative SBP operator, 1st order accurate at first 6 boundary points 42 % First derivative SBP operator, 1st order accurate at first 6 boundary points
43 43
44 q1=1/2; 44 q1=1/2;
45 Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); 45 % Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
46 stencil = [-q1,0,q1];
47 d = (length(stencil)-1)/2;
48 diags = -d:d;
49 Q = stripeMatrix(stencil, diags, m);
46 50
47 %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)); 51 %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));
48 52
49 53
50 e_1=zeros(m,1);e_1(1)=1; 54 e_1=sparse(m,1);e_1(1)=1;
51 e_m=zeros(m,1);e_m(m)=1; 55 e_m=sparse(m,1);e_m(m)=1;
52 56
53 57
54 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 58 D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
55 59
56 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57 61
58 62
59 63
67 71
68 % Below for variable coefficients 72 % Below for variable coefficients
69 % Require a vector c with the koeffients 73 % Require a vector c with the koeffients
70 74
71 S_U=[-3/2 2 -1/2]/h; 75 S_U=[-3/2 2 -1/2]/h;
72 S_1=zeros(1,m); 76 S_1=sparse(1,m);
73 S_1(1:3)=S_U; 77 S_1(1:3)=S_U;
74 S_m=zeros(1,m); 78 S_m=sparse(1,m);
75 S_m(m-2:m)=fliplr(-S_U); 79 S_m(m-2:m)=fliplr(-S_U);
76 80
77 S_1 = S_1'; 81 S_1 = S_1';
78 S_m = S_m'; 82 S_m = S_m';
79 83
113 117
114 118
115 % Third derivative, 1st order accurate at first 6 boundary points 119 % Third derivative, 1st order accurate at first 6 boundary points
116 120
117 q2=1/2;q1=-1; 121 q2=1/2;q1=-1;
118 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)); 122 % 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));
123 stencil = [-q2,-q1,0,q1,q2];
124 d = (length(stencil)-1)/2;
125 diags = -d:d;
126 Q3 = stripeMatrix(stencil, diags, m);
119 127
120 %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)); 128 %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));
121 129
122 130
123 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;]; 131 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;];
124 Q3(1:4,1:4)=Q3_U; 132 Q3(1:4,1:4)=Q3_U;
125 Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) ); 133 Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 );
126 Q3=Q3/h^2; 134 Q3=Q3/h^2;
127 135
128 136
129 137
130 S2_U=[1 -2 1;]/h^2; 138 S2_U=[1 -2 1;]/h^2;
131 S2_1=zeros(1,m); 139 S2_1=sparse(1,m);
132 S2_1(1:3)=S2_U; 140 S2_1(1:3)=S2_U;
133 S2_m=zeros(1,m); 141 S2_m=sparse(1,m);
134 S2_m(m-2:m)=fliplr(S2_U); 142 S2_m(m-2:m)=fliplr(S2_U);
135 S2_1 = S2_1'; 143 S2_1 = S2_1';
136 S2_m = S2_m'; 144 S2_m = S2_m';
137 145
138 146
139 147
140 D3=HI*(Q3 - e_1*S2_1' + e_m*S2_m' +1/2*S_1*S_1' -1/2*S_m*S_m' ) ; 148 D3=HI*(Q3 - e_1*S2_1' + e_m*S2_m' +1/2*(S_1*S_1') -1/2*(S_m*S_m') ) ;
141 149
142 % Fourth derivative, 0th order accurate at first 6 boundary points (still 150 % Fourth derivative, 0th order accurate at first 6 boundary points (still
143 % yield 4th order convergence if stable: for example u_tt=-u_xxxx 151 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
144 152
145 m2=1;m1=-4;m0=6; 153 m2=1;m1=-4;m0=6;
146 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); 154 % 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);
155 stencil = [m2,m1,m0,m1,m2];
156 d = (length(stencil)-1)/2;
157 diags = -d:d;
158 M4 = stripeMatrix(stencil, diags, m);
147 159
148 %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)); 160 %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));
149 161
150 M4_U=[0.13e2 / 0.10e2 -0.12e2 / 0.5e1 0.9e1 / 0.10e2 0.1e1 / 0.5e1; -0.12e2 / 0.5e1 0.26e2 / 0.5e1 -0.16e2 / 0.5e1 0.2e1 / 0.5e1; 0.9e1 / 0.10e2 -0.16e2 / 0.5e1 0.47e2 / 0.10e2 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;]; 162 M4_U=[0.13e2 / 0.10e2 -0.12e2 / 0.5e1 0.9e1 / 0.10e2 0.1e1 / 0.5e1; -0.12e2 / 0.5e1 0.26e2 / 0.5e1 -0.16e2 / 0.5e1 0.2e1 / 0.5e1; 0.9e1 / 0.10e2 -0.16e2 / 0.5e1 0.47e2 / 0.10e2 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;];
151 163
152 164
153 M4(1:4,1:4)=M4_U; 165 M4(1:4,1:4)=M4_U;
154 166
155 M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) ); 167 M4(m-3:m,m-3:m)=rot90( M4_U ,2 );
156 M4=M4/h^3; 168 M4=M4/h^3;
157 169
158 S3_U=[-1 3 -3 1;]/h^3; 170 S3_U=[-1 3 -3 1;]/h^3;
159 S3_1=zeros(1,m); 171 S3_1=sparse(1,m);
160 S3_1(1:4)=S3_U; 172 S3_1(1:4)=S3_U;
161 S3_m=zeros(1,m); 173 S3_m=sparse(1,m);
162 S3_m(m-3:m)=fliplr(-S3_U); 174 S3_m(m-3:m)=fliplr(-S3_U);
163 S3_1 = S3_1'; 175 S3_1 = S3_1';
164 S3_m = S3_m'; 176 S3_m = S3_m';
165 177
166 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 178 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');