Mercurial > repos > public > sbplib

comparison diracDiscr.m @ 1237:6e4cc4b66de0 feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove overloading of norm() and fnorm
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 20 Nov 2019 00:13:56 +0100
parents 3722c2579818
children dea852e85b77
comparison
equal deleted inserted replaced
1236:3722c2579818 1237:6e4cc4b66de0
39 if(x_s < x(1) || x_s > x(end) ) 39 if(x_s < x(1) || x_s > x(end) )
40 40
41 ret = zeros(size(x)); 41 ret = zeros(size(x));
42 42
43 else 43 else
44 fnorm = diag(H);
45 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for 44 tot_order = m_order+s_order; %This is equiv. to the number of equations solved for
46 S = []; 45 S = [];
47 M = []; 46 M = [];
48 47
49 % Get interior grid spacing 48 % Get interior grid spacing
52 51
53 index = sourceIndecies(x_s, x, tot_order, h) 52 index = sourceIndecies(x_s, x, tot_order, h)
54 53
55 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1))); 54 polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1)));
56 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1))); 55 x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1)));
57 norm = fnorm(index)/h; 56
57 quadrature = diag(H);
58 quadrature_weights = quadrature(index)/h;
58 59
59 h_polynomial = polynomial(2)-polynomial(1); 60 h_polynomial = polynomial(2)-polynomial(1);
60 b = zeros(tot_order,1); 61 b = zeros(tot_order,1);
61 62
62 for i = 1:m_order 63 for i = 1:m_order
63 b(i,1) = x_0^(i-1); 64 b(i,1) = x_0^(i-1);
64 end 65 end
65 66
66 for i = 1:tot_order 67 for i = 1:tot_order
67 for j = 1:m_order 68 for j = 1:m_order
68 M(j,i) = polynomial(i)^(j-1)*h_polynomial*norm(i); 69 M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i);
69 end 70 end
70 end 71 end
71 72
72 for i = 1:tot_order 73 for i = 1:tot_order
73 for j = 1:s_order 74 for j = 1:s_order