In [9]:
L = 10.0
N = 5
Nbasis = 2*N + 1
grid_x = Array{Float64}(Nbasis)
for α = 1:Nbasis
grid_x[α] = L/Nbasis*(α-N-1)
end
Δ = L/Nbasis
grid_x
Out[9]:
In [21]:
function eval_LF1d_p( α::Int64, L::Float64, grid_x::Array{Float64}, x)
Nbasis = size(grid_x)[1]
N = (Nbasis-1)/2
f = 0.0
for l = 1:Nbasis
k = -N + (l-1)
f = f + sqrt(1/L/Nbasis)*cos(2*pi*k*(x-grid_x[α])/L)
end
return f
end
Out[21]:
In [24]:
f1 = eval_LF1d_p( 1, L, grid_x, grid_x[2] )
f2 = eval_LF1d_p( 2, L, grid_x, grid_x[2] )
println(f1, " ", f2)
In [25]:
import PyPlot
const plt = PyPlot
Out[25]:
In [35]:
plt.clf()
NptsPlot = 200
x = Array{Float64}(linspace(-L,L,NptsPlot));
# Plot for basis function 1, 4, and 7
for α = [1,4,7]
y = Array{Float64}(NptsPlot)
for i = 1:size(x)[1]
y[i] = eval_LF1d_p(α, L, grid_x, x[i])
end
plt.grid()
plt.plot(x,y)
end
In [39]:
# Transformation matrix
T = Array{Complex128}(Nbasis,Nbasis)
for α = 1:Nbasis
for l = 1:Nbasis
k = -N + (l-1)
x = grid_x[α]
T[α,l] = sqrt(1/Nbasis)*exp(2*pi*im*k*x/L)
end
end
T
Out[39]:
In [42]:
# Test unitary property of transformation matrix
T' - inv(T)
Out[42]:
In [45]:
X_dvr = diagm(grid_x)
Out[45]:
In [49]:
X = T' * X_dvr * T
println(size(X))
println(real(X))
In [61]:
real(X)
imag(T'*X_dvr*T)
Out[61]:
In [66]:
g1 = T*(T'*grid_x)
Out[66]:
In [69]:
abs.(g1 - grid_x)
Out[69]:
In [70]:
T*grid_x
Out[70]:
In [73]:
sort( imag(T'*grid_x) )
Out[73]:
In [75]:
xx = sort( imag(T*grid_x) );
xx[2] - xx[1]
Out[75]:
In [76]:
grid_x[2] - grid_x[1]
Out[76]:
In [78]:
abs.(T'*T)
Out[78]:
In [112]:
L = 10.0
N = 10
Nbasis = 2*N + 1
grid_x = Array{Float64}(Nbasis)
for α = 1:Nbasis
grid_x[α] = L/Nbasis*(α-N-1)
end
Δ = L/Nbasis;
In [96]:
function myfunc(L::Float64, x::Float64)
ω = 2*pi/L
f = cos(ω*x)*sin(2*ω*x)
return f
end
Out[96]:
In [97]:
NptsPlot = 200
x = Array{Float64}(linspace(-L,L,NptsPlot));
y = Array{Float64}(NptsPlot)
for i = 1:NptsPlot
y[i] = myfunc(L,x[i])
end
plt.clf()
plt.grid()
plt.plot(x,y)
Out[97]:
In [113]:
ex_coefs = Array{Float64}(Nbasis)
for i = 1:Nbasis
ex_coefs[i] = myfunc(L, grid_x[i])
end
In [114]:
function eval_from_ex_coefs( ex_coef::Array{Float64}, L, grid_x, x )
f = 0.0
for i = 1:Nbasis
f = f + ex_coef[i]*eval_LF1d_p( i, L, grid_x, x )
end
return f
end
Out[114]:
In [116]:
eval_from_ex_coefs( ex_coefs, L, grid_x, grid_x[1] ), myfunc(L, grid_x[1])
Out[116]:
In [101]:
NptsPlot = 200
x = Array{Float64}(linspace(-L,L,NptsPlot));
y1 = Array{Float64}(NptsPlot)
y2 = Array{Float64}(NptsPlot)
for i = 1:NptsPlot
y1[i] = myfunc(L,x[i])
y2[i] = eval_from_ex_coefs( ex_coefs, L, grid_x, x[i])
end
plt.clf()
plt.grid()
plt.plot(x,y1)
plt.plot(x,y2)
In [ ]: