Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6.m @ 316:203afa156f59 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Collected boundary operators.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 23:10:44 +0200
parents 9230c056a574
children c7ac7e12de8a
comparison
equal deleted inserted replaced
315:297d2cbfbe15 316:203afa156f59
13 %%% %%% 13 %%% %%%
14 %%% RR ?r dissipation) %%% 14 %%% RR ?r dissipation) %%%
15 %%% Dissipationen uppbyggd av D4: %%% 15 %%% Dissipationen uppbyggd av D4: %%%
16 %%% DI=D4*B*H*D4 %%% 16 %%% DI=D4*B*H*D4 %%%
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18
19
20 %m=10; %problemstorlek
21 %h=1/(m-1);
22
23 % Variable koefficicients are stored in vector: c, size m,
24 % with the unknown stored as c(1), c(2), ..., c_m
25 % x=1:h:m*h;x=x';
26 % c=x.^0;
27 18
28 BP = 8; 19 BP = 8;
29 if(m<2*BP) 20 if(m<2*BP)
30 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 21 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
31 end 22 end
74 S_1(1:5)=S_U; 65 S_1(1:5)=S_U;
75 S_m=sparse(1,m); 66 S_m=sparse(1,m);
76 S_m(m-4:m)=fliplr(-S_U); 67 S_m(m-4:m)=fliplr(-S_U);
77 S_1 = S_1'; 68 S_1 = S_1';
78 S_m = S_m'; 69 S_m = S_m';
79 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';
80 88
81 89
82 %DS=sparse(m,m); 90 %DS=sparse(m,m);
83 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4]; 91 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
84 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]); 92 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
88 H=h*H; 96 H=h*H;
89 HI=inv(H); 97 HI=inv(H);
90 98
91 99
92 M=sparse(m,m); 100 M=sparse(m,m);
93 e_1 = sparse(e_1);
94 e_m = sparse(e_m);
95 S_1 = sparse(S_1);
96 S_m = sparse(S_m);
97 101
98 scheme_width = 7; 102 scheme_width = 7;
99 scheme_radius = (scheme_width-1)/2; 103 scheme_radius = (scheme_width-1)/2;
100 r = (1+scheme_radius):(m-scheme_radius); 104 r = (1+scheme_radius):(m-scheme_radius);
101 105
143 147
144 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m'); 148 D2=HI*(-M-diag(c)*e_1*S_1'+diag(c)*e_m*S_m');
145 end 149 end
146 D2 = @D2_fun; 150 D2 = @D2_fun;
147 151
148 S2_U=[0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2; 152
149 S2_1=sparse(1,m);
150 S2_1(1:5)=S2_U;
151 S2_m=sparse(1,m);
152 S2_m(m-4:m)=fliplr(S2_U);
153 S2_1 = S2_1';
154 S2_m = S2_m';
155 153
156 154
157 155
158 156
159 157
187 M4(1:8,1:8) = M4_U; 185 M4(1:8,1:8) = M4_U;
188 186
189 M4(m-7:m,m-7:m) = rot90( M4_U ,2 ); 187 M4(m-7:m,m-7:m) = rot90( M4_U ,2 );
190 M4 = M4/h^3; 188 M4 = M4/h^3;
191 189
192 S3_U = [-5/2 9 -12 7 -3/2]/h^3; 190
193 S3_1 = sparse(1,m);
194 S3_1(1:5)=S3_U;
195 S3_m = sparse(1,m);
196 S3_m(m-4:m) = fliplr(-S3_U);
197 S3_1 = S3_1';
198 S3_m = S3_m';
199 191
200 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 192 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');
201 193
202 end 194 end