Prototype phi functions (scalar and matrix versions)

Phi functions.md

Phi functions

The phi functions can be defined using the recurrence relation

$$ \varphi_{k+1}(z) = \frac{\varphi_k(z) - \varphi_k(0)}{z},\quad \varphi_0(z) = e^z $$

The zero point values are $\varphi_k(0) = 1/k!$.

The recurrence relation naturally lends itself to a dynamic programming style implementation, where it is more efficient to calculate all $\varphi_i(z)$ for $0 \le i \le k$.

Scalar version

The scalar version comes in handy when $A$ is known to be diagonalizable (in particular, when it is symmetric/Hermitian), in which case computing $\varphi_k(A)$ using diagonalization is better.

"""
    phi(z,k) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the value of the phi functions of z for all orders up to k.
"""
function phi(z::T, k) where {T <: Number}
    # Compute zero points
    zero_points = Vector{T}(k+1)
    zero_points[1] = one(T)
    for i = 1:k
        zero_points[i+1] = zero_points[i] / i
    end
    
    if iszero(z)
        return zero_points
    else
        # Compute phi(z) using the recurrence relation
        out = Vector{T}(k+1)
        out[1] = exp(z) # phi_0(z) = exp(z)
        for i = 1:k
            out[i+1] = (out[i] - zero_points[i]) / z
        end
        return out
    end
end;

phi(1e-3, 4)

5-element Array{Float64,1}:
 1.001    
 1.0005   
 0.500167 
 0.166708 
 0.0417179

phi(1e-3 + 0.0im, 4)

5-element Array{Complex{Float64},1}:
     1.001+0.0im
    1.0005+0.0im
  0.500167+0.0im
  0.166708+0.0im
 0.0417179+0.0im

phi(big(1e-3), 4)

5-element Array{BigFloat,1}:
 1.001000500166708341688893262798302888729404554195754965197362071336734197553476    
 1.000500166708341668066169275410855560811798459447420824073778226865822743593049    
 5.001667083416680557574642405075253788758590184408141904359913822670546237956644e-01
 1.667083416680557539939260203325124970459099317152964308056605779086203103308218e-01
 4.167500138908732639181842659364040804292368276439509042852030958052561993062031e-02

Matrix version

"""
    phim(A,k) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the value of the matrix phi functions of A for all orders up to k.
"""
function phim(A::AbstractMatrix{T}, k) where {T <: Number}
    @assert size(A,1) == size(A,2) "Dimension mismatch"
    N = size(A,1)
    # Compute zero points
    zero_points = Vector{Matrix{T}}(k+1)
    zero_points[1] = eye(T, N)
    for i = 1:k
        zero_points[i+1] = zero_points[i] ./ i
    end
    
    if iszero(A)
        return zero_points
    else
        # Compute phi(A) using the recurrence relation
        out = Vector{Matrix{T}}(k+1)
        out[1] = expm(A) # phi_0(A) = exp(A)
        for i = 1:k
            out[i+1] = A \ (out[i] - zero_points[i])
        end
        return out
    end
end;

A = 1e-3 * rand(2,2)
phim(A, 4)

5-element Array{Array{Float64,2},1}:
 [1.00006 0.000588689; 0.000689637 1.00047]    
 [1.00003 0.000294318; 0.000344788 1.00024]    
 [0.50001 9.81015e-5; 0.000114924 0.500079]    
 [0.166669 2.52552e-5; 2.88701e-5 0.166686]    
 [0.0420821 -0.0016183; -0.000304539 0.0430815]

A = 1e-3 * (rand(2,2) + rand(2,2)*1im)
phim(A, 4)

5-element Array{Array{Complex{Float64},2},1}:
 Complex{Float64}[1.00039+0.000450626im 0.000520239+8.26741e-5im; 0.000320235+0.000313714im 1.00089+0.00051044im]   
 Complex{Float64}[1.00019+0.000225268im 0.000260071+4.12866e-5im; 0.000160109+0.000156798im 1.00044+0.000255129im]  
 Complex{Float64}[0.500065+7.50822e-5im 8.66823e-5+1.37539e-5im; 5.33681e-5+5.22559e-5im 0.500148+8.50276e-5im]     
 Complex{Float64}[0.166684+1.84696e-5im 2.22459e-5+3.20231e-6im; 1.25508e-5+1.28725e-5im 0.166703+2.12101e-5im]     
 Complex{Float64}[0.0428557-0.00442375im 0.000211594-0.00200878im; -0.00174884+0.00197766im 0.0410412+0.000971505im]

expm is not defined for bigfloats. Here's an ad-hoc "fix":

import Base.expm
expm(A::AbstractMatrix{BigFloat}) = big.(expm(Float64.(A)))
expm(A::AbstractMatrix{Complex{BigFloat}}) = big.(expm(Complex128.(A)));

A = big.(1e-3 * rand(2,2))
phim(A, 4)

5-element Array{Array{BigFloat,2},1}:
 BigFloat[1.000648620898963336856013484066352248191833496093750000000000000000000000000000 5.590826613970026207414698582454093411797657608985900878906250000000000000000000e-04; 3.812815994498952912837574924509453921928070485591888427734375000000000000000000e-04 1.000145260150193449177891125145833939313888549804687500000000000000000000000000]        
 BigFloat[1.000324257650573401023248814035382104095563377955766960688538597840504584291828 2.795043602293500341094181357164102806190025069698612443347187501150603706442567e-04; 1.906155872486168273175545883699931231257441055190328240108753739666022743828566e-04 1.000072610562138804623510220201247452656182503593195562498606688804247676307962]        
 BigFloat[5.001080769881851237921546670051368454383762623535359481625508441782698991585526e-01 9.316741313655890454408016671208263836838843828753663400030418260122785842195713e-05; 6.353445173850038429334428629540216804734181530434842282811616006123778802656736e-05 5.000241924789780013552074709453902874036213746194673184309782522546594590672645e-01]
 BigFloat[1.666943865694994700727724839726121342955424713490380447491979504465781127296786e-01 -1.997939609939074958560523547000702648224290086744619129920824974422941408992846e-05; 1.489629969582110076009763044215445411731242312229191150594263941068256997643818e-05 1.66732669024940836022963906795015072700215229985993657052777716672466269417994e-01]
 BigFloat[3.617775796709503704785806302955147447781247290265358660654993069352994590764836e-02 3.34611820498423694297245008676604670842389745637363934520837712569930951092457e-01; 7.632869347691420280285586767218748841171527948531562825578127144034041930028434e-03 -4.239152146977518381974488318427355946195113788924899703904484598515053361477154e-01]