Merge pull request #101 from ErikQQY/qqy/refactor_mtfde

Refactor multi-terms FODE solvers

docs/src/algorithms.md

@@ -16,7 +16,7 @@ AtanganaSeda
 
 ### Multi-Term FODE
 
-FODEMatrixDiscrete
+MatrixDiscrete
 ClosedForm
 ClosedFormHankelM
 PIPECE

src/FractionalDiffEq.jl

@@ -67,6 +67,10 @@ include("mlfun.jl")
 include("utils.jl")
 include("auxiliary.jl")
 
+function __solve(prob::MultiTermsFODEProblem, alg::MultiTermsFODEAlgorithm, args...; kwargs...)
+    cache = init(prob, alg, args...; kwargs...)
+    return solve!(cache)
+end
 
 function __solve(prob::FODEProblem, alg::FODESystemAlgorithm, args...; kwargs...)
     cache = init(prob, alg, args...; kwargs...)
@@ -94,8 +98,8 @@ export FODESolution, FDifferenceSolution, DODESolution, FFMODESolution
 export FODESystemSolution, FDDESystemSolution, FFMODESystem
 
 # FODE solvers
-export PIPECE, PIRect, PITrap
-export PECE, FODEMatrixDiscrete, ClosedForm, ClosedFormHankelM, GL
+export PIPECE, PIRect, PITrap, MTPIEX
+export PECE, MatrixDiscrete, GL
 export AtanganaSedaAB
 
 # System of FODE solvers

src/fodesystem/bdf.jl

@@ -77,11 +77,9 @@ function SciMLBase.__init(prob::FODEProblem, alg::BDF;
                             omega, w, s, dt, reltol, abstol, maxiters, kwargs)
 end
 
-function SciMLBase.solve!(cache::BDFCache)
+function SciMLBase.solve!(cache::BDFCache{iip, T}) where {iip, T}
     @unpack prob, alg, mesh, u0, order, halpha, y, fy, zn, p, m_alpha, m_alpha_factorial, r, N, Nr, Q, NNr, omega, w, s, dt, reltol, abstol, maxiters, kwargs = cache
     t0 = mesh[1]; tfinal = mesh[end]
-    iip = isinplace(prob)
-    T = eltype(cache)
     problem_size = length(u0)
     # generate jacobian of input function
     Jfdefun(t, u) = jacobian_of_fdefun(prob.f, t, u, p)

src/multitermsfode/explicit_pi.jl

@@ -1,79 +1,128 @@
-function solve(prob::MultiTermsFODEProblem, h, ::PIEX)
-    @unpack parameters, orders, rightfun, u0, tspan = prob
-    t0 = tspan[1]; T = tspan[2]
+@concrete mutable struct MTPIEXCache{T}
+    prob
+    alg
+    mesh
+    u0
+    bet
+    lam_rat_i
+    gamma_val
+
+    highest_order_parameter
+    highest_order_ceiling
+    other_orders_ceiling
+
+    y
+    fy
+    p
+
+    zn
+
+    r
+    N
+    Nr
+    Qr
+    NNr
+    bn
+
+
+    kwargs
+end
+
+Base.eltype(::MTPIEXCache{T}) where {T} = T
+
+function SciMLBase.__init(prob::MultiTermsFODEProblem, alg::MTPIEX; dt = 0.0, abstol = 1e-6, kwargs...)
+    dt ≤ 0 ? throw(ArgumentError("dt must be positive")) : nothing
+    @unpack parameters, orders, f, u0, tspan, p = prob
+    t0 = tspan[1]; tfinal = tspan[2]
+    T = eltype(u0)
     u0 = u0[:]'
-    Q = length(orders)
+    orders_length = length(orders)
     orders= sort(orders)
-    i_al = sortperm(orders, rev=true)
-    parameters = parameters[i_al]
-    al_Q = orders[end]
-    al_i = orders[1:end-1]
-    lam_Q = parameters[end]
-    lam_rat_i = parameters[1:end-1]/lam_Q
-    m_Q = ceil(Int64, al_Q)
-    m_i = ceil.(Int64, orders[1:end-1])
-    bet = [al_Q .- al_i; al_Q]
+    sorted_parameters_index = sortperm(orders, rev=true)
+    parameters = parameters[sorted_parameters_index]
+    highest_order = orders[end]
+    other_orders = orders[1:end-1]
+    highest_order_parameter = parameters[end]
+    lam_rat_i = parameters[1:end-1]/highest_order_parameter
+    highest_order_ceiling = ceil(Int64, highest_order)
+    other_orders_ceiling = ceil.(Int64, orders[1:end-1])
+    bet = [highest_order .- other_orders; highest_order]
 
-    gamma_val = zeros(Q, m_Q)
-    for i = 1 : Q-1
-        k = collect(Int, 0:m_i[i]-1)
+    gamma_val = zeros(orders_length, highest_order_ceiling)
+    for i = 1 : orders_length-1
+        k = collect(Int, 0:other_orders_ceiling[i]-1)
         gamma_val[i, k.+1] = gamma.(k.+bet[i].+1)
     end
-    k = collect(0:m_Q-1)
-    gamma_val[Q, :] = factorial.(k)
+    k = collect(0:highest_order_ceiling-1)
+    gamma_val[orders_length, :] = factorial.(k)
 
 
     problem_size = size(u0, 1)
 
     r::Int = 16
-    N::Int = ceil(Int, (T-t0)/h)
+    N::Int = ceil(Int, (tfinal-t0)/dt)
     Nr::Int = ceil(Int, (N+1)/r)*r
     Qr::Int = ceil(Int, log2((Nr)/r))-1
     NNr::Int = 2^(Qr+1)*r
 
     y = zeros(Float64, problem_size, N+1)
     fy = zeros(Float64, problem_size, N+1)
-    zn = zeros(Float64, problem_size, NNr+1, Q)
+    zn = zeros(Float64, problem_size, NNr+1, orders_length)
 
     nvett = collect(0:NNr+1)
-    bn = zeros(Q, NNr+1)
-    for i = 1:Q
+    bn = zeros(orders_length, NNr+1)
+    for i = 1:orders_length
         nbeta = nvett.^bet[i]
-        bn[i, :] = (nbeta[2:end] - nbeta[1:end-1]) * h^bet[i]/gamma(bet[i]+1)
+        bn[i, :] = (nbeta[2:end] - nbeta[1:end-1]) * dt^bet[i]/gamma(bet[i]+1)
     end
     C = 0
-    for i = 1:Q-1
+    for i = 1:orders_length-1
         C = C + lam_rat_i[i]*bn[i, 1]
     end
 
-    t = collect(0:N)*h
+    mesh = collect(0:N)*dt
     y[:, 1] = u0[:, 1]
-    fy[:, 1] .= rightfun(t0, u0[:, 1])
-    (y, fy) = PIEX_multiterms_triangolo(1, r-1, t, y, fy, zn, N, bn, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q)
+    fy[:, 1] .= f(u0[:], p, t0)
+
+    return MTPIEXCache{T}(prob, alg, mesh, u0, bet, lam_rat_i, gamma_val,
+                               highest_order_parameter, highest_order_ceiling, other_orders_ceiling,
+                               y, fy, p, zn,
+                               r, N, Nr, Qr, NNr, bn, kwargs)
+end
+function SciMLBase.solve!(cache::MTPIEXCache{T}) where {T}
+    @unpack prob, alg, mesh, u0, bet, lam_rat_i, gamma_val,
+            highest_order_parameter, highest_order_ceiling, other_orders_ceiling,
+            y, fy, p, zn,
+            r, N, Nr, Qr, NNr, bn, kwargs = cache
+    t0 = mesh[1]; tfinal = mesh[end]
+    PIEX_multiterms_triangolo(cache, 1, r-1, t0)
 
     ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0 2]; card_ff = 2
     nx0 = 0; nu0 = 0
     for qr = 0 : Qr
         L = 2^qr
-        (y, fy) = PIEX_multiterms_disegna_blocchi(L, ff, r, Nr, nx0+L*r, nu0, t, y, fy, zn, N, bn, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q) ;
+        PIEX_multiterms_disegna_blocchi(cache, L, ff, nx0+L*r, nu0, t0)
         ff[1:2*card_ff] = [ff[1:card_ff] ff[1:card_ff]]
         card_ff = 2*card_ff
         ff[card_ff] = 4*L
     end
 
-    if T<t[N+1]
-        c = (T - t[N])/h
-        t[N+1] = tfinal
+    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
 
-    t = t[1:N+1]; y = y[:, 1:N+1]
-    return FODESolution(t, y[:])
+    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
 
 
-function PIEX_multiterms_disegna_blocchi(L, ff, r, Nr, nx0, nu0, t, y, fy, zn, N, bn, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q)
-
+function PIEX_multiterms_disegna_blocchi(cache::MTPIEXCache{T}, L, ff, nx0, nu0, t0) where {T}
+    @unpack prob, alg, mesh, u0, bet, lam_rat_i, gamma_val,
+            highest_order_parameter, highest_order_ceiling, other_orders_ceiling,
+            p, zn, r, N, Nr, Qr, NNr, bn, kwargs = cache
     nxi::Int = nx0
     nxf::Int = nx0 + L*r - 1
     nyi::Int = nu0
@@ -92,9 +141,9 @@ function PIEX_multiterms_disegna_blocchi(L, ff, r, Nr, nx0, nu0, t, y, fy, zn, N
     while stop == false
         stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1 >= Nr-1)
 
-        zn = PIEX_multiterms_quadrato(nxi, nxf, nyi, nyf, y, fy, zn, bn, problem_size, Q)
+        PIEX_multiterms_quadrato(cache, nxi, nxf, nyi, nyf)
 
-        (y, fy) = PIEX_multiterms_triangolo(nxi, nxi+r-1, t, y, fy, zn, N, bn, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q) ;
+        PIEX_multiterms_triangolo(cache, nxi, nxi+r-1, t0)
         i_triangolo = i_triangolo + 1
 
         if stop==false
