Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_variable_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 f18142c1530b
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
5 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 5 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
6 end 6 end
7 7
8 N = m; 8 N = m;
9 9
10 H = speye(N); 10 H = spspeye(N);
11 H(1,1) = 17/48; H(2,2) = 59/48; H(3,3) = 43/48; H(4,4) = 49/48; 11 H(1,1) = 17/48; H(2,2) = 59/48; H(3,3) = 43/48; H(4,4) = 49/48;
12 H(N,N) = 17/48; H(N-1,N-1) = 59/48; H(N-2,N-2) = 43/48; H(N-3,N-3) = 49/48; 12 H(N,N) = 17/48; H(N-1,N-1) = 59/48; H(N-2,N-2) = 43/48; H(N-3,N-3) = 49/48;
13 H = h*H; 13 H = h*H;
14 14
15 HI = inv(H); 15 HI = inv(H);
36 Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; 36 Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
37 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_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;];
38 Q(1:4,1:4)=Q_U; 38 Q(1:4,1:4)=Q_U;
39 Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); 39 Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
40 40
41 D1=HI*(Q-1/2*e_l*e_l'+1/2*e_r*e_r') ; 41 D1=HI*(Q-1/2*(e_l*e_l')+1/2*(e_r*e_r')) ;
42 42
43 function D2 = D2_fun(c) 43 function D2 = D2_fun(c)
44 M = 78+(N-12)*5; 44 M = 78+(N-12)*5;
45 %h = 1/(N-1); 45 %h = 1/(N-1);
46 46
60 48/(59)*(-0.8673e4 / 0.2904112e7 * c(N-2) - 0.8673e4 / 0.2904112e7 * c(N-1) + 0.8673e4 / 0.1452056e7 * c(N-3)) 48/(59)*(-0.1328188692663e13 / 0.37594290333616e14 * c(N-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(N-1) - 0.8673e4 / 0.363014e6 * c(N-3)) 48/(59)*(0.260297319232891e15 / 0.2556411742685888e16 * c(N-2) - 0.59e2 / 0.1088e4 * c(N) - 0.106641839640553e15 / 0.1278205871342944e16 * c(N-1) + 0.26019e5 / 0.726028e6 * c(N-3)) 48/(59)*(-0.59e2 / 0.408e3 * c(N) - 0.29294615794607e14 / 0.29725717938208e14 * c(N-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(N-1) - 0.8673e4 / 0.363014e6 * c(N-3)) 48/(59)*(0.9258282831623875e16 / 0.7669235228057664e16 * c(N-2) + 0.3481e4 / 0.3264e4 * c(N) + 0.228389721191751e15 / 0.1278205871342944e16 * c(N-1) + 0.8673e4 / 0.1452056e7 * c(N-3)) 48/(59)*(-0.6025413881e10 / 0.21126554976e11 * c(N-2) - 0.59e2 / 0.68e2 * c(N) - 0.537416663e9 / 0.7042184992e10 * c(N-1));... 60 48/(59)*(-0.8673e4 / 0.2904112e7 * c(N-2) - 0.8673e4 / 0.2904112e7 * c(N-1) + 0.8673e4 / 0.1452056e7 * c(N-3)) 48/(59)*(-0.1328188692663e13 / 0.37594290333616e14 * c(N-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(N-1) - 0.8673e4 / 0.363014e6 * c(N-3)) 48/(59)*(0.260297319232891e15 / 0.2556411742685888e16 * c(N-2) - 0.59e2 / 0.1088e4 * c(N) - 0.106641839640553e15 / 0.1278205871342944e16 * c(N-1) + 0.26019e5 / 0.726028e6 * c(N-3)) 48/(59)*(-0.59e2 / 0.408e3 * c(N) - 0.29294615794607e14 / 0.29725717938208e14 * c(N-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(N-1) - 0.8673e4 / 0.363014e6 * c(N-3)) 48/(59)*(0.9258282831623875e16 / 0.7669235228057664e16 * c(N-2) + 0.3481e4 / 0.3264e4 * c(N) + 0.228389721191751e15 / 0.1278205871342944e16 * c(N-1) + 0.8673e4 / 0.1452056e7 * c(N-3)) 48/(59)*(-0.6025413881e10 / 0.21126554976e11 * c(N-2) - 0.59e2 / 0.68e2 * c(N) - 0.537416663e9 / 0.7042184992e10 * c(N-1));...
61 48/(17)*(-c(N-3) / 0.392e3 + c(N-2) / 0.784e3 + c(N-1) / 0.784e3) 48/(17)*(c(N-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(N-2) - 0.256702175e9 / 0.10149031312e11 * c(N-1)) 48/(17)*(-0.1244724001e10 / 0.21126554976e11 * c(N-2) + 0.3e1 / 0.68e2 * c(N) + 0.752806667e9 / 0.21126554976e11 * c(N-1)) 48/(17)*(-0.59e2 / 0.3136e4 * c(N-3) - 0.13249937023e11 / 0.48148892736e11 * c(N-1) + 0.2e1 / 0.17e2 * c(N) + 0.2083938599e10 / 0.8024815456e10 * c(N-2)) 48/(17)*(-0.6025413881e10 / 0.21126554976e11 * c(N-2) - 0.59e2 / 0.68e2 * c(N) - 0.537416663e9 / 0.7042184992e10 * c(N-1)) 48/(17)*(0.3e1 / 0.3136e4 * c(N-3) + 0.27010400129e11 / 0.345067064608e12 * c(N-2) + 0.234566387291e12 / 0.690134129216e12 * c(N-1) + 0.12e2 / 0.17e2 * c(N))]; 61 48/(17)*(-c(N-3) / 0.392e3 + c(N-2) / 0.784e3 + c(N-1) / 0.784e3) 48/(17)*(c(N-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(N-2) - 0.256702175e9 / 0.10149031312e11 * c(N-1)) 48/(17)*(-0.1244724001e10 / 0.21126554976e11 * c(N-2) + 0.3e1 / 0.68e2 * c(N) + 0.752806667e9 / 0.21126554976e11 * c(N-1)) 48/(17)*(-0.59e2 / 0.3136e4 * c(N-3) - 0.13249937023e11 / 0.48148892736e11 * c(N-1) + 0.2e1 / 0.17e2 * c(N) + 0.2083938599e10 / 0.8024815456e10 * c(N-2)) 48/(17)*(-0.6025413881e10 / 0.21126554976e11 * c(N-2) - 0.59e2 / 0.68e2 * c(N) - 0.537416663e9 / 0.7042184992e10 * c(N-1)) 48/(17)*(0.3e1 / 0.3136e4 * c(N-3) + 0.27010400129e11 / 0.345067064608e12 * c(N-2) + 0.234566387291e12 / 0.690134129216e12 * c(N-1) + 0.12e2 / 0.17e2 * c(N))];
62 62
63 63
64 64
65 R = zeros(M,1); 65 R = sparse(M,1);
66 R(1:24) = reshape(U(1:4,:)',24,1); 66 R(1:24) = reshape(U(1:4,:)',24,1);
67 R(25:30) = U(5,:); 67 R(25:30) = U(5,:);
68 R(31) = -c(5+1) / 0.6e1 + c(5) / 0.8e1 + c(5+2) / 0.8e1; 68 R(31) = -c(5+1) / 0.6e1 + c(5) / 0.8e1 + c(5+2) / 0.8e1;
69 R(32:37) = U(6,:); 69 R(32:37) = U(6,:);
70 R(38:39) = [-c(6-1) / 0.6e1 - c(6+2) / 0.6e1 - c(6) / 0.2e1 - c(6+1) / 0.2e1;... 70 R(38:39) = [-c(6-1) / 0.6e1 - c(6+2) / 0.6e1 - c(6) / 0.2e1 - c(6+1) / 0.2e1;...
90 D2(1:4) = -48/17/h/h*[c(1)*(-11/6);c(1)*3;c(1)*(-3/2);c(1)*1/3] + D2(1:4); 90 D2(1:4) = -48/17/h/h*[c(1)*(-11/6);c(1)*3;c(1)*(-3/2);c(1)*1/3] + D2(1:4);
91 D2(M-3:M) = -48/17/h/h*[c(N)*1/3;c(N)*(-3/2);c(N)*3;c(N)*(-11/6)] + D2(M-3:M); 91 D2(M-3:M) = -48/17/h/h*[c(N)*1/3;c(N)*(-3/2);c(N)*3;c(N)*(-11/6)] + D2(M-3:M);
92 92
93 93
94 94
95 BS = sparse(N,N); 95 % BS = sparse(N,N);
96 BS(1,1:4) = -c(1)*1/h*[(-11/6);3;(-3/2);1/3]; 96 % BS(1,1:4) = -c(1)*1/h*[(-11/6);3;(-3/2);1/3];
97 BS(N,N-3:N) = c(N)*1/h*[(-1/3);3/2;(-3);11/6]; 97 % BS(N,N-3:N) = c(N)*1/h*[(-1/3);3/2;(-3);11/6];
98 BS = sparse(BS); 98 % BS = sparse(BS);
99 99
100 % %%Row and column indices%% 100 % %%Row and column indices%%
101 M = 78+(N-12)*5; 101 M = 78+(N-12)*5;
102 rows = [kron([1;2;3;4],ones(6,1));... 102 rows = [kron([1;2;3;4],ones(6,1));...
103 5*ones(7,1);... 103 5*ones(7,1);...
105 kron((7:N-6)',ones(5,1));... 105 kron((7:N-6)',ones(5,1));...
106 (N-5)*ones(8,1);... 106 (N-5)*ones(8,1);...
107 (N-4)*ones(7,1);... 107 (N-4)*ones(7,1);...
108 kron([N-3;N-2;N-1;N],ones(6,1))]; 108 kron([N-3;N-2;N-1;N],ones(6,1))];
109 109
110 cols = zeros(M,1); 110 cols = sparse(M,1);
111 cols(1:24) = kron(ones(4,1),[1;2;3;4;5;6]); 111 cols(1:24) = kron(ones(4,1),[1;2;3;4;5;6]);
112 cols(25:39) = [(1:7)';(1:8)']; 112 cols(25:39) = [(1:7)';(1:8)'];
113 cols(M-23:M) = kron(ones(4,1),[N-5;N-4;N-3;N-2;N-1;N]); 113 cols(M-23:M) = kron(ones(4,1),[N-5;N-4;N-3;N-2;N-1;N]);
114 cols(M-38:M-24) = [(N-7:N)';(N-6:N)']; 114 cols(M-38:M-24) = [(N-7:N)';(N-6:N)'];
115 for i=7:N-6 115 for i=7:N-6