Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6.m @ 317:c7ac7e12de8a feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up 6th order.
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 25 Sep 2016 16:52:26 +0200
parents 203afa156f59
children 99005a80b4c2
comparison
equal deleted inserted replaced
316:203afa156f59 317:c7ac7e12de8a
19 BP = 8; 19 BP = 8;
20 if(m<2*BP) 20 if(m<2*BP)
21 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 21 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
22 end 22 end
23 23
24 24 % Norm
25 H=speye(m,m); 25 Hv = ones(m,1);
26 H(1:6,1:6)=diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]); 26 Hv(1:6) = [13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200];
27 H(m-5:m,m-5:m)=rot90(diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640,43801/43200]),2); 27 Hv(m-5:m) = rot90(Hv(1:6),2);
28 Hv = h*Hv;
29 H = spdiag(Hv, 0);
30 HI = spdiag(1./Hv, 0);
28 31
29 32
33 % Boundary operators
34 e_l = sparse(m,1);
35 e_l(1) = 1;
36 e_r = rot90(e_l, 2);
37
38 d1_l = sparse(m,1);
39 d1_l(1:5) = [-25/12, 4, -3, 4/3, -1/4]/h;
40 d1_r = -rot90(d1_l);
41
42 d2_l = sparse(m,1);
43 d2_l(1:5) = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
44 d2_r = rot90(d2_l, 2);
45
46 d3_l = sparse(m,1);
47 d3_l(1:5) = [-5/2 9 -12 7 -3/2]/h^3;
48 d3_r = -rot90(d3_l, 2);
49
50
51 % First derivtive
30 % x1=0.70127127127127; 52 % x1=0.70127127127127;
31 53
32 54
33 % 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)); 55 % 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));
34 % 56 %
53 % -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ... 75 % -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ...
54 % -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801]; 76 % -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801];
55 % D1(m-5:m,m-8:m)=rot90( -D1(1:6,1:9),2); 77 % D1(m-5:m,m-8:m)=rot90( -D1(1:6,1:9),2);
56 % D1=D1/h; 78 % D1=D1/h;
57 79
58 e_1=sparse(m,1);
59 e_1(1)=1;
60 e_m=sparse(m,1);
61 e_m(m)=1;
62 80
63 S_U=[-25/12, 4, -3, 4/3, -1/4]/h; 81 % Second derivative
64 S_1=sparse(1,m);
65 S_1(1:5)=S_U;
66 S_m=sparse(1,m);
67 S_m(m-4:m)=fliplr(-S_U);
68 S_1 = S_1';
69 S_m = S_m';
70 e_1 = sparse(e_1);
71 e_m = sparse(e_m);
72 S_1 = sparse(S_1);
73 S_m = sparse(S_m);
74 S2_U=[0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
75 S2_1=sparse(1,m);
76 S2_1(1:5)=S2_U;
77 S2_m=sparse(1,m);
78 S2_m(m-4:m)=fliplr(S2_U);
79 S2_1 = S2_1';
80 S2_m = S2_m';
81 S3_U = [-5/2 9 -12 7 -3/2]/h^3;
82 S3_1 = sparse(1,m);
83 S3_1(1:5)=S3_U;
84 S3_m = sparse(1,m);
85 S3_m(m-4:m) = fliplr(-S3_U);
86 S3_1 = S3_1';
87 S3_m = S3_m';
88
89
90 %DS=sparse(m,m);
91 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
92 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
93 %DS=diag(c)*DS/h;
94
95
96 H=h*H;
97 HI=inv(H);
98
99
100 M=sparse(m,m); 82 M=sparse(m,m);
101 83
102 scheme_width = 7; 84 scheme_width = 7;
103 scheme_radius = (scheme_width-1)/2; 85 scheme_radius = (scheme_width-1)/2;
104 r = (1+scheme_radius):(m-scheme_radius); 86 r = (1+scheme_radius):(m-scheme_radius);
148 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m'); 130 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m');
149 end 131 end
150 D2 = @D2_fun; 132 D2 = @D2_fun;
151 133
152 134
153
154
155
156
157
158 % Fourth derivative, 1th order accurate at first 8 boundary points (still 135 % Fourth derivative, 1th order accurate at first 8 boundary points (still
159 % yield 5th order convergence if stable: for example u_tt=-u_xxxx 136 % yield 5th order convergence if stable: for example u_tt=-u_xxxx
160 137 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240];
161 m4=7/240; 138 diags = -4:4;
162 m3=-2/5;
163 m2=169/60;
164 m1=-122/15;
165 m0=91/8;
166 % 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);
167 stencil = [m4,m3,m2,m1,m0,m1,m2,m3,m4];
168 d = (length(stencil)-1)/2;
169 diags = -d:d;
170 M4 = stripeMatrix(stencil, diags, m); 139 M4 = stripeMatrix(stencil, diags, m);
171
172 %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));
173 140
174 M4_U = [ 141 M4_U = [
175 0.1394226315049e13/0.367201486080e12 -0.1137054563243e13/0.114750464400e12 0.16614189027367e14/0.1836007430400e13 -0.1104821700277e13/0.306001238400e12 0.1355771086763e13/0.1836007430400e13 -0.27818686453e11/0.459001857600e12 -0.40671054239e11/0.1836007430400e13 0.5442887371e10/0.306001238400e12; 142 0.1394226315049e13/0.367201486080e12 -0.1137054563243e13/0.114750464400e12 0.16614189027367e14/0.1836007430400e13 -0.1104821700277e13/0.306001238400e12 0.1355771086763e13/0.1836007430400e13 -0.27818686453e11/0.459001857600e12 -0.40671054239e11/0.1836007430400e13 0.5442887371e10/0.306001238400e12;
176 -0.1137054563243e13/0.114750464400e12 0.70616795535409e14/0.2570410402560e13 -0.173266854731041e15/0.6426026006400e13 0.28938615291031e14/0.2570410402560e13 -0.146167361863e12/0.71400288960e11 0.2793470836571e13/0.12852052012800e14 0.6219558097e10/0.428401733760e12 -0.7313844559e10/0.166909766400e12; 143 -0.1137054563243e13/0.114750464400e12 0.70616795535409e14/0.2570410402560e13 -0.173266854731041e15/0.6426026006400e13 0.28938615291031e14/0.2570410402560e13 -0.146167361863e12/0.71400288960e11 0.2793470836571e13/0.12852052012800e14 0.6219558097e10/0.428401733760e12 -0.7313844559e10/0.166909766400e12;
177 0.16614189027367e14/0.1836007430400e13 -0.173266854731041e15/0.6426026006400e13 0.378613061504779e15/0.12852052012800e14 -0.9117069604217e13/0.642602600640e12 0.632177582849e12/0.233673672960e12 -0.1057776382577e13/0.6426026006400e13 0.443019868399e12/0.4284017337600e13 -0.3707981e7/0.2318191200e10; 144 0.16614189027367e14/0.1836007430400e13 -0.173266854731041e15/0.6426026006400e13 0.378613061504779e15/0.12852052012800e14 -0.9117069604217e13/0.642602600640e12 0.632177582849e12/0.233673672960e12 -0.1057776382577e13/0.6426026006400e13 0.443019868399e12/0.4284017337600e13 -0.3707981e7/0.2318191200e10;
181 -0.40671054239e11/0.1836007430400e13 0.6219558097e10/0.428401733760e12 0.443019868399e12/0.4284017337600e13 -0.13731270505e11/0.64260260064e11 0.310830296467e12/0.171360693504e12 -0.14432772918527e14/0.2142008668800e13 0.27102479467823e14/0.2570410402560e13 -0.1216032192203e13/0.153000619200e12; 148 -0.40671054239e11/0.1836007430400e13 0.6219558097e10/0.428401733760e12 0.443019868399e12/0.4284017337600e13 -0.13731270505e11/0.64260260064e11 0.310830296467e12/0.171360693504e12 -0.14432772918527e14/0.2142008668800e13 0.27102479467823e14/0.2570410402560e13 -0.1216032192203e13/0.153000619200e12;
182 0.5442887371e10/0.306001238400e12 -0.7313844559e10/0.166909766400e12 -0.3707981e7/0.2318191200e10 0.2933596129e10/0.40800165120e11 -0.55284274391e11/0.183600743040e12 0.58102695589e11/0.22666758400e11 -0.1216032192203e13/0.153000619200e12 0.20799922829107e14/0.1836007430400e13; 149 0.5442887371e10/0.306001238400e12 -0.7313844559e10/0.166909766400e12 -0.3707981e7/0.2318191200e10 0.2933596129e10/0.40800165120e11 -0.55284274391e11/0.183600743040e12 0.58102695589e11/0.22666758400e11 -0.1216032192203e13/0.153000619200e12 0.20799922829107e14/0.1836007430400e13;
183 ]; 150 ];
184 151
185 M4(1:8,1:8) = M4_U; 152 M4(1:8,1:8) = M4_U;
186
187 M4(m-7:m,m-7:m) = rot90( M4_U ,2 ); 153 M4(m-7:m,m-7:m) = rot90( M4_U ,2 );
188 M4 = M4/h^3; 154 M4 = M4/h^3;
189 155
190 156
191 157
192 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 158 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
193
194 end 159 end