Skip to content

Commit

Permalink
Initial improvement for FODE solvers
Browse files Browse the repository at this point in the history 
Signed-off-by: Qingyu Qu <2283984853@qq.com>
  • Loading branch information
@ErikQQY
ErikQQY committed Oct 9, 2024
1 parent acf5643 commit 74382af
Show file tree
Hide file tree
Showing 11 changed files with 253 additions and 523 deletions.
6 changes: 5 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,20 @@
name = "FractionalDiffEq"
uuid = "c7492dd8-6170-483b-af64-cefb6c377d9a"
authors = ["Qingyu Qu <erikqqy123@gmail.com>"]
version = "0.4.0"
version = "0.5.0"

[deps]
ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
InvertedIndices = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Expand All @@ -23,12 +25,14 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
[compat]
ConcreteStructs = "0.2.2"
DiffEqBase = "6.156"
FastClosures = "0.3.2"
FFTW = "1.8.0"
ForwardDiff = "0.10.36"
HypergeometricFunctions = "0.3.23"
InvertedIndices = "1.3"
Polynomials = "3, 4"
RecipesBase = "1.3.4"
RecursiveArrayTools = "3.27"
Reexport = "1.2.2"
SciMLBase = "2.55"
SparseArrays = "1.10"
Expand Down
6 changes: 2 additions & 4 deletions src/FractionalDiffEq.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module FractionalDiffEq

using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, ToeplitzMatrices,
FFTW, RecipesBase, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs
using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, ToeplitzMatrices, RecursiveArrayTools,
FFTW, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs, FastClosures

import SciMLBase: __solve
import DiffEqBase: solve
Expand All @@ -16,8 +16,6 @@ include("types/problems.jl")
include("types/algorithms.jl")
include("types/solutions.jl")

include("types/problem_utils.jl")

# Multi-terms fractional ordinary differential equations
include("multitermsfode/matrix.jl")
include("multitermsfode/hankelmatrix.jl")
Expand Down
175 changes: 90 additions & 85 deletions src/fode/bdf.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
y
fy
zn
Jfdefun
jac
p
problem_size

Expand Down Expand Up @@ -65,55 +65,63 @@ function SciMLBase.__init(prob::FODEProblem, alg::BDF; dt = 0.0, reltol = 1e-6,
NNr = 2^(Q + 1) * r

# Preallocation of some variables
y = zeros(T, problem_size, N + 1)
fy = zeros(T, problem_size, N + 1)
y = [Vector{T}(undef, problem_size) for _ in 1:(N + 1)]
fy = similar(y)
zn = zeros(T, problem_size, NNr + 1)

# generate jacobian of the input function
Jfdefun(t, u) = jacobian_of_fdefun(prob.f, t, u, p)
if prob.f.jac === nothing
if iip
jac = (df, u, p, t) -> begin
_du = similar(u)
prob.f(_du, u, p, t)
_f = @closure (du, u) -> prob.f(du, u, p, t)
ForwardDiff.jacobian!(df, _f, _du, u)
return
end
else
jac = (df, u, p, t) -> begin
_du = prob.f(u, p, t)
_f = @closure (du, u) -> (du .= prob.f(u, p, t))
ForwardDiff.jacobian!(df, _f, _du, u)
return
end
end
else
jac = prob.f.jac
end

# Evaluation of convolution and starting BDF_weights of the FLMM
(omega, w, s) = BDF_weights(alpha, NNr + 1)
halpha = dt^alpha

# Initializing solution and proces of computation
mesh = t0 .+ collect(0:N) * dt
y[:, 1] = high_order_prob ? u0[1, :] : u0
y[1] .= high_order_prob ? u0[1, :] : u0
temp = high_order_prob ? similar(u0[1, :]) : similar(u0)
f(temp, u0, p, t0)
fy[:, 1] = temp
prob.f(temp, u0, p, t0)
fy[1] = temp

return BDFCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, Jfdefun, p,
return BDFCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, jac, prob.p,
problem_size, m_alpha, m_alpha_factorial, r, N, Nr, Q, NNr, omega,
w, s, dt, reltol, abstol, maxiters, high_order_prob, kwargs)
end

function SciMLBase.solve!(cache::BDFCache{iip, T}) where {iip, T}
(; prob, alg, mesh, y, r, N, Q, s, dt) = cache
tfinal = mesh[end]

BDF_first_approximations(cache)
BDF_triangolo(cache, s + 1, r - 1, 0)

# Main process of computation by means of the FFT algorithm
nx0 = 0
ny0 = 0
ff = zeros(T, 1, 2^(Q + 2), 1)
ff[1:2] = [0 2]
ff = zeros(T, 1, 2^(Q + 2))
ff[1:2] = [0; 2]
for q in 0:Q
L = 2^q
BDF_disegna_blocchi(cache, L, ff, nx0 + L * r, ny0)
ff[1:(4 * L)] = [ff[1:(2 * L)]; ff[1:(2 * L - 1)]; 4 * L]
end
# Evaluation solution in TFINAL when TFINAL is not in the mesh
if tfinal < mesh[N + 1]
c = (tfinal - mesh[N]) / dt
mesh[N + 1] = tfinal
y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]
end
mesh = mesh[1:(N + 1)]
y = y[:, 1:(N + 1)]
y = collect(Vector{eltype(y)}, eachcol(y))

