In [3]:
from numpy import *
from matplotlib.pyplot import *
from scipy import io
import numpy.polynomial.legendre as leg
%matplotlib inline
In [4]:
def vsf(th):
return .5 * np.exp(-th/2) / (1 - np.exp(-np.pi/2))
# Smallest difference between two angles
def angle_diff(th,thp):
return np.pi - abs(abs(th - thp) - pi)
In [5]:
# Check to see whether sum of scattering coefs. = 2
# nth = number of theta values
# kk = integer index num. of theta to diff other angles wrt.
def check_scat_sum(nth,kk):
th_un, ww_un = leg.leggauss(nth)
th_vals = pi * (th_un + 1)
weights = pi * ww_un
th = th_vals[kk]
thp = th_vals[th_vals!=th]
ww = weights[th_vals!=th]
return sum(vsf(angle_diff(th,thp)) * ww)
In [6]:
figure(figsize=[10,6])
for nth in range(10,160,10):
res = zeros(nth)
for kk in range(nth):
res[kk] = check_scat_sum(nth,kk)
plot(linspace(0,1,nth),res,label=r'$n_\theta = {}$'.format(nth))
legend(bbox_to_anchor=[1.01,1],loc='upper left')
xlabel(r'$k / n_\theta$')
ylabel(r"Sum of scattering coefficients")
title(r"$\sum_{l \neq k} \beta(|\theta_k - \theta_l|) w_l \
< \sum_{l=1}^{n_\theta} \beta(|\theta_k - \theta_l|) w_l \
\approx \int_0^{2\pi}\beta(|\theta_k-\theta'|) \,d\theta' = 2$"+'\n',fontsize=16)
tight_layout()
While the coefficients are not exactly 2, they do seem to be bounded above by 2 Adding in the l=k term to the sum reduces the error significantly, so all sums are much closer to 2 ~5e-2. However, 2 is no longer an upper bound then.
In [11]:
# Plot diagonal dominance quantity (D - Q) as a function of row number
# Add quantity mu to D to get D' = D + Q. If D' > Q <=> D' - Q = D + mu - Q > 0 (plot is positive),
# then diagonal dominance is achieved
def plot_dd(A,mu):
# Diagonal dominance quantity, D - Q (pre-shift)
ddq = zeros(nrows)
row_list = arange(nrows)
for k in row_list:
if k%100 == 0:
print('k={}'.format(k))
D = abs(A[k,k])
Q = abs(A[k,:k]).sum() + abs(A[k,k+1:]).sum()
ddq[k] = D - Q
# Post-shift DD quantity, D' - Q
ddq_p = ddq + mu
figure(figsize=[10,6])
plot(row_list,ddq,label="$D - Q$ (pre)")
plot(row_list,ddq_p,label="$D' - Q$ (post)")
print("Min (pre) : {:.3e}".format(min(ddq)))
print("Min (post): {:.3e}".format(min(ddq_p)))
title('Diagonal dominance plot')
legend(bbox_to_anchor=(1.01,1),loc='upper left')
In [12]:
# Grid parameters
nx = 20
ny = 20
nth = 24
nrows = nx * ny * nth
# Load matrix to test
dct = io.loadmat('../mat/ddom_{}x{}x{}_012.mat'.format(nx,ny,nth))
A = dct['A']
dx = 1/nx
dy = 1/ny
# IOPS
a_w = 1
b_w = 1
a_k = 5
b_k = 1
# Maximum number density of individuals
ind_up = 2
# Max/min abs. coef.
a_up = a_k * ind_up
a_dn = a_w
# Max/min scat. coef.
b_up = a_k * ind_up
b_dn = b_w
In [13]:
# Term to add to create diagonal dominance
# mu = a_up - b_dn * (2*pi - 1) - 1/dx - 4/dy # (technically correct)
# mu = a_up - b_dn * (2*pi - 1) - nx - 4*ny # (equivalent for normailized length scale)
mu = 1/dx + 4/dy + b_up * (2*pi-1) - a_dn
In [14]:
plot_dd(A,mu)
In [ ]:
# Try to solve DDOM w/ spsolve & gauss_seidel