Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 647:46c40711830f feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add test for diracDiscr.m
author Martin Almquist <malmquist@stanford.edu>
date Tue, 14 Nov 2017 15:35:28 -0800
parents
children 4ee7d15bd8e6
comparison
equal deleted inserted replaced
646:0990765e3e4d 647:46c40711830f
1 function tests = diracDiscrTest()
2 tests = functiontests(localfunctions);
3 end
4
5 function testLeftGP(testCase)
6
7 orders = [2, 4, 6];
8 mom_conds = orders;
9
10 for o = 1:length(orders)
11 order = orders(o);
12 mom_cond = mom_conds(o);
13 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
14
15 % Test left boundary grid points
16 x0s = xl + [0, h, 2*h];
17
18 for j = 1:length(fs)
19 f = fs{j};
20 fx = f(x);
21 for i = 1:length(x0s)
22 x0 = x0s(i);
23 delta = diracDiscr(x0, x, mom_cond, 0, H);
24 integral = delta'*H*fx;
25 err = abs(integral - f(x0));
26 testCase.verifyLessThan(err, 1e-12);
27 end
28 end
29 end
30 end
31
32 function testLeftRandom(testCase)
33
34 orders = [2, 4, 6];
35 mom_conds = orders;
36
37 for o = 1:length(orders)
38 order = orders(o);
39 mom_cond = mom_conds(o);
40 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
41
42 % Test random points near left boundary
43 x0s = xl + 2*h*rand(1,10);
44
45 for j = 1:length(fs)
46 f = fs{j};
47 fx = f(x);
48 for i = 1:length(x0s)
49 x0 = x0s(i);
50 delta = diracDiscr(x0, x, mom_cond, 0, H);
51 integral = delta'*H*fx;
52 err = abs(integral - f(x0));
53 testCase.verifyLessThan(err, 1e-12);
54 end
55 end
56 end
57 end
58
59 function testRightGP(testCase)
60
61 orders = [2, 4, 6];
62 mom_conds = orders;
63
64 for o = 1:length(orders)
65 order = orders(o);
66 mom_cond = mom_conds(o);
67 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
68
69 % Test right boundary grid points
70 x0s = xr-[0, h, 2*h];
71
72 for j = 1:length(fs)
73 f = fs{j};
74 fx = f(x);
75 for i = 1:length(x0s)
76 x0 = x0s(i);
77 delta = diracDiscr(x0, x, mom_cond, 0, H);
78 integral = delta'*H*fx;
79 err = abs(integral - f(x0));
80 testCase.verifyLessThan(err, 1e-12);
81 end
82 end
83 end
84 end
85
86 function testRightRandom(testCase)
87
88 orders = [2, 4, 6];
89 mom_conds = orders;
90
91 for o = 1:length(orders)
92 order = orders(o);
93 mom_cond = mom_conds(o);
94 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
95
96 % Test random points near right boundary
97 x0s = xr - 2*h*rand(1,10);
98
99 for j = 1:length(fs)
100 f = fs{j};
101 fx = f(x);
102 for i = 1:length(x0s)
103 x0 = x0s(i);
104 delta = diracDiscr(x0, x, mom_cond, 0, H);
105 integral = delta'*H*fx;
106 err = abs(integral - f(x0));
107 testCase.verifyLessThan(err, 1e-12);
108 end
109 end
110 end
111 end
112
113 function testInteriorGP(testCase)
114
115 orders = [2, 4, 6];
116 mom_conds = orders;
117
118 for o = 1:length(orders)
119 order = orders(o);
120 mom_cond = mom_conds(o);
121 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
122
123 % Test interior grid points
124 m_half = round(m/2);
125 x0s = xl + (m_half-1:m_half+1)*h;
126
127 for j = 1:length(fs)
128 f = fs{j};
129 fx = f(x);
130 for i = 1:length(x0s)
131 x0 = x0s(i);
132 delta = diracDiscr(x0, x, mom_cond, 0, H);
133 integral = delta'*H*fx;
134 err = abs(integral - f(x0));
135 testCase.verifyLessThan(err, 1e-12);
136 end
137 end
138 end
139 end
140
141 function testInteriorRandom(testCase)
142
143 orders = [2, 4, 6];
144 mom_conds = orders;
145
146 for o = 1:length(orders)
147 order = orders(o);
148 mom_cond = mom_conds(o);
149 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
150
151 % Test random points in interior
152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
153
154 for j = 1:length(fs)
155 f = fs{j};
156 fx = f(x);
157 for i = 1:length(x0s)
158 x0 = x0s(i);
159 delta = diracDiscr(x0, x, mom_cond, 0, H);
160 integral = delta'*H*fx;
161 err = abs(integral - f(x0));
162 testCase.verifyLessThan(err, 1e-12);
163 end
164 end
165 end
166 end
167
168 % x0 outside grid should yield 0 integral!
169 function testX0OutsideGrid(testCase)
170
171 orders = [2, 4, 6];
172 mom_conds = orders;
173
174 for o = 1:length(orders)
175 order = orders(o);
176 mom_cond = mom_conds(o);
177 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
178
179 % Test points outisde grid
180 x0s = [xl-1.1*h, xr+1.1*h];
181
182 for j = 1:length(fs)
183 f = fs{j};
184 fx = f(x);
185 for i = 1:length(x0s)
186 x0 = x0s(i);
187 delta = diracDiscr(x0, x, mom_cond, 0, H);
188 integral = delta'*H*fx;
189 err = abs(integral - 0);
190 testCase.verifyLessThan(err, 1e-12);
191 end
192 end
193 end
194 end
195
196 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
197
198 % Grid
199 xl = -3;
200 xr = 2;
201 m = 21;
202 h = (xr-xl)/(m-1);
203 g = grid.equidistant(m, {xl, xr});
204 x = g.points();
205
206 % Quadrature
207 ops = sbp.D2Standard(m, {xl, xr}, order);
208 H = ops.H;
209
210 % Moment conditions
211 fs = cell(mom_cond,1);
212 for p = 0:mom_cond-1
213 fs{p+1} = @(x) x.^p;
214 end
215
216 end
217
218