Mercurial > repos > public > sbplib
comparison diracDiscr.m @ 645:a5307adc6177 feature/d1_staggered
Add diracDiscr.m. Noticed bug if source is on grid point.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 14 Nov 2017 14:28:35 -0800 |
| parents | |
| children | 0990765e3e4d |
comparison
equal
deleted
inserted
replaced
| 644:dc2918fb104d | 645:a5307adc6177 |
|---|---|
| 1 function ret = diracDiscr(x_0in , x , m_order, s_order, H) | |
| 2 | |
| 3 fnorm = diag(H); | |
| 4 eta = abs(x-x_0in); | |
| 5 h = x(2)-x(1); | |
| 6 tot = m_order+s_order; | |
| 7 S = []; | |
| 8 M = []; | |
| 9 poss = find(tot*h/2>eta); | |
| 10 | |
| 11 % Use first tot grid points | |
| 12 if length(poss)<tot && eta(end)>eta(1) | |
| 13 index=1:tot; | |
| 14 pol=(x(1:tot)-x(1))/(x(tot)-x(1)); | |
| 15 x_0=(x_0in-x(1))/(x(tot)-x(1)); | |
| 16 norm=fnorm(1:tot)/h; | |
| 17 | |
| 18 % Use last tot grid points | |
| 19 elseif length(poss)<tot && eta(end)<eta(1) | |
| 20 index = length(x)-tot:tot; | |
| 21 pol = (x(end-tot:end)-x(end-tot))/(x(end)-x(end-tot)); | |
| 22 norm = fnorm(end-tot:end)/h; | |
| 23 x_0 = (x_0in-x(end-tot))/(x(end)-x(end-tot)); | |
| 24 | |
| 25 % Interior | |
| 26 else | |
| 27 pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); | |
| 28 x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); | |
| 29 norm = fnorm(poss)/h; | |
| 30 index = poss; | |
| 31 end | |
| 32 | |
| 33 h_pol = pol(2)-pol(1); | |
| 34 b = zeros(m_order+s_order,1); | |
| 35 | |
| 36 for i = 1:m_order | |
| 37 b(i,1) = x_0^(i-1); | |
| 38 end | |
| 39 | |
| 40 for i = 1:(m_order+s_order) | |
| 41 for j = 1:m_order | |
| 42 M(j,i) = pol(i)^(j-1)*h_pol*norm(i); | |
| 43 end | |
| 44 end | |
| 45 | |
| 46 for i = 1:(m_order+s_order) | |
| 47 for j = 1:s_order | |
| 48 S(j,i) = (-1)^(i-1)*pol(i)^(j-1); | |
| 49 end | |
| 50 end | |
| 51 | |
| 52 | |
| 53 A = [M;S]; | |
| 54 | |
| 55 d = A\b; | |
| 56 ret = x*0; | |
| 57 ret(index) = d/h*h_pol; | |
| 58 end | |
| 59 | |
| 60 | |
| 61 | |
| 62 | |
| 63 | |
| 64 |