return DiffEqBase.build_solution(prob, alg, mesh, y)
end
Expand All @@ -122,28 +130,28 @@ function BDF_disegna_blocchi(
cache::BDFCache{iip, T}, L::P, ff, nx0::P, ny0::P) where {P <: Integer, iip, T}
(; r, Nr, N) = cache

nxi::Int = copy(nx0)
nxf::Int = copy(nx0 + L * r - 1)
nyi::Int = copy(ny0)
nyf::Int = copy(ny0 + L * r - 1)
nxi::Int = nx0
nxf::Int = nx0 + L * r - 1
nyi::Int = ny0
nyf::Int = ny0 + L * r - 1
is = 1
s_nxi = zeros(T, N)
s_nxf = zeros(T, N)
s_nyi = zeros(T, N)
s_nyf = zeros(T, N)
s_nxi[is] = nxi
s_nxf[is] = nxf
s_nyi[is] = nyi
s_nyf[is] = nyf
s_nxi = Vector{T}(undef, N)
s_nxf = similar(s_nxi)
s_nyi = similar(s_nxi)
s_nyf = similar(s_nxi)
s_nxi[1] = nxi
s_nxf[1] = nxf
s_nyi[1] = nyi
s_nyf[1] = nyf
i_triangolo = 0
stop = false
while ~stop
while !stop
stop = (nxi + r - 1 == nx0 + L * r - 1) || (nxi + r - 1 >= Nr - 1)
BDF_quadrato(cache, nxi, nxf, nyi, nyf)
BDF_triangolo(cache, nxi, nxi + r - 1, nxi)
i_triangolo = i_triangolo + 1

