Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 649:1bdbe026abbc feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make diracDiscr return zeros if x0 is outside grid.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 14 Nov 2017 15:40:06 -0800
parents 9e5dd0d3cf60
children 4ee7d15bd8e6
comparison
equal deleted inserted replaced
648:9e5dd0d3cf60 649:1bdbe026abbc
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 2
3 fnorm = diag(H); 3 % Return zeros if x0 is outside grid
4 eta = abs(x-x_0in); 4 if(x_0in < x(1) || x_0in > x(end) )
5 h = x(2)-x(1);
6 tot = m_order+s_order;
7 S = [];
8 M = [];
9 poss = find(tot*h/2 >= eta);
10 5
11 % Ensure that poss is not too long 6 ret = zeros(size(x));
12 if length(poss) == (tot + 2) 7
13 poss = poss(2:end-1); 8 else
14 elseif length(poss) == (tot + 1) 9
15 poss = poss(1:end-1); 10 fnorm = diag(H);
11 eta = abs(x-x_0in);
12 h = x(2)-x(1);
13 tot = m_order+s_order;
14 S = [];
15 M = [];
16 poss = find(tot*h/2 >= eta);
17
18 % Ensure that poss is not too long
19 if length(poss) == (tot + 2)
20 poss = poss(2:end-1);
21 elseif length(poss) == (tot + 1)
22 poss = poss(1:end-1);
23 end
24
25 % Use first tot grid points
26 if length(poss)<tot && eta(end)>eta(1)
27 index=1:tot;
28 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
29 x_0=(x_0in-x(1))/(x(tot)-x(1));
30 norm=fnorm(1:tot)/h;
31
32 % Use last tot grid points
33 elseif length(poss)<tot && eta(end)<eta(1)
34 index = length(x)-tot+1:length(x);
35 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
36 norm = fnorm(end-tot+1:end)/h;
37 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
38
39 % Interior
40 else
41 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)));
43 norm = fnorm(poss)/h;
44 index = poss;
45 end
46
47 h_pol = pol(2)-pol(1);
48 b = zeros(m_order+s_order,1);
49
50 for i = 1:m_order
51 b(i,1) = x_0^(i-1);
52 end
53
54 for i = 1:(m_order+s_order)
55 for j = 1:m_order
56 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
57 end
58 end
59
60 for i = 1:(m_order+s_order)
61 for j = 1:s_order
62 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
63 end
64 end
65
66 A = [M;S];
67
68 d = A\b;
69 ret = x*0;
70 ret(index) = d/h*h_pol;
16 end 71 end
17 72
18 % Use first tot grid points
19 if length(poss)<tot && eta(end)>eta(1)
20 index=1:tot;
21 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
22 x_0=(x_0in-x(1))/(x(tot)-x(1));
23 norm=fnorm(1:tot)/h;
24
25 % Use last tot grid points
26 elseif length(poss)<tot && eta(end)<eta(1)
27 index = length(x)-tot+1:length(x);
28 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
29 norm = fnorm(end-tot+1:end)/h;
30 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
31
32 % Interior
33 else
34 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
35 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
36 norm = fnorm(poss)/h;
37 index = poss;
38 end
39
40 h_pol = pol(2)-pol(1);
41 b = zeros(m_order+s_order,1);
42
43 for i = 1:m_order
44 b(i,1) = x_0^(i-1);
45 end
46
47 for i = 1:(m_order+s_order)
48 for j = 1:m_order
49 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
50 end
51 end
52
53 for i = 1:(m_order+s_order)
54 for j = 1:s_order
55 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
56 end
57 end
58
59 A = [M;S];
60
61 d = A\b;
62 ret = x*0;
63 ret(index) = d/h*h_pol;
64 end 73 end
65 74
66 75
67 76
68 77