@@ -116,78 +165,69 @@ function PIEX_multiterms_disegna_blocchi(L, ff, r, Nr, nx0, nu0, t, y, fy, zn, N
                 s_nyf[is] = nyf
             end
         end
-
     end
-    return y, fy
 end
 
-function PIEX_multiterms_quadrato(nxi, nxf, nyi, nyf, y, fy, zn, bn,  problem_size, Q)
+function PIEX_multiterms_quadrato(cache::MTPIEXCache{T}, nxi, nxf, nyi, nyf) where {T}
+    @unpack prob, alg, mesh, r, N, Nr, Qr, NNr, bn, kwargs = cache
+    orders_length = length(prob.orders)
+    problem_size = size(prob.u0, 2)
     coef_beg = nxi-nyf; coef_end = nxf-nyi+1
     funz_beg = nyi+1; funz_end = nyf+1
 
-    for i = 1:Q
+    for i = 1:orders_length
         vett_coef = bn[i, coef_beg:coef_end]
-        if i < Q
-            vett_funz = [y[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
+        if i < orders_length
+            vett_funz = [cache.y[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
         else
-            vett_funz = [fy[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
+            vett_funz = [cache.fy[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
         end
         zzn = real.(fast_conv(vett_coef, vett_funz))
-        zn[:, nxi+1:nxf+1, i] = zn[:, nxi+1:nxf+1, i] + zzn[:, nxf-nyf:end-1]
+        cache.zn[:, nxi+1:nxf+1, i] = cache.zn[:, nxi+1:nxf+1, i] + zzn[:, nxf-nyf:end-1]
     end
-    return zn
 end
-
-
-
 
-function PIEX_multiterms_triangolo(nxi, nxf, t, y, fy, zn, N, bn, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q)
-
+function PIEX_multiterms_triangolo(cache::MTPIEXCache{T}, nxi, nxf, t0) where {T}
+    @unpack prob, alg, mesh, u0, bet, lam_rat_i, gamma_val,
+            highest_order_parameter, highest_order_ceiling, other_orders_ceiling,
+            p, zn, r, N, Nr, Qr, NNr, bn, kwargs = cache
+    problem_size = size(prob.u0, 2)
+    orders_length = length(prob.orders)
     for n = nxi:min(N, nxf)
-        Phi_n = PIEX_multiterms_starting_term(t[n+1], t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val)
+        Phi_n = PIEX_multiterms_starting_term(mesh[n+1], t0, problem_size, u0, orders_length, highest_order_ceiling, other_orders_ceiling, bet, lam_rat_i, gamma_val)
         if nxi == 1
             j_beg = 0
         else
             j_beg = nxi
         end
 
-        for i = 1:Q-1
+        for i = 1:orders_length-1
             temp = zn[:, n+1, i]
             for j = j_beg:n-1
-                temp = temp + bn[i, n-j]*y[:, j+1]
+                temp = temp + bn[i, n-j]*cache.y[:, j+1]
             end
             Phi_n = Phi_n - lam_rat_i[i]*temp
         end
-        temp = zn[:, n+1, Q]
+        temp = zn[:, n+1, orders_length]
         for j = j_beg : n-1
-            temp = temp + bn[Q, n-j]*fy[:, j+1]
+            temp = temp + bn[orders_length, n-j]*cache.fy[:, j+1]
         end
-        Phi_n = Phi_n + temp/lam_Q
-
-        y[:, n+1] = Phi_n
-        fy[:, n+1] .= rightfun(t[n+1], y[:, n+1]) 
-
-    end
-    return y, fy
-end
+        Phi_n = Phi_n + temp/highest_order_parameter
 
-function process_rightfun(t, y, rightfun)
-    if typeof(rightfun) <: Function
-        return rightfun(t, y)
-    else
-        return rightfun*ones(length(y))
+        cache.y[:, n+1] = Phi_n
+        cache.fy[:, n+1] .= prob.f(cache.y[:, n+1], p, mesh[n+1])
     end
 end
 
-function PIEX_multiterms_starting_term(t, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val)
+function PIEX_multiterms_starting_term(t, t0, problem_size, u0, orders_length, highest_order_ceiling, other_orders_ceiling, bet, lam_rat_i, gamma_val)
     ys = zeros(problem_size)
 
-    for k = 0:m_Q-1
-        ys = ys .+ (t-t0)^k./gamma_val[Q, k+1]*u0[:, k+1]
+    for k = 0:highest_order_ceiling-1
+        ys = ys .+ (t-t0)^k./gamma_val[orders_length, k+1]*u0[:, k+1]
     end
-    for i = 1 : Q-1
+    for i = 1 : orders_length-1
         temp = zeros(problem_size)
-        for k = 0 : m_i[i]-1
+        for k = 0 : other_orders_ceiling[i]-1
             temp = temp .+ (t-t0)^(k+bet[i])/gamma_val[i, k+1]*u0[:, k+1]
         end
         ys = ys + lam_rat_i[i]*temp