Examples
Linear-recursive Sequence
Consider the recurrence
\[x_1=x_2=x_3=\sqrt{2}, x_{n+3}-\frac{10}{3}x_{n+2}+3x_{n+1}-\frac{2}{3}x_n=0.\]
First, we shall recognise that the linear recursive sequence is basically a 1D stencil recurrence. The recurrence needs to be converted to an assignment first, that is, $x_{n}=\frac{10}{3}x_{n-1}-3x_{n-2}+\frac{2}{3}x_{n-3}$. Such a recurrence can be defined by the stencil with the coefficients
stencil = (CartesianIndex(-3), CartesianIndex(-2), CartesianIndex(-1))
coefs = (n -> 2//3, n -> -3, n -> 10//3)To define a linear inhomogeneous recurrence, the coefficient associated with the zero cartesian index is used. For example, the recurrence
\[x_{n}=\frac{10}{3}x_{n-1}-3x_{n-2}+\frac{2}{3}x_{n-3}+\frac{1}{n}\]
should be defined by
stencil = (CartesianIndex(-3), CartesianIndex(-2), CartesianIndex(-1), CartesianIndex(0))
coefs = (n -> 2//3, n -> -3, n -> 10//3, n -> 1//n)We then need to define the initial values
f_init(T) = [sqrt(T(2)), sqrt(T(2)), sqrt(T(2)), T(0)]where the place for the next step should be reserved. This is a function that can generate values based on type T.
Now we only need to further provide the size of the recurrence and we are ready to go.
P = StencilRecurrencePlan{Real}(stencil, coefs, f_init, (100,)) # 'Real' specifies the domain of the entries, as opposed to 'Complex', etc.
stable_recurrence(P) # defaults to stable_recurrence(P, Float64)100-element Vector{Array{Float64, 0}}:
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
⋮
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)
fill(1.4142135623730951)A definite integral
Consider the integral
\[I(m,n) = \int_0^\pi\left(\frac{\cos x}{1+\sin x}\right)^me^{inx}\mathrm{d}x\]
which has the recurrence
\[I(m,n)=\begin{cases} 2i/n,& m=0, n\text{ odd}\\ \displaystyle(-1)^{m/2}\left(\pi-4\sum_{k=1}^{m/2}\frac{(-1)^k}{2k-1}\right),& m\text{ even}, n=0\\ 0,& \!\!\!\begin{array}{l} m=0\text{ or }n=0,\\ \text{excluding above cases} \end{array}\\ \displaystyle\frac{m+n-1}{ni}I(m,n-1)+\frac{m}{n}I(m-1,n-1)+\frac{2}{n}\begin{cases} i,\\ -1, \end{cases}&\!\!\!\begin{array}{l} m+n\text{ odd}\\ m+n\text{ even}. \end{array} \end{cases}\]
Our goal is to obtain the numerical integrals for a range of $m$ and $n$. We may use a matrix to store the results, with m+1 and n+1 being the row and column indices.
The problem is again a stencil recurrence, where
stencil = (CartesianIndex(-1,-1), CartesianIndex(0,-1), CartesianIndex(0,0))and
coef1(m,n) = (m-1)//(n-1)
coef2(m,n) = (m+n-3)//(n-1)//im
coef3(m,n) = 2//(n-1)*ifelse(isodd(m+n), im, -1)
coef = (coef1, coef2, coef3)Note that Julia indices starts at 1, hence the m-1, n-1, etc.. Recall that the CartesianIndex(0,0) and coef3 corresponds to the inhomogeneous term $\displaystyle\frac{2}{n}\begin{cases}i, & m+n\text{ odd}\\ -1, & m+n\text{ even}\end{cases}$.
We need to initialise the first column as well.
function f_init(T, m)
A = zeros(Complex{T},m,2)
A[1,1] = π
for mm in 3:2:m
A[mm,1] = -A[mm-2,1] + 4//(mm-2)
end
A
endHere, m specifies the size of the first dimension. For an $n$-dimensional stencil recurrence, each slice has $n-1$ dimensions, so the initialization function has to take $n-1$ size arguments.
PseudostableRecurrences.jl doesn't support initialising the first row in the backend. However, if a stencil exceeds the boundary, the external entries will be ignored. This allows one to modify the coefficients to achieve the initialization of the first row. In this example, the recurrence happens to initialise the first row correctly, so there is no need to modify it.
Then we are ready to run the recurrence.
P = StencilRecurrencePlan{Complex}(stencil, coef, f_init, (101,101), (1,2))
ret = hcat(stable_recurrence(P)...)101×101 Matrix{ComplexF64}:
3.14159+0.0im 0.0+2.0im … 3.68182e-91+0.0im
0.0+0.0im 1.14159+0.0im 0.0+0.019998im
0.858407-0.0im 0.0+0.283185im 0.0003998+0.0im
0.0+0.0im 0.575222+0.0im 0.0+0.019982im
0.474926+0.0im 0.0+0.100296im 0.000798643+0.0im
0.0+0.0im 0.37463+0.0im … 0.0+0.0199501im
0.325074-0.0im 0.0+0.0495559im 0.00119558+0.0im
0.0+0.0im 0.275518+0.0im 0.0+0.0199025im
0.246355+0.0im 0.0+0.0291635im 0.00158968+0.0im
0.0+0.0im 0.217191+0.0im 0.0+0.0198394im
⋮ ⋱ ⋮
0.0217366+0.0im 0.0+0.000236156im 0.00996587+0.0im
0.0+0.0im 0.0215004+0.0im 0.0+0.0107245im
0.0212742-0.0im 0.0+0.000226219im 0.00998141+0.0im
0.0+0.0im 0.021048+0.0im … 0.0+0.0105125im
0.0208311+0.0im 0.0+0.000216896im 0.00999219+0.0im
0.0+0.0im 0.0206142+0.0im 0.0+0.0103045im
0.020406-0.0im 0.0+0.000208138im 0.00999847+0.0im
0.0+0.0im 0.0201979+0.0im 0.0+0.0101005im
0.019998+0.0im 0.0+0.0001999im … 0.0100005+0.0imBy design, stable_recurrence returns a vector of slices, so hcat is used to re-assemble the result. The last argument of StencilRecurrencePlan specifies the index where the recurrence starts at. By default, it will start from where the stencil doesn't go out of bound. In this example, however, we do want to apply the stencil to the first row in order to initialise it.
Now we compare the results against that from QuadGK.jl to verify the correctness:
maximum(abs, ret[92:end, 92:end] - [quadgk(x->(cos(x)/(1+sin(x)))^m*exp(im*n*x), 0, π)[1] for m in 91:100, n in 91:100])1.0048920100568746e-15The recurrence is faster than generic numerical integration when the intermediate results are useful. In this example, when m is small and n is large, the function is oscillatory which QuadGK.jl has difficulty to handle:
julia> @btime quadgk(x->(cos(x)/(1+sin(x)))^100*exp(im*100*x), 0, π);19.908 μs (3 allocations: 2.19 KiB)julia> @btime quadgk(x->(cos(x)/(1+sin(x)))^2*exp(im*100*x), 0, π);118.411 μs (4 allocations: 9.06 KiB)julia> @btime quadgk(x->exp(im*100*x), 0, π);269.570 ms (13 allocations: 23.76 MiB)
This is why we don't compare the whole matrix.
Fractional integral operator
Reference
Pu, Tianyi, and Marco Fasondini. "The numerical solution of fractional integral equations via orthogonal polynomials in fractional powers." Advances in Computational Mathematics 49, no. 1 (2023): 7.
As a preparation, we load some functions for operators:
include("../../test/fractional_integral_operator.jl")Suppose we know a banded integral operator $L$:
# 0,0,0,2 are parameters of the object space. See the reference for details.
N = 100 # truncation size
L = OpI(0,0,0,2,N)100×100 BandedMatrix{Float64} with bandwidths (2, 2):
0.666667 0.0 -0.133333 ⋅ … ⋅ ⋅
1.0 0.2 -0.2 -0.0857143 ⋅ ⋅
0.333333 0.333333 0.047619 -0.142857 ⋅ ⋅
⋅ 0.133333 0.2 0.0222222 ⋅ ⋅
⋅ ⋅ 0.0857143 0.142857 ⋅ ⋅
⋅ ⋅ ⋅ 0.0634921 … ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ … ⋅ ⋅
⋅ ⋅ ⋅ ⋅ -0.00255109 ⋅
⋅ ⋅ ⋅ ⋅ -0.00507614 -0.00252532
⋅ ⋅ ⋅ ⋅ 2.57699e-5 -0.00502513
⋅ ⋅ ⋅ ⋅ 0.00507614 2.52544e-5and we want to obtain $\sqrt{L}$ which is upper Hessenberg. We know how to compute each entry of $\sqrt{L}$:
# 0.5 is the order for sqrt.
# the last argument is truncation size.
sqrtL = OpI11(0,0,0,2,0.5,5)7×6 Matrix{Float64}:
0.797885 0.0376582 -0.112975 -0.00353121 -0.011373 -0.00177763
0.797885 0.45543 -0.0896734 -0.120037 -0.0117627 -0.0181185
0.0 0.417771 0.342455 -0.100871 -0.109833 -0.0154641
0.0 0.0 0.319154 0.284202 -0.0981327 -0.10092
0.0 0.0 0.0 0.268567 0.247721 -0.0931577
0.0 0.0 0.0 0.0 0.23641 0.222272
0.0 0.0 0.0 0.0 0.0 0.213633but it's very expensive. Fortunately, we have the banded Sylvester equation $\sqrt{L}L = L\sqrt{L}$ which is basically a stencil recurrence, so we can just compute the first few entries and use recurrence to get the rest. Further, PseudostableRecurrences.jl has BandedSylvesterRecurrence implemented, so we have a shortcut.
The BandedSylvesterRecurrence solves the equation $AX+XB+C=O$ through recurrence. In this example, $A=L$, $B=-L$, $C=O$ and $X=\sqrt{L}$.
A(T) = OpI(T(0), T(0), T(0), T(2), N+4) # computes a big enough truncation of L
B(T) = -A(T)
C(T) = Zeros{T}(∞,∞)We then set the initial condition. The stencil is a 5x5 cross so we need 5 columns for the buffer.
function init(T, N) # N is the height
ret = zeros(T, N, 5)
ret[1:3, 1:2] = OpI11(T(0), T(0), T(0), T(2), T(0.5), 1) # just need the very first results
ret
endNow we are ready to go.
P = BandedSylvesterRecurrencePlan(A, B, C, init, (N+2,N+1), (2, 2))
sqrtL = hcat(stable_recurrence(P)...)102×101 Matrix{Float64}:
0.797885 0.0376582 -0.112975 -0.00353121 … -5.56587e-7 -5.71584e-7
0.797885 0.45543 -0.0896734 -0.120037 -1.76799e-6 -1.62137e-6
0.0 0.417771 0.342455 -0.100871 -2.78554e-6 -2.8604e-6
0.0 0.0 0.319154 0.284202 -4.1314e-6 -3.78899e-6
0.0 0.0 0.0 0.268567 -5.02495e-6 -5.15916e-6
0.0 0.0 0.0 0.0 … -6.50949e-6 -5.97056e-6
0.0 0.0 0.0 0.0 -7.28327e-6 -7.47594e-6
0.0 0.0 0.0 0.0 -8.91084e-6 -8.1742e-6
0.0 0.0 0.0 0.0 -9.56916e-6 -9.81899e-6
0.0 0.0 0.0 0.0 -1.13443e-5 -1.04083e-5
⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 -0.00309099 -0.00190603
0.0 0.0 0.0 0.0 -0.00308256 -0.00307602
0.0 0.0 0.0 0.0 … -0.0062166 -0.00306771
0.0 0.0 0.0 0.0 -0.0062002 -0.0061861
0.0 0.0 0.0 0.0 -0.0249988 -0.00616995
0.0 0.0 0.0 0.0 -0.0249349 -0.0248748
0.0 0.0 0.0 0.0 0.0502511 -0.0248118
0.0 0.0 0.0 0.0 … 0.0501258 0.0499992
0.0 0.0 0.0 0.0 0.0 0.0498758where (N+2,N+1) is the size of sqrtL and (2, 2) is the bandwidth of B.
We can test that sqrtL is indeed the square root of L:
maximum(abs, sqrtL[1:N,1:N+1]*sqrtL[1:N+1,1:N] - L)1.1102230246251565e-16