Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_halfvariable_6.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
24 % Variable koefficicients are stored in vector: c, size m, 24 % Variable koefficicients are stored in vector: c, size m,
25 % with the unknown stored as c(1), c(2), ..., c_m 25 % with the unknown stored as c(1), c(2), ..., c_m
26 % x=1:h:m*h;x=x'; 26 % x=1:h:m*h;x=x';
27 % c=x.^0; 27 % c=x.^0;
28 28
29 BP = 6; 29 BP = 8;
30 if(m<2*BP) 30 if(m<2*BP)
31 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 31 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
32 end 32 end
33 33
34 34
35 H=diag(ones(m,1),0); 35 H=speye(m,m);
36 H(1:6,1:6)=diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, ... 36 H(1:6,1:6)=diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, ...
37 43801/43200]); 37 43801/43200]);
38 H(m-5:m,m-5:m)=fliplr(flipud(diag([13649/43200,12013/8640, ... 38 H(m-5:m,m-5:m)=rot90(diag([13649/43200,12013/8640, ...
39 2711/4320,5359/4320,7877/8640,43801/43200]))); 39 2711/4320,5359/4320,7877/8640,43801/43200]),2);
40 40
41 41
42 x1=0.70127127127127; 42 % x1=0.70127127127127;
43 43
44 44
45 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)); 45 % 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));
46 %
47 %
48 %
49 % D1(1:6,1:9)=[-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ...
50 % 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ...
51 % x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; ...
52 % -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, ...
53 % -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, ...
54 % -34560/12013*x1+240569/120130, 0, 0, 0; ...
55 % 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, ...
56 % 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, ...
57 % 25920/2711*x1-182261/27110, 0, 0, 0; ...
58 % -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, ...
59 % -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, ...
60 % -17280/5359*x1+710473/321540, 72/5359, 0, 0; ...
61 % 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, ...
62 % 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, ...
63 % 0, 8640/7877*x1, -1296/7877, 144/7877, 0; ...
64 % -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602,...
65 % -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ...
66 % -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801];
67 % D1(m-5:m,m-8:m)=rot90( -D1(1:6,1:9),2);
68 % D1=D1/h;
46 69
47 70 e_1=sparse(m,1);e_1(1)=1;
48 71 e_m=sparse(m,1);e_m(m)=1;
49 D1(1:6,1:9)=[-21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ ...
50 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* ...
51 x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; ...
52 -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, ...
53 -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, ...
54 -34560/12013*x1+240569/120130, 0, 0, 0; ...
55 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, ...
56 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, ...
57 25920/2711*x1-182261/27110, 0, 0, 0; ...
58 -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, ...
59 -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, ...
60 -17280/5359*x1+710473/321540, 72/5359, 0, 0; ...
61 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, ...
62 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, ...
63 0, 8640/7877*x1, -1296/7877, 144/7877, 0; ...
64 -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602,...
65 -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, ...
66 -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801];
67 D1(m-5:m,m-8:m)=flipud( fliplr(-D1(1:6,1:9)));
68 D1=D1/h;
69
70 e_1=zeros(m,1);e_1(1)=1;
71 e_m=zeros(m,1);e_m(m)=1;
72 72
73 S_U=[-25/12, 4, -3, 4/3, -1/4]/h; 73 S_U=[-25/12, 4, -3, 4/3, -1/4]/h;
74 S_1=zeros(1,m); 74 S_1=sparse(1,m);
75 S_1(1:5)=S_U; 75 S_1(1:5)=S_U;
76 S_m=zeros(1,m); 76 S_m=sparse(1,m);
77 S_m(m-4:m)=fliplr(-S_U); 77 S_m(m-4:m)=fliplr(-S_U);
78 S_1 = S_1'; 78 S_1 = S_1';
79 S_m = S_m'; 79 S_m = S_m';
80 80
81 81
82 82
83 %DS=zeros(m,m); 83 %DS=sparse(m,m);
84 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4]; 84 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
85 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]); 85 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
86 %DS=diag(c)*DS/h; 86 %DS=diag(c)*DS/h;
87 87
88 88
125 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m'); 125 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m');
126 end 126 end
127 D2 = @D2_fun; 127 D2 = @D2_fun;
128 128
129 S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2; 129 S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2;
130 S2_1=zeros(1,m); 130 S2_1=sparse(1,m);
131 S2_1(1:5)=S2_U; 131 S2_1(1:5)=S2_U;
132 S2_m=zeros(1,m); 132 S2_m=sparse(1,m);
133 S2_m(m-4:m)=fliplr(S2_U); 133 S2_m(m-4:m)=fliplr(S2_U);
134 S2_1 = S2_1'; 134 S2_1 = S2_1';
135 S2_m = S2_m'; 135 S2_m = S2_m';
136 136
137 137
140 140
141 % Fourth derivative, 1th order accurate at first 8 boundary points (still 141 % Fourth derivative, 1th order accurate at first 8 boundary points (still
142 % yield 5th order convergence if stable: for example u_tt=-u_xxxx 142 % yield 5th order convergence if stable: for example u_tt=-u_xxxx
143 143
144 m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8; 144 m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8;
145 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); 145 % 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);
146 stencil = [m4,m3,m2,m1,m0,m1,m2,m3,m4];
147 d = (length(stencil)-1)/2;
148 diags = -d:d;
149 M4 = stripeMatrix(stencil, diags, m);
146 150
147 %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)); 151 %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));
148 152
149 M4_U=[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; -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; 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; -0.1104821700277e13 / 0.306001238400e12 0.28938615291031e14 / 0.2570410402560e13 -0.9117069604217e13 / 0.642602600640e12 0.5029150721885e13 / 0.514082080512e12 -0.5209119714341e13 / 0.1285205201280e13 0.12235427457469e14 / 0.12852052012800e14 -0.13731270505e11 / 0.64260260064e11 0.2933596129e10 / 0.40800165120e11; 0.1355771086763e13 / 0.1836007430400e13 -0.146167361863e12 / 0.71400288960e11 0.632177582849e12 / 0.233673672960e12 -0.5209119714341e13 / 0.1285205201280e13 0.14871726798559e14 / 0.2570410402560e13 -0.7504337615347e13 / 0.1606506501600e13 0.310830296467e12 / 0.171360693504e12 -0.55284274391e11 / 0.183600743040e12; -0.27818686453e11 / 0.459001857600e12 0.2793470836571e13 / 0.12852052012800e14 -0.1057776382577e13 / 0.6426026006400e13 0.12235427457469e14 / 0.12852052012800e14 -0.7504337615347e13 / 0.1606506501600e13 0.106318657014853e15 / 0.12852052012800e14 -0.14432772918527e14 / 0.2142008668800e13 0.58102695589e11 / 0.22666758400e11; -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; 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;]; 153 M4_U=[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; -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; 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; -0.1104821700277e13 / 0.306001238400e12 0.28938615291031e14 / 0.2570410402560e13 -0.9117069604217e13 / 0.642602600640e12 0.5029150721885e13 / 0.514082080512e12 -0.5209119714341e13 / 0.1285205201280e13 0.12235427457469e14 / 0.12852052012800e14 -0.13731270505e11 / 0.64260260064e11 0.2933596129e10 / 0.40800165120e11; 0.1355771086763e13 / 0.1836007430400e13 -0.146167361863e12 / 0.71400288960e11 0.632177582849e12 / 0.233673672960e12 -0.5209119714341e13 / 0.1285205201280e13 0.14871726798559e14 / 0.2570410402560e13 -0.7504337615347e13 / 0.1606506501600e13 0.310830296467e12 / 0.171360693504e12 -0.55284274391e11 / 0.183600743040e12; -0.27818686453e11 / 0.459001857600e12 0.2793470836571e13 / 0.12852052012800e14 -0.1057776382577e13 / 0.6426026006400e13 0.12235427457469e14 / 0.12852052012800e14 -0.7504337615347e13 / 0.1606506501600e13 0.106318657014853e15 / 0.12852052012800e14 -0.14432772918527e14 / 0.2142008668800e13 0.58102695589e11 / 0.22666758400e11; -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; 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;];
150 154
151 M4(1:8,1:8)=M4_U; 155 M4(1:8,1:8)=M4_U;
152 156
153 M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) ); 157 M4(m-7:m,m-7:m)=rot90( M4_U ,2 );
154 M4=M4/h^3; 158 M4=M4/h^3;
155 159
156 S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3; 160 S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3;
157 S3_1=zeros(1,m); 161 S3_1=sparse(1,m);
158 S3_1(1:5)=S3_U; 162 S3_1(1:5)=S3_U;
159 S3_m=zeros(1,m); 163 S3_m=sparse(1,m);
160 S3_m(m-4:m)=fliplr(-S3_U); 164 S3_m(m-4:m)=fliplr(-S3_U);
161 S3_1 = S3_1'; 165 S3_1 = S3_1';
162 S3_m = S3_m'; 166 S3_m = S3_m';
163 167
164 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 168 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');