Mercurial > repos > public > sbplib

comparison +util/calc_borrowing.m @ 327:d24869abc7cd feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Calculated borrowin for D4Lonely.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 27 Sep 2016 09:46:58 +0200
parents 8b4993d53663
children
comparison
equal deleted inserted replaced
326:b19e142fcae1 327:d24869abc7cd
1 1 function calc_borrowing(m, h)
2 m = 100; 2 default_arg('m',100);
3 h = 1; 3 default_arg('h',1);
4 4
5 5 operators = {
6 %% 4th order non-compatible 6 {
7 [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4(m,h); 7 'd4_lonely', getM4_lonely, {
8 S1 = S_1*S_1' + S_m*S_m'; 8 {4, 'min_boundary_points'},
9 S2 = S2_1*S2_1' + S2_m*S2_m'; 9 {6, 'min_boundary_points'},
10 S3 = S3_1*S3_1' + S3_m*S3_m'; 10 {6, '2'},
11 11 {6, '3'},
12 alpha_I = util.matrixborrow(M4, h^-1*S1 ); 12 {8, 'min_boundary_points'},
13 alpha_II = util.matrixborrow(M4, h*S2 ); 13 {8, 'higher_boundary_order'},
14 alpha_III = util.matrixborrow(M4, h^3*S3); 14 }
15 fprintf('4th order non-compatible\n') 15 }, {
16 fprintf('alpha_I1: %.10f\n',alpha_I) 16 'd4_variable', {
17 fprintf('alpha_II: %.10f\n',alpha_II) 17 {2},
18 fprintf('alpha_III: %.10f\n',alpha_III) 18 {4},
19 fprintf('\n') 19 {6},
20 20 }
21 21 }
22 %% 6th order non-compatible 22 % BORKEN BAD IDEA
23 [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6(m,h); 23 }
24 S1 = S_1*S_1' + S_m*S_m'; 24
25 S2 = S2_1*S2_1' + S2_m*S2_m'; 25
26 S3 = S3_1*S3_1' + S3_m*S3_m'; 26 for i = 1:operators
27 27 baseName = operators{i}{1};
28 alpha_II = util.matrixborrow(M4, h*S2 ); 28 postFixes = operators{i}{2};
29 alpha_III = util.matrixborrow(M4, h^3*S3); 29 for pf = postFixes
30 fprintf('6th order non-compatible\n') 30 [a2, a3] = borrowFromD4(m, h, l{:});
31 fprintf('alpha_II: %.10f\n',alpha_II) 31 end
32 fprintf('alpha_III: %.10f\n',alpha_III) 32 end
33 fprintf('\n') 33
34 34
35 35
36 %% 2nd order compatible 36 lonely = {
37 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible2(m,h); 37 {4, 'min_boundary_points'},
38 S1 = S_1*S_1' + S_m*S_m'; 38 {6, 'min_boundary_points'},
39 S2 = S2_1*S2_1' + S2_m*S2_m'; 39 {6, '2'},
40 S3 = S3_1*S3_1' + S3_m*S3_m'; 40 {6, '3'},
41 41 {8, 'min_boundary_points'},
42 alpha_II = util.matrixborrow(M4, h*S2 ); 42 {8, 'higher_boundary_order'},
43 alpha_III = util.matrixborrow(M4, h^3*S3); 43 };
44 fprintf('2nd order compatible\n') 44
45 fprintf('alpha_II: %.10f\n',alpha_II) 45 for i = 1:length(lonely)
46 fprintf('alpha_III: %.10f\n',alpha_III) 46 l = lonely{i};
47 fprintf('\n') 47 [a2, a3] = d4_lonely(m, h, l{:});
48 48 fprintf('d4_lonely %d %s\n', l{:})
49 49 fprintf('\t alpha_II = %f\n', a2)
50 %% 4th order compatible 50 fprintf('\t alpha_III = %f\n', a3)
51 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible4(m,h); 51 fprintf('\n')
52 S1 = S_1*S_1' + S_m*S_m'; 52 end
53 S2 = S2_1*S2_1' + S2_m*S2_m'; 53
54 S3 = S3_1*S3_1' + S3_m*S3_m'; 54 variable = {
55 55 {2},
56 alpha_II = util.matrixborrow(M4, h*S2 ); 56 {4},
57 alpha_III = util.matrixborrow(M4, h^3*S3); 57 {6},
58 fprintf('4th order compatible\n') 58 };
59 fprintf('alpha_II: %.10f\n',alpha_II) 59
60 fprintf('alpha_III: %.10f\n',alpha_III) 60 for i = 1:length(variable)
61 fprintf('\n') 61 l = variable{i};
62 62 [a2, a3] = d4_variable(m, h, l{:});
63 %% 6th order compatible 63 fprintf('d4_variable %d\n', l{:})
64 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible6(m,h); 64 fprintf('\t alpha_II = %f\n', a2)
65 S1 = S_1*S_1' + S_m*S_m'; 65 fprintf('\t alpha_III = %f\n', a3)
66 S2 = S2_1*S2_1' + S2_m*S2_m'; 66 fprintf('\n')
67 S3 = S3_1*S3_1' + S3_m*S3_m'; 67 end
68 68
69 alpha_II = util.matrixborrow(M4, h*S2 ); 69
70 alpha_III = util.matrixborrow(M4, h^3*S3); 70 %% 4th order non-compatible
71 fprintf('6th order compatible\n') 71 [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4(m,h);
72 fprintf('alpha_II: %.10f\n',alpha_II) 72 S1 = S_1*S_1' + S_m*S_m';
73 fprintf('alpha_III: %.10f\n',alpha_III) 73 S2 = S2_1*S2_1' + S2_m*S2_m';
74 fprintf('\n') 74 S3 = S3_1*S3_1' + S3_m*S3_m';
75 75
76 76 alpha_I = util.matrixborrow(M4, h^-1*S1 );
77 77 alpha_II = util.matrixborrow(M4, h*S2 );
78 78 alpha_III = util.matrixborrow(M4, h^3*S3);
79 79 fprintf('4th order non-compatible\n')
80 % Ordinary 80 fprintf('alpha_I1: %.10f\n',alpha_I)
81 81 fprintf('alpha_II: %.10f\n',alpha_II)
82 for order = [2 4 6 8 10] 82 fprintf('alpha_III: %.10f\n',alpha_III)
83 op = sbp.Ordinary(m,h, order); 83 fprintf('\n')
84 84
85 S_1 = op.boundary.S_1; 85
86 S_m = op.boundary.S_m; 86 %% 6th order non-compatible
87 87 [H, HI, D1, D2, D3, D4, e_1, e_m, M, M4,Q, Q3, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6(m,h);
88 M = op.norms.M; 88 S1 = S_1*S_1' + S_m*S_m';
89 89 S2 = S2_1*S2_1' + S2_m*S2_m';
90 S1 = S_1*S_1' + S_m*S_m'; 90 S3 = S3_1*S3_1' + S3_m*S3_m';
91 alpha = util.matrixborrow(M, h*S1); 91
92 fprintf('%dth order Ordinary\n', order) 92 alpha_II = util.matrixborrow(M4, h*S2 );
93 fprintf('alpha: %.10f\n', alpha) 93 alpha_III = util.matrixborrow(M4, h^3*S3);
94 fprintf('\n') 94 fprintf('6th order non-compatible\n')
95 end 95 fprintf('alpha_II: %.10f\n',alpha_II)
96 96 fprintf('alpha_III: %.10f\n',alpha_III)
97 97 fprintf('\n')
98
99
100 %% 2nd order compatible
101 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible2(m,h);
102 S1 = S_1*S_1' + S_m*S_m';
103 S2 = S2_1*S2_1' + S2_m*S2_m';
104 S3 = S3_1*S3_1' + S3_m*S3_m';
105
106 alpha_II = util.matrixborrow(M4, h*S2 );
107 alpha_III = util.matrixborrow(M4, h^3*S3);
108 fprintf('2nd order compatible\n')
109 fprintf('alpha_II: %.10f\n',alpha_II)
110 fprintf('alpha_III: %.10f\n',alpha_III)
111 fprintf('\n')
112
113
114 %% 4th order compatible
115 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible4(m,h);
116 S1 = S_1*S_1' + S_m*S_m';
117 S2 = S2_1*S2_1' + S2_m*S2_m';
118 S3 = S3_1*S3_1' + S3_m*S3_m';
119
120 alpha_II = util.matrixborrow(M4, h*S2 );
121 alpha_III = util.matrixborrow(M4, h^3*S3);
122 fprintf('4th order compatible\n')
123 fprintf('alpha_II: %.10f\n',alpha_II)
124 fprintf('alpha_III: %.10f\n',alpha_III)
125 fprintf('\n')
126
127 %% 6th order compatible
128 [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible6(m,h);
129 S1 = S_1*S_1' + S_m*S_m';
130 S2 = S2_1*S2_1' + S2_m*S2_m';
131 S3 = S3_1*S3_1' + S3_m*S3_m';
132
133 alpha_II = util.matrixborrow(M4, h*S2 );
134 alpha_III = util.matrixborrow(M4, h^3*S3);
135 fprintf('6th order compatible\n')
136 fprintf('alpha_II: %.10f\n',alpha_II)
137 fprintf('alpha_III: %.10f\n',alpha_III)
138 fprintf('\n')
139
140
141
142
143
144 % Ordinary
145
146 for order = [2 4 6 8 10]
147 op = sbp.Ordinary(m,h, order);
148
149 S_1 = op.boundary.S_1;
150 S_m = op.boundary.S_m;
151
152 M = op.norms.M;
153
154 S1 = S_1*S_1' + S_m*S_m';
155 alpha = util.matrixborrow(M, h*S1);
156 fprintf('%dth order Ordinary\n', order)
157 fprintf('alpha: %.10f\n', alpha)
158 fprintf('\n')
159 end
160
161
162
163
164 end
165
166 function [alpha_II, alpha_III] = d4_lonely(m, h, order, modifier)
167 default_arg('modifier', [])
168 func = sprintf('sbp.implementations.d4_lonely_%d', order);
169 if ~isempty(modifier)
170 func = sprintf('%s_%s', func, modifier);
171 end
172 funcCall = sprintf('%s(%s,%s)', func, toString(m), toString(h));
173 [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall);
174
175 d2d2 = d2_l*d2_l' + d2_r*d2_r';
176 alpha_II = util.matrixborrow(M4, h*d2d2);
177
178 d3d3 = d3_l*d3_l' + d3_r*d3_r';
179 alpha_III = util.matrixborrow(M4, h^3*d3d3);
180 end
181
182 function [alpha_II, alpha_III] = d4_variable(m, h, order)
183 default_arg('modifier', [])
184 func = sprintf('sbp.implementations.d4_variable_%d', order);
185
186 funcCall = sprintf('%s(%s,%s)', func, toString(m), toString(h));
187 [H, HI, D1, D2, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall);
188
189 d2d2 = d2_l*d2_l' + d2_r*d2_r';
190 alpha_II = util.matrixborrow(M4, h*d2d2);
191
192 d3d3 = d3_l*d3_l' + d3_r*d3_r';
193 alpha_III = util.matrixborrow(M4, h^3*d3d3);
194 end
195
196 function [d2_l, d2_r, d3_l, d3_r, M4] = getM4_lonely(m, h, order, modifier)
197 fStr = getFunctionCallStr('d4_lonely', {order, modifier}, {m ,h});
198 [H, HI, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = eval(funcCall);
199 end
200
201
202 % Calculates the borrowing constants for a D4 operator.
203 % getM4 is a function handle on the form
204 % [d2_l, d2_r, d3_l, d3_r, M4] = getM4(m,h)
205 function [a2, a3] = borrowFromD4(m, h, getM4)
206 [d2_l, d2_r, d3_l, d3_r, M4] = getM4(m, h);
207
208 d2d2 = d2_l*d2_l' + d2_r*d2_r';
209 a2 = util.matrixborrow(M4, h*d2d2);
210
211 d3d3 = d3_l*d3_l' + d3_r*d3_r';
212 a3 = util.matrixborrow(M4, h^3*d3d3);
213 end
214
215
216 function funcCallStr = getFunctionCallStr(baseName, postFix, parameters)
217 default_arg('postFix', [])
218 default_arg('parameters', [])
219
220 funcCallStr = sprintf('sbp.implementations.%s', baseName);
221
222 for i = 1:length(postFix)
223 if ischar(postFix{i})
224 funcCallStr = [funcCallStr '_' postFix{i}];
225 else
226 funcCallStr = [funcCallStr '_' toString(postFix{i})];
227 end
228 end
229
230 if isempty(parameters)
231 return
232 end
233
234 funcCallStr = [funcCallStr '(' toString(parameters{1})];
235
236 for i = 2:length(parameters)
237 funcCallStr = [funcCallStr ', ' toString(parameters{i})];
238 end
239
240 funcCallStr = [funcCallStr ')';
241 end