Mercurial > repos > public > sbplib

comparison diracPrimDiscr1D.m @ 1130:99fd66ffe714 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add derivative of delta functions and corresponding tests, tested for 1D.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 May 2019 18:44:01 -0700
parents
children
comparison
equal deleted inserted replaced
1129:b29892853daf 1130:99fd66ffe714
1 % Generates discretized derivative of delta function in 1D
2 function ret = diracPrimDiscr1D(x_0in, x, m_order, s_order, H)
3
4 % diracPrim satisfies one more moment condition than dirac
5 m_order = m_order + 1;
6
7 m = length(x);
8
9 % Return zeros if x0 is outside grid
10 if(x_0in < x(1) || x_0in > x(end) )
11
12 ret = zeros(size(x));
13
14 else
15
16 fnorm = diag(H);
17 eta = abs(x-x_0in);
18 tot = m_order+s_order;
19 S = [];
20 M = [];
21
22 % Get interior grid spacing
23 middle = floor(m/2);
24 h = x(middle+1) - x(middle);
25
26 poss = find(tot*h/2 >= eta);
27
28 % Ensure that poss is not too long
29 if length(poss) == (tot + 2)
30 poss = poss(2:end-1);
31 elseif length(poss) == (tot + 1)
32 poss = poss(1:end-1);
33 end
34
35 % Use first tot grid points
36 if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
37 index=1:tot;
38 pol=(x(1:tot)-x(1))/(x(tot)-x(1));
39 x_0=(x_0in-x(1))/(x(tot)-x(1));
40 norm=fnorm(1:tot)/h;
41
42 % Use last tot grid points
43 elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
44 index = length(x)-tot+1:length(x);
45 pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
46 norm = fnorm(end-tot+1:end)/h;
47 x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
48
49 % Interior, compensate for round-off errors.
50 elseif length(poss) < tot
51 if poss(end)<m
52 poss = [poss; poss(end)+1];
53 else
54 poss = [poss(1)-1; poss];
55 end
56 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
57 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
58 norm = fnorm(poss)/h;
59 index = poss;
60
61 % Interior
62 else
63 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
64 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
65 norm = fnorm(poss)/h;
66 index = poss;
67 end
68
69 h_pol = pol(2)-pol(1);
70 b = zeros(m_order+s_order,1);
71
72 b(1) = 0;
73 for i = 2:m_order
74 b(i) = -(i-1)*x_0^(i-2);
75 end
76
77 for i = 1:(m_order+s_order)
78 for j = 1:m_order
79 M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
80 end
81 end
82
83 for i = 1:(m_order+s_order)
84 for j = 1:s_order
85 S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
86 end
87 end
88
89 A = [M;S];
90
91 d = A\b;
92 ret = x*0;
93 ret(index) = d*(h_pol/h)^2;
94 end
95
96 end