In [1]:
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
# Imports
import scipy as sp
import scipy.linalg as la
In [41]:
# Make a chebychev domain for the unit disk
Nr = 10
Ntheta_twos = 10
## Domain in radius
Nx = 2 * (Nr + 1)
jx = sp.arange(Nx)
x = sp.cos(jx * sp.pi / (Nx-1))
jr = sp.arange(1, Nr+1)
r = x[jr]
plt.plot(x, sp.ones_like(x), 'o', color='white')
plt.plot(r, sp.ones_like(r), 'x', color='red', markersize=11)
## Domain in angle
Ntheta = 2 * Ntheta_twos
jtheta = sp.arange(Ntheta) # includes theta = 0 = 2*pi
theta = 2 * sp.pi * jtheta / Ntheta
## Domain in both
rr, tt = sp.meshgrid(r, theta, indexing='ij')
xx = rr * sp.cos(tt)
yy = rr * sp.sin(tt)
# Make a figure of the nodes
f0 = plt.figure(0)
f0.clear()
a0 = f0.add_axes([.1, .1, .8, .8], polar=True)
p0 = a0.plot(tt, rr, 'o', color='white', markersize=2)
# Plot a function
fpolar = (rr-1)**2 * (rr)**2 * sp.sin(3*tt)
plt.contourf(tt, rr, fpolar, polar=True)
Out[41]:
In [42]:
# Assemble the differentiation matrices
#N = Nx
#j = sp.arange(N+1)
#x = _r
D = sp.zeros((Nx, Nx))
N = Nx - 1
## fill in the middle values (eveywhere)
Dmid = (-1.0)**(jx[1:-1, None] + jx[None, 1:-1]) / (x[1:-1, None] - x[None, 1:-1])
D[1:-1,1:-1] = Dmid#[1:-1, 1:-1]
## Fill in the diagonals
Ddiag = 0.0 * x #- x / (2 * (1 - x**2))
Idiag = sp.diag_indices_from(D)
D[Idiag] = Ddiag
## Fill in the top
top = 2 * (-1)**(jx[1:-1]) / (1 - x[1:-1])
bot = -2 * (-1)**(N + jx[1:-1]) / (1 + x[1:-1])
left = -.5 * (-1)**(jx[1:-1]) / (1 - x[1:-1])
right = .5 * (-1)**(N + jx[1:-1]) / (1 + x[1:-1])
D[0, 1:-1] = top
D[-1, 1:-1] = bot
D[1:-1, 0] = left
D[1:-1, -1] = right
## Fill in the corners
#D[0, 0] = (2.0 * N**2 + 1) / 6.0
D[0, -1] = .5 * (-1)**N
D[-1, 0] = -.5 * (-1)**N
#D[-1, -1] = - (2.0 * N**2 + 1) / 6.0
## Fill diagonals by summing over non-diagonal entries
for i in range(N+1):
D[i, i] = -1.0 * D[i, :].sum()
## remove the edges
#D = D[1:-1, 1:-1]
# Store the matrix as Dr to make distinct from Dtheta
Dr = D[1:-1, 1:-1]
Dr2 = sp.dot(D, D)[1:-1, 1:-1]
## Plot the differentiation matrix
#plt.matshow(Dr)
#plt.title("The Chebyshev differentiation matrix")
# show the matrix
f1 = plt.figure(1)
a1 = f1.add_axes([.1, .1, .8, .8])
#Dr2 = sp.dot(Dr, Dr)
plotD = a1.matshow(Dr)
a1.set_title("Chebychev differentiation matrix")
f1.colorbar(plotD)
## compute a function and differentiate it
_r = x[1:-1]
fvals = 2 * (_r-1)**2 * (_r+1)**2
df = sp.dot(Dr, fvals) # take the second derivative
df_true = 2 * (2 * (_r-1)*(_r+1)**2 + 2 *(_r-1)**2 * (_r+1)) # second derivative
df2 = sp.dot(Dr2, fvals)
df2_true = 2.0 * (2 * (_r+1)**2
+ 4 * (_r-1) * (_r+1)
+ 4 * (_r-1) * (_r+1)
+ 2 * (_r-1)**2)
## Plot the test case for differentiation
f2 = plt.figure(2)
a2 = f2.add_axes([.1, .1, .8, .8])
a2.plot(_r, df_true-df_true, 'o', color='white')
a2.plot(_r, df2_true-df2_true, 'o', color='white')
a2.plot(_r, df-df_true, 'x', color='red', label='D')
a2.plot(_r, df2-df2_true, 'x', color='blue', label="D2")
a2.legend(loc='best')
#a2.set_ybound((.5, 1.1))
a2.set_xbound([-1.1, 1.1])
a2.set_title(r"Relative error in $\partial_{r}^{2}f(r)$ using $D_{r}^{2}$")
## Compute the eigenvalues
print(-1 * la.eigvals(D))
In [43]:
# Compute the differentiation matrix for theta (periodic domain)
## Compute the periodic Sinc function derivatives
v = sp.zeros(Ntheta)
h = 2 * sp.pi / Ntheta
v[1:] = (-1)**jtheta[1:] / (2.0 * sp.tan(jtheta[1:] * h / 2))
## Store these values in the differentiation array by rolling
Dtheta = sp.zeros((Ntheta, Ntheta))
for i in range(Ntheta):
Dtheta[i, :] = sp.roll(v, i)
_Dtheta2 = sp.dot(Dtheta, Dtheta)
## Compute the second order derivative
Dtheta2 = sp.zeros((Ntheta, Ntheta))
v = sp.zeros(Ntheta)
v[0] = -sp.pi**2 / (3 * h**2) - 1.0 / 6.0
v[1:] = -1 * (-1)**(jtheta[1:]) / (2 * sp.sin(jtheta[1:]*h/2)**2)
for i in range(Ntheta):
Dtheta2[i, :] = sp.roll(v, i)
print("error = {0}".format(abs(Dtheta2-_Dtheta2).mean() / abs(Dtheta2).mean()))
print(la.eigvalsh(Dtheta2))
## Show the differentiation arrays
plt.matshow(Dtheta)
plt.colorbar()
plt.title("Periodic differentiation array")
plt.matshow(sp.dot(Dtheta, Dtheta))
plt.colorbar()
plt.title("Periodic laplacian array")
## Show the error
fvals = -sp.sin(theta)
df2_true = sp.sin(theta)
df2 = sp.dot(Dtheta2, fvals)
f3 = plt.figure(3)
a3 = f3.add_axes([.1, .1, .8, .8])
a3.plot(theta, df2_true - df2_true, 'o', color='white')
a3.plot(theta, df2 - df2_true, 'x', color='red')
Out[43]:
In [44]:
# Compute the matricies needed for the polar coordinate problem
#Dr2 = sp.dot(Dr, Dr)
_D1 = Dr2[:Nr, :Nr][:, None, :, None]
_D2 = Dr2[:Nr, Nr:][:, None, :, None]
_E1 = Dr[:Nr, :Nr][:, None, :, None]
_E2 = Dr[:Nr, Nr:][:, None, :, None]
R = sp.diag(1.0 / r)[:, None, :, None]
R2 = sp.diag(1.0 / r**2)[:, None, :, None]
_Dtheta2 = Dtheta2[None, :, None, :]
Itheta = (sp.eye(Ntheta)[None, :, None, :])
Itheta_roll = sp.roll(sp.eye(Ntheta), Ntheta//2, 1)[None, :, None, :]
plt.matshow(Itheta.reshape(Ntheta, -1))
plt.matshow(Itheta_roll.reshape(Ntheta, -1))
## Define the laplacian in polar coordinates shape ~ (Nr, Ntheta, Nr, Ntheta)
L = sp.zeros((Nr, Ntheta, Nr, Ntheta))
L += (_D1 * Itheta + _D2 * Itheta_roll
+ R*_E1 * Itheta + R * _E2 * Itheta_roll
+ R2 * _Dtheta2)
L = L.reshape(Nr*Ntheta, Nr*Ntheta)
print(L.shape)
plt.matshow(L)
plt.colorbar()
#la.eigvals(L)
Out[44]:
In [45]:
(L<1e-10).sum() / float(L.size)
Out[45]:
In [46]:
sp.sqrt(-1 * la.eigvalsh(L))
Out[46]:
In [47]:
# Plot a function
f10 = plt.figure(10)
a10 = f10.add_axes([.1, .1, .8, .8], polar=True)
fpolar = (rr-1)**3 * (rr)**3 * sp.sin(2*tt)
cf = a10.contourf(tt, rr, fpolar, polar=True)
f10.colorbar(cf)
# Plot a function
f11 = plt.figure(11)
a11 = f11.add_axes([.1, .1, .8, .8], polar=True)
#compute laplacian
Lf = sp.dot(L, fpolar.flatten()).reshape([Nr, Ntheta])
cf = a11.contourf(tt, rr, Lf, polar=True)
f11.colorbar(cf)
Out[47]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: