Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 651:4ee7d15bd8e6 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix roundoff bug in diracDiscr and improve test.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 15 Nov 2017 14:58:13 -0800
parents 1bdbe026abbc
children 5264ce57b573
comparison
equal deleted inserted replaced
650:8e55298657b9 651:4ee7d15bd8e6
1 function ret = diracDiscr(x_0in , x , m_order, s_order, H) 1 function ret = diracDiscr(x_0in , x , m_order, s_order, H)
2
3 m = length(x);
2 4
3 % Return zeros if x0 is outside grid 5 % Return zeros if x0 is outside grid
4 if(x_0in < x(1) || x_0in > x(end) ) 6 if(x_0in < x(1) || x_0in > x(end) )
5 7
6 ret = zeros(size(x)); 8 ret = zeros(size(x));
12 h = x(2)-x(1); 14 h = x(2)-x(1);
13 tot = m_order+s_order; 15 tot = m_order+s_order;
14 S = []; 16 S = [];
15 M = []; 17 M = [];
16 poss = find(tot*h/2 >= eta); 18 poss = find(tot*h/2 >= eta);
17 19
18 % Ensure that poss is not too long 20 % Ensure that poss is not too long
19 if length(poss) == (tot + 2) 21 if length(poss) == (tot + 2)
20 poss = poss(2:end-1); 22 poss = poss(2:end-1);
21 elseif length(poss) == (tot + 1) 23 elseif length(poss) == (tot + 1)
22 poss = poss(1:end-1); 24 poss = poss(1:end-1);
23 end 25 end
24 26
25 % Use first tot grid points 27 % Use first tot grid points
26 if length(poss)<tot && eta(end)>eta(1) 28 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
27 index=1:tot; 29 index=1:tot;
28 pol=(x(1:tot)-x(1))/(x(tot)-x(1)); 30 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
29 x_0=(x_0in-x(1))/(x(tot)-x(1)); 31 x_0=(x_0in-x(1))/(x(tot)-x(1));
30 norm=fnorm(1:tot)/h; 32 norm=fnorm(1:tot)/h;
31 33
32 % Use last tot grid points 34 % Use last tot grid points
33 elseif length(poss)<tot && eta(end)<eta(1) 35 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
34 index = length(x)-tot+1:length(x); 36 index = length(x)-tot+1:length(x);
35 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1)); 37 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
36 norm = fnorm(end-tot+1:end)/h; 38 norm = fnorm(end-tot+1:end)/h;
37 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1)); 39 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
38 40
41 % Interior, compensate for round-off errors.
42 elseif length(poss) < tot
43 if poss(end)<m
44 poss = [poss; poss(end)+1];
45 else
46 poss = [poss(1)-1; poss];
47 end
48 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
49 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
50 norm = fnorm(poss)/h;
51 index = poss;
52
39 % Interior 53 % Interior
40 else 54 else
41 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); 55 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
42 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); 56 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
43 norm = fnorm(poss)/h; 57 norm = fnorm(poss)/h;