Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_12.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
5 BP = 15; 5 BP = 15;
6 if(m<2*BP) 6 if(m<2*BP)
7 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 7 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
8 end 8 end
9 9
10 H=diag(ones(m,1),0); 10 H=speye(m,m);
11 H(1:15,1:15)=diag([2.880607858916397e-01,... 11 H(1:15,1:15)=diag([2.880607858916397e-01,...
12 1.560376162339675e+00,... 12 1.560376162339675e+00,...
13 2.403139445479289e-01,... 13 2.403139445479289e-01,...
14 1.531292896754208e+00,... 14 1.531292896754208e+00,...
15 1.245107365551116e+00,... 15 1.245107365551116e+00,...
27 H(m-14:m,m-14:m)=rot90(H(1:15,1:15),2); 27 H(m-14:m,m-14:m)=rot90(H(1:15,1:15),2);
28 28
29 H=H*h; 29 H=H*h;
30 HI = inv(H); 30 HI = inv(H);
31 31
32 a1 = 6/7; a2 = -15/56; a3 = 5/63; a4 = -1/56; a5 = 1/385; a6 = -1/5544; 32 diags = -6:6;
33 33 stencil = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544];
34 D1=a6*(diag(ones(m-6,1),6)-diag(ones(m-6,1),-6))+... 34 D1 = stripeMatrix(stencil, diags, m);
35 a5*(diag(ones(m-5,1),5)-diag(ones(m-5,1),-5))+...
36 a4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+...
37 a3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+...
38 a2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+...
39 a1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
40 35
41 36
42 D1(1, 1)= -1.735744761135539e+00; 37 D1(1, 1)= -1.735744761135539e+00;
43 D1(1 , 2 )= 2.307154400697361e+00; 38 D1(1 , 2 )= 2.307154400697361e+00;
44 D1(1, 3)= -8.049541548921420e-02; 39 D1(1, 3)= -8.049541548921420e-02;
279 274
280 % Wide 2nd derivative 275 % Wide 2nd derivative
281 276
282 D2=D1*D1; 277 D2=D1*D1;
283 278
284 e_1 = zeros(m,1); 279 e_1 = sparse(m,1);
285 e_1(1)= 1; 280 e_1(1)= 1;
286 e_m = zeros(m,1); 281 e_m = sparse(m,1);
287 e_m(end)= 1; 282 e_m(end)= 1;
288 283
289 S_1 = (e_1'*D1)'; 284 S_1 = (e_1'*D1)';
290 S_m = (e_m'*D1)'; 285 S_m = (e_m'*D1)';
291 286
292 Q = H*D1-(-e_1*e_1' + e_m*e_m'); 287 Q = H*D1-(-(e_1*e_1') + (e_m*e_m'));
293 M = -(H*D2-(-e_1*S_1' + e_m*S_m')); 288 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
294 end 289 end