if ~stop
if !stop
if nxi + r - 1 == nxf
i_Delta = ff[i_triangolo]
Delta = i_Delta * r
Expand Down Expand Up @@ -178,81 +186,73 @@ function BDF_quadrato(cache::BDFCache, nxi::P, nxf::P, nyi::P, nyf::P) where {P
funz_beg = nyi + 1
funz_end = nyf + 1
vett_coef = omega[(coef_beg + 1):(coef_end + 1)]
vett_funz = [cache.fy[:, funz_beg:funz_end] zeros(problem_size, funz_end - funz_beg + 1)]
vett_funz = [reduce(hcat, cache.fy[funz_beg:funz_end]) zeros(problem_size, funz_end - funz_beg + 1)]
zzn = real(fast_conv(vett_coef, vett_funz))
cache.zn[:, (nxi + 1):(nxf + 1)] = cache.zn[:, (nxi + 1):(nxf + 1)] +
zzn[:, (nxf - nyf):(end - 1)]
end

function BDF_triangolo(
cache::BDFCache{iip, T}, nxi::P, nxf::P, j0) where {P <: Integer, iip, T}
(; prob, mesh, problem_size, zn, Jfdefun, N, abstol, maxiters, s, w, omega, halpha, p) = cache
(; prob, mesh, problem_size, zn, jac, N, abstol, maxiters, s, w, omega, halpha, p) = cache
Jfn0 = Matrix{T}(undef, problem_size, problem_size)
for n in nxi:min(N, nxf)
n1 = n + 1
St = ABM_starting_term(cache, mesh[n1])

Phi = zeros(T, problem_size, 1)
for j in 0:s
Phi = Phi + w[j + 1, n1] * cache.fy[:, j + 1]
Phi = Phi + w[j + 1, n1] * cache.fy[j + 1]
end
for j in j0:(n - 1)
Phi = Phi + omega[n - j + 1] * cache.fy[:, j + 1]
Phi = Phi + omega[n - j + 1] * cache.fy[j + 1]
end
Phi_n = St + halpha * (zn[:, n1] + Phi)

yn0 = cache.y[:, n]
temp = zeros(T, problem_size)
prob.f(temp, yn0, p, mesh[n1])
fn0 = temp
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
yn0 = cache.y[n]
fn0 = similar(yn0)
prob.f(fn0, yn0, p, mesh[n1])
@views jac(Jfn0, yn0, p, mesh[n1])
Gn0 = yn0 - halpha * omega[1] * fn0 - Phi_n
stop = false
it = 0
yn1 = similar(yn0)
fn1 = similar(yn0)
while ~stop
while !stop
JGn0 = zeros(T, problem_size, problem_size) + I - halpha * omega[1] * Jfn0
yn1 = yn0 - JGn0 \ Gn0
yn1 = yn0 - vec(JGn0 \ Gn0)
prob.f(fn1, yn1, p, mesh[n1])
Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
it = it + 1

stop = (norm(yn1 - yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
if it > maxiters && ~stop
if it > maxiters && !stop
@warn "Non Convergence"
stop = true
end

yn0 = yn1
Gn0 = Gn1
if ~stop
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
if !stop
@views jac(Jfn0, yn0, p, mesh[n1])
end
end
cache.y[:, n1] = yn1
cache.fy[:, n1] = fn1
cache.y[n1] = yn1
cache.fy[n1] = fn1
end
end

function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
(; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, Jfdefun, p) = cache
(; prob, mesh, abstol, problem_size, maxiters, s, halpha, omega, w, jac, p) = cache

Im = zeros(problem_size, problem_size) + I
Ims = zeros(problem_size * s, problem_size * s) + I
Y0 = zeros(T, s * problem_size, 1)
F0 = copy(Y0)
B0 = copy(Y0)
Y0 = VectorOfArray([cache.y[1] for _ in 1:s])
F0 = similar(Y0)
B0 = similar(Y0)
for j in 1:s
Y0[((j - 1) * problem_size + 1):(j * problem_size), 1] = cache.y[:, 1]
temp = zeros(length(cache.y[:, 1]))
prob.f(temp, cache.y[:, 1], p, mesh[j + 1])
F0[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
prob.f(F0.u[j], cache.y[1], p, mesh[j + 1])
St = ABM_starting_term(cache, mesh[j + 1])
B0[((j - 1) * problem_size + 1):(j * problem_size), 1] = St +
halpha *
(omega[j + 1] +
w[1, j + 1]) *
cache.fy[:, 1]
B0.u[j] = St + halpha * (omega[j + 1] + w[1, j + 1]) * cache.fy[1]
end
W = zeros(T, s, s)
for i in 1:s
Expand All @@ -265,32 +265,29 @@ function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
end
end
W = halpha * kron(W, Im)
G0 = Y0 - B0 - W * F0
G0 = vec(Y0 - B0) - W * vec(F0)
JF = zeros(T, s * problem_size, s * problem_size)
J_temp = Matrix{T}(undef, problem_size, problem_size)
for j in 1:s
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
mesh[j + 1], cache.y[:, 1], Jfdefun)
jac(J_temp, cache.y[1], p, mesh[j + 1])
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
end
stop = false
it = 0
F1 = zeros(T, s * problem_size, 1)
F1 = similar(F0)
Y1 = similar(Y0)
while ~stop
while !stop
JG = Ims - W * JF
Y1 = Y0 - JG \ G0
recursive_unflatten!(Y1, vec(Y0) - JG \ G0)

for j in 1:s
temp = zeros(T, length(Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]))
prob.f(temp, Y1[((j - 1) * problem_size + 1):(j * problem_size), 1],
p, mesh[j + 1])
F1[((j - 1) * problem_size + 1):(j * problem_size), 1] = temp
prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
end
G1 = Y1 - B0 - W * F1
G1 = vec(Y1 - B0) - W * vec(F1)

it = it + 1

stop = (norm(Y1 - Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > maxiters && ~stop
if it > maxiters && !stop
@warn "Non Convergence"
stop = 1
end
Expand All @@ -299,15 +296,14 @@ function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
G0 = G1
if ~stop
for j in 1:s
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
mesh[j + 1],
Y1[((j - 1) * problem_size + 1):(j * problem_size), 1], Jfdefun)
jac(J_temp, Y1.u[j], p, mesh[j+1])
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
end
end
end
for j in 1:s
cache.y[:, j + 1] = Y1[((j - 1) * problem_size + 1):(j * problem_size), 1]
cache.fy[:, j + 1] = F1[((j - 1) * problem_size + 1):(j * problem_size), 1]
cache.y[j + 1] = Y1.u[j]
cache.fy[j + 1] = F1.u[j]
end
end

Expand Down Expand Up @@ -347,7 +343,7 @@ function BDF_weights(alpha, N)
nn .^ (nu + alpha)
else
if i == 0
nn_nu_alpha[i + 1, :] = zeros(1, N + 1)
nn_nu_alpha[i + 1, :] = zeros(N + 1)
else
nn_nu_alpha[i + 1, :] = gamma(nu + 1) / gamma(nu + 1 + alpha) *
nn .^ (nu + alpha)
Expand All @@ -369,7 +365,7 @@ function ABM_starting_term(cache::BDFCache{iip, T}, t) where {iip, T}
(; u0, m_alpha, mesh, m_alpha_factorial, high_order_prob) = cache
t0 = mesh[1]
u0 = high_order_prob ? reshape(u0, 1, length(u0)) : u0
ys = zeros(size(u0, 1), 1)
ys = zeros(size(u0, 1))
for k in 1:m_alpha
ys = ys + (t - t0)^(k - 1) / m_alpha_factorial[k] * u0[:, k]
end
Expand Down Expand Up @@ -405,3 +401,12 @@ function _convert_single_term_to_vectorized_prob!(prob::FODEProblem)
return new_prob
end
end

@views function recursive_unflatten!(y::VectorOfArray, x::AbstractArray)
i = 0
for yᵢ in y
copyto!(yᵢ, x[(i + 1):(i + length(yᵢ))])
i += length(yᵢ)
end
return nothing
end
Loading

0 comments on commit 74382af

Please sign in to comment.