Grid

Grid(ID, name, T, N, r, r′, r′′, h, r0, epn, epw, k, p, polynom)

Type with fields:

• .ID: grid identifer name (::Int)
• .name: grid identifer name (::String)
• .T: gridtypename (::Type)
• .N: number of grid points (::Int)
• .r: tabulated grid function (::Vector{T})
• .r′: tabulated first derivative of grid function (::Vector{T})
• .r′′: tabulated second derivative of grid function (::Vector{T})
• .h : grid step multiplyer (::T)
• .r0: grid scale factor (::T)
• .epn: number of endpoints used for trapezoidal endpoint correction (must be odd) (::Int)
• .epw: trapezoidal endpoint weights for n=1:epn (::Vector{Vector{T}})
• .k: finite-difference order (::Int)
• .p: only for truncated-exponential grid; truncation power (::Int)
• .polynom: only for polynomial grid: polynom (::Vector{T})

The object Grid is best created with the function castGrid.

castGrid(ID::Int, N::Int, T::Type; h=1, rmax=10, p=5, polynom=[0,1], epn=5, k=5, msg=false)
castGrid(name::String, N::Int, T::Type; h=1, rmax=10, p=5, polynom=[0,1], epn=5, k=5, msg=false)

Method to create a Grid object that covers the radial range [0, rmax] with N points.

ID = 1: exponential, ID = 2: truncated-exponential, ID = 3: linear (uniform) ID = 4: polynomial

Examples:

julia> grid = castGrid(1, 1000, Float64; h = 0.005, rmax = 10, msg=true);
Grid: exponential, Float64, rmax = 10.0, N = 1000, h = 0.005, r0 = 0.0681789

julia> grid = castGrid("exponential", 1000, Float64; h = 0.005, rmax = 10, msg=true);
Grid: exponential, Float64, rmax = 10.0, N = 1000, h = 0.005, r0 = 0.0681789

julia> grid = castGrid(2, 1000, Float64; h = 0.005, rmax = 10, p=5, msg=true);
Grid: truncated-exponential, Float64, rmax = 10.0, N = 1000, p = 5, h = 0.005, r0 = 0.111

julia> grid = castGrid(3, 1000, Float64; h = 0.1, rmax = 10, msg=true);
Grid: linear (uniform), Float64, rmax = 10.0, N = 1000, p = 1, h = 0.1, r0 = 0.1001

julia> grid = castGrid(4, 1000, Float64; h = 0.1, rmax = 10, polynom=[0,0,1], msg=true);
Grid: polynomial of degree 2, Float64, rmax = 10.0, N = 1000, polynom = [0.0, 0.0, 1.0], h = 0.1, r0 = 0.001002

julia> r = grid.r[1:4]; println("r = ", r)
r = [0.0, 1.002003004005006e-5, 4.008012016020024e-5, 9.018027036045053e-5]

julia> r = grid.r[997:1000]; println("r = ", r)
r = [9.9400301202103, 9.96000004008012, 9.979990000010021, 10.0] # note the end of the range (r = rmax) 

julia> r′= grid.r′[1:4]; println("r′ = ", r′)
r′ = [0.0, 2.0040060080100123e-5, 4.008012016020025e-5, 6.012018024030035e-5]

julia> r′′= grid.r′′[1:4]; println("r′′ = ", r′′)
r′′ = [2.004006008010012e-5, 2.004006008010012e-5, 2.004006008010012e-5, 2.004006008010012e-5]

gridfunction(ID::Int, n::Int, T::Type; h=1, p=5, polynom=[0,1], deriv=0)

CamiDiff offers three internal grid functions:

• ID = 1: exponential grid function,

\[ g(t) = \rm{exp}(t) - 1\]

• ID = 2: truncated-exponential grid function of degree p (linear grid for p = 1),

\[ g(t) = t + \frac{1}{2}t^2 + ⋯ + \frac{1}{p!}t^p\]

• ID = 3: linear grid function,

\[ g(t) = t\]

• ID = 4: polynomial grid function of degree p = length(c)-1 defined by its CamiMath.polynom vector $c = [c_0, c_1,c_2,⋯\ c_p]$,

\[ g(t) = c_0 + c_1 t + c_2 t^2 + ⋯ + c_p t^p,\]

with $c_0 ≡ 0$ because, by definition, all gridfunctions run through the origin, $g(0) = 0$.

The actual Grid is given by

\[ r[n] = r_0 * g(t[n]),\]

where $t[n] = (n-1) * h$ is the ticks function for the unit-based indexing of Julia.

NB. All gridfunctions satisfy the properties $t[1] = 0$ and $r[1] = 0$.

Examples:

julia> h = 0.1; r0=1.0; N=4; T=Float64;

julia> r = [r0*gridfunction(1, n, T; h) for n=1:N]; println("r = ", r)
r = [0.0, 0.10517091807564771, 0.22140275816016985, 0.3498588075760032]

julia> r′= r0 .* [r0*gridfunction(1, n, T; h, deriv=1) for n=1:N]; println("r′= ", r′)
r′ = [0.1, 0.11051709180756478, 0.122140275816017, 0.13498588075760032]

julia> r′′= [r0*gridfunction(1, n, T; h, deriv=2) for n=1:N]; println("r′′= ", r′′)
r′′ = [0.010000000000000002, 0.011051709180756479, 0.012214027581601701, 0.013498588075760034]

julia> r = [r0*gridfunction(4, n, T; h, polynom=[0,0,1]) for n=1:N]; println("r = ", r)
r = [0.0, 0.010000000000000002, 0.04000000000000001, 0.09000000000000002]

julia> r′= [r0*gridfunction(4, n, T; h, polynom=[0,0,1], deriv=1) for n=1:N]; println("r′= ", r′)
r′= [0.0, 0.020000000000000004, 0.04000000000000001, 0.06000000000000001]

julia> r′′= [r0*gridfunction(4, n, T; h, polynom=[0,0,1], deriv=2) for n=1:N]; println("r′′= ", r′′)
r′′= [0.020000000000000004, 0.020000000000000004, 0.020000000000000004, 0.020000000000000004]

gridPos(rval::T, grid::Grid{T}) where T<:Number

The approximate grid position defined as the largest integer n satisfying the condition grid.r[n] < rval on the Grid.

Example:

Consider the exponential grid of 4 points defined by

julia> grid = castGrid("exponential", 4, Float64; h = 0.1, rmax = 2.0);

julia> println(grid.r)
[0.0, 0.6012192107114549, 1.265669197778149, 2.0]

The approximate grid position of the point $r = 1.0$ is $2$ ($n = 2$); i.e., $r = 1.0$ is larger than $0.6012192107114549$ ($n = 2$) but smaller than $1.265669197778149$ ($n = 3$).

julia> r = 1.0;

julia> n = gridPos(r, grid)
2

fracPos(n::Int, rval::T, grid::Grid{T}; ϵ = 1e-8, k = 7) where T<:Real

Fractional grid offset with respect to Grid position n.

Example:

Consider the exponential grid of 4 points defined by

julia> grid = castGrid("exponential", 4, Float64; h = 0.1, rmax = 2.0);

julia> println(grid.r)
[0.0, 0.6012192107114549, 1.265669197778149, 2.0]

The point $r = 1.0$ is located just above $n = 2$, with fractional position $Δn = 0.6120806373655796$.

julia> r = 1.0;

julia> n = gridPos(r, grid)
2

julia> Δn = fracPos(n, r, grid);

julia> t = (n+Δn-1)*grid.h;

julia> grid.r0 * (exp(t)-1) ≈ r

regularize!(f::Vector{T} [; k=3]) where T<:Real

Regularization of a function with a non-analytic point in the origin.

Example:

julia> r = [0,1,2,3,4,5];

julia> f = r ./ r; println(f)
[NaN, 1.0, 1.0, 1.0, 1.0, 1.0]

julia> regularize!(f); ; println(f)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

gridtypename(ID::Int)

Name corresponding to the Grid ID.

Example:

julia> gridtypename(2)
"truncated-exponential"

gridtypeID(name::String)

ID corresponding to the gridtypename.

Example:

julia> gridtypeID("truncated-exponential")
2

grid_interpolation(f::Vector{T}, grid::Grid{T}, rv::T, notation=fwd; k=5) where T<:Real

$k^{th}$-order lagrangian interpolation of the function $f(r)$ at position r

• f[1:N] : the function f(r) tabulated in forward order on a Grid of N points
• fwd using fwd-difference notation
• bwd using bwd-difference notation

Example:

julia> grid = castGrid(1, 1000, Float64; h = 0.01, rmax=25, msg=false);

julia> r = grid.r;

julia> f = [exp(-r[n]) for n=1:grid.N];

julia> rv = 1.0;

julia> grid_interpolation(f, grid, rv, fwd) ≈ exp(-rv)
true

grid_differentiation(f::Vector{T}, grid::Grid{T}; k=5) where T<:real

$k^{th}$-order lagrangian derivative of the function $f(r)$ over the range $0 ≤ r ≤ rmax$

• f[1:N] : the function f(r) tabulated in forward order on a Grid of N points

Example:

julia> grid = castGrid(3, 1001, Float64; h=2π/1000.0, r0=1.0, msg=false);

julia> f = [sin(grid.r[i]) for i=1:grid.N]

julia> f′ = [cos(grid.r[i]) for i=1:grid.N]

julia> grid_differentiation(f, grid) ≈ f′
true

grid_differentiation(f::Vector{T}, grid::Grid{T}, rv::T, notation=fwd; k=5) where T<:Real

$k^{th}$-order lagrangian derivative of the function $f(r)$ at position rv

• f[1:N] : the function f(r) tabulated in forward order on a Grid of N points
• fwd using fwd-difference notation
• bwd using bwd-difference notation

Example:

julia> grid = castGrid(3, 1001, Float64; h=2π/1000.0, r0=1.0, msg=false);

julia> f = [sin(grid.r[i]) for i=1:grid.N];

julia> grid_differentiation(f, grid, float(π), fwd) ≈ cos(π)
true

grid_differentiation(f::Vector{T}, grid::Grid{T}, n::Int, notation=fwd; k=5) where T<:Real

$k^{th}$-order lagrangian derivative of the function $f(r)$ at grid position n

• f[1:N] : the function f(r) tabulated in forward order on a Grid of N points
• fwd using fwd-difference notation
• bwd using bwd-difference notation

Example:

julia> grid = castGrid(3, 1001, Float64; h=2π/1000.0, r0=1.0, msg=false);

julia> f = [sin(grid.r[i]) for i=1:grid.N];

julia> f′ = [cos(grid.r[i]) for i=1:grid.N];

julia> grid_differentiation(f, grid, 500, fwd) ≈ f′[500]
true

grid_differentiation(f::Vector{T}, grid::Grid{T}, n1::Int, n2::Int; k=5) where T<:Real
grid_differentiation(f::Vector{T}, grid::Grid{T}, itr::UnitRange; k=5) where T<:Real

$k^{th}$-order lagrangian derivative of the function $f(r)$ over the grid range n1 ≤ n ≤ n2

• f[1:N] : the function f(r) tabulated in forward order on a Grid of N points
• n1=itr.start, n2=itr.stop.

grid_integration(f::Vector{T}, grid::Grid{T}) where T<:Real
grid_integration(f::Vector{T}, grid::Grid{T}, n1::Int, n2::Int) where T<:Real
grid_integration(f::Vector{T}, grid::Grid{T}, itr::UnitRange) where T<:Real

Integral of the analytic function $f(r)$ tabulated on a generally nonlinear Grid and evaluated with the generalized trapezoidal rule optimized with endpoint correction by the weightsvector grid.epw,

\[ ∫_{0}^{r_n} f(r) dr = ∫_{0}^{n} f(x) r^{\prime}(x) dx.\]

Here, the latter integral corresponds to the optimized trapezoidal rule for a uniform grid (see trapezoidal_integration). The rule is exact for polynomial functions of degree $d=0,\ 1,⋯\ k-1$, where $k=$ grid.epn. For $k=1$ the rule reduces to the ordinary trapezoidal rule (weights = [1/2]).

Examples:

julia> ftest(r) = sqrt(2.0/π) * exp(-r^2/2.0);

julia> grid1 = castGrid(1, 1000, Float64; h = 0.005, r0 = 0.1, msg=true);
Grid created: exponential, Float64, rmax = 14.7413, N = 1000, h = 0.005, r0 = 0.1

julia> grid2 = castGrid(2, 1000, Float64; h = 0.005, r0 = 0.1, p=5, msg=true);
Grid created: truncated-exponential, Float64, rmax = 9.04167, N = 1000, p = 5, h = 0.005, r0 = 0.1

julia> grid3 = castGrid(3, 1000, Float64; h = 0.1, r0 = 0.1, msg=true);
Grid created: linear (uniform), Float64, rmax = 10.0, N = 1000, p = 1, h = 0.1, r0 = 0.1

julia> grid4 = castGrid(4, 1000, Float64; h = 0.1, r0 = 0.001, polynom=[0,0,1], msg=true);
Grid created: polynomial, Float64, rmax = 10.0, N = 1000, polynom = [0.0, 0.0, 1.0], h = 0.1, r0 = 0.001

julia> r1 = grid1.r;
julia> r2 = grid2.r;
julia> r3 = grid3.r;
julia> r4 = grid4.r;
julia> f1 = [ftest(r1[n]) for n=1:grid1.N];
julia> f2 = [ftest(r2[n]) for n=1:grid2.N];
julia> f3 = [ftest(r3[n]) for n=1:grid3.N];
julia> f4 = [ftest(r4[n]) for n=1:grid4.N];
julia> o1 = grid_integration(f1, grid1);
julia> o2 = grid_integration(f2, grid2);
julia> o3 = grid_integration(f3, grid3, 1:900);
julia> o4 = grid_integration(f4, grid4, 1:900);

julia> println("integral on " * grid1.name * " grid = ", o1)
integral on exponential grid = 1.0

julia> println("integral on " * grid2.name * " grid = ", o2)
integral on truncated-exponential grid: 1.0

julia> println("integral on " * grid3.name * " grid = ", o3)
integral on linear (uniform) grid = 1.000000000000003

julia> println("integral on " * grid3.name * " grid = ", o4)
integral on polynomial grid = 1.0000000000000013