FODE solvers overhaul

Project.toml

@@ -1,7 +1,7 @@
 name = "FractionalDiffEq"
 uuid = "c7492dd8-6170-483b-af64-cefb6c377d9a"
 authors = ["Qingyu Qu <[email protected]>"]
-version = "0.3.4"
+version = "0.3.5"
 
 [deps]
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"

docs/src/dode.md

@@ -14,7 +14,7 @@ We can write the general form of distributed order differential equations as:
 \int_0^m \mathscr{A}(r,\ D_*^r u(t))dr = f(t)
 ```
 
-Similar with what we have learned about single-term and multi-term fractional differential equations in linear fractional differential equations, we can also write the single-term distributed order differential equations:
+Similar to what we have learned about single-term and multi-term fractional differential equations in linear fractional differential equations, we can also write the single-term distributed order differential equations:
 
 ```math
 D_*^ru(t)=f(t,\ u(t))
@@ -28,20 +28,20 @@ And multi-term distributed order differential equations
 
 ## Example1: Distributed order relaxation
 
-The distributed order relaxation equation is similar with fractional relaxation equation, only the order is changed to a distributed function. Let's see an example here, the distributed order relaxation:
+The distributed order relaxation equation is similar to fractional relaxation equation, only the order is changed to a distributed function. Let's see an example here, the distributed order relaxation:
 
 ```math
 {_0D_t^{\omega(\alpha)} u(t)}+bu(t)=f(t),\quad x(0)=1
 ```
 
-With distribution of orders ``\alpha``: ``\omega(\alpha)=6\alpha(1-\alpha)``
+With a distribution of orders ``\alpha``: ``\omega(\alpha)=6\alpha(1-\alpha)``
 
 By using the ```DOMatrixDiscrete``` method to solve this problem:
 
 !!! info
-    The usage of ```DOMatrixDiscrete``` method is similiar with the ```FODEMatrixDiscrete``` method, all we need to do is to pass the parameters array and orders array to the problem difinition and solve the problem.
+    The usage of ```DOMatrixDiscrete``` method is similar to the ```FODEMatrixDiscrete``` method, all we need to do is to pass the parameters array and order array to the problem definition and solve the problem.
     !!! tip
-        Pass the weight function and other orders to the order array is the right way:
+        Pass the weight function and other orders to the order array in the right way:
         ```julia-repl
         julia> orders = [x->x*(1-x), 1.2, 3]
         3-element Vector{Any}:

docs/src/get_started.md

@@ -25,17 +25,17 @@ We can solve this problem by the following code using FractionalDiffEq.jl:
 ```julia
 using FractionalDiffEq, Plots
 fun(u, p, t) = 1-u
-α=1.8; h=0.01; tspan = (0, 20); u0 = [0, 0]
-prob = SingleTermFODEProblem(fun, α, u0, tspan)
-sol = solve(prob, h, PECE())
+α=1.8; tspan = (0, 20); u0 = [0, 0]
+prob = FODEProblem(fun, α, u0, tspan)
+sol = solve(prob, PIEX(), dt=0.01)
 plot(sol)
 ```
 
 By plotting the numerical result, we can get the approximation result:
 
 ![Relaxation Oscillation](./assets/example.png)
 
-To provide users a simple way to solve fractional differential equations, we follow the design pattern of [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl)
+To provide users with a simple way to solve fractional differential equations, we follow the design pattern of [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl)
 
 ## Step 1: Defining a Problem
 
@@ -45,10 +45,10 @@ First, we need to specify the problem we want to solve. Just by passing the para
 using FractionalDiffEq
 fun(u, p, t) = 1-u
 α = 1.8; u0 = [0, 0]; tspan = (0, 20); h = 0.01;
-prob = SingleTermFODEProblem(fun, α, u0, tspan)
+prob = FODEProblem(fun, α, u0, tspan)
 ```
 
-The ```SingleTermFODEProblem``` is a class of fractional differential equation, describing equations with ``D^{\alpha}u=f(t, u)`` pattern. For other patterns and classes of fractional differential equation, please refer to [Problem types](@ref problems)
+The ```FODEProblem``` is a class of fractional differential equations, describing equations with ``D^{\alpha}u=f(t, u)`` pattern. For other patterns and classes of fractional differential equations, please refer to [Problem types](@ref problems)
 
 ## Step 2: Solving a Problem
 

src/FractionalDiffEq.jl

@@ -36,6 +36,8 @@ include("fode/bdf.jl")
 include("fode/newton_gregory.jl")
 include("fode/trapezoid.jl")
 include("fode/explicit_pi.jl")
+include("fode/implicit_pi_rectangle.jl")
+include("fode/implicit_pi_trapzoid.jl")
 include("fode/grunwald_letnikov.jl")
 include("fode/NonLinear.jl")
 include("fode/newton_polynomials.jl")
@@ -103,12 +105,14 @@ export FODESolution, FDifferenceSolution, DODESolution, FFMODESolution
 export FODESystemSolution, FDDESystemSolution, FFMODESystem
 
 # FODE solvers
-export PIPECE, PIRect, PITrap, MTPIEX
-export PECE, MatrixDiscrete, GL
+export PIRect, PITrap, PECE, PIEX
+export MatrixDiscrete, GL
 export AtanganaSedaAB
 
+export MTPIRect, MTPITrap, MTPECE, MTPIEX
+
 # System of FODE solvers
-export NonLinearAlg, BDF, NewtonGregory, Trapezoid, PIEX, NewtonPolynomial
+export NonLinearAlg, BDF, NewtonGregory, Trapezoid, NewtonPolynomial
 export AtanganaSedaCF
 export AtanganaSeda
 

src/fode/bdf.jl

@@ -30,6 +30,7 @@
     reltol
     abstol
     maxiters
+    high_order_prob
 
     kwargs
 end
@@ -47,13 +48,14 @@
 
     all(x->x==order[1], order) ? nothing : throw(ArgumentError("BDF method is only for commensurate order FODE"))
     alpha = order[1] # commensurate ordre FODE
-    (alpha > 1.0) && throw(ArgumentError("BDF method is only for order <= 1.0"))
-
+
 
     m_alpha = ceil.(Int, alpha)
     m_alpha_factorial = factorial.(collect(0:m_alpha-1))
-    problem_size = length(u0)
-
+    problem_size = length(order)
+    u0_size = length(u0)
+    high_order_prob = problem_size !== u0_size
+
     # Number of points in which to evaluate the solution or the BDF_weights
     r = 16
     N = ceil(Int, (tfinal-t0)/dt)
@@ -66,7 +68,7 @@
     fy = zeros(T, problem_size, N+1)
     zn = zeros(T, problem_size, NNr+1)
 
-    # generate jacobian of input function
+    # generate jacobian of the input function
     Jfdefun(t, u) = jacobian_of_fdefun(prob.f, t, u, p)
 
     # Evaluation of convolution and starting BDF_weights of the FLMM
@@ -75,13 +77,14 @@
 
     # Initializing solution and proces of computation
     mesh = t0 .+ collect(0:N)*dt
-    y[:, 1] .= u0
-    temp = similar(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
+
     return BDFCache{iip, T}(prob, alg, mesh, u0, alpha, halpha, y, fy, zn, Jfdefun,
                             p, problem_size, m_alpha, m_alpha_factorial, r, N, Nr, Q, NNr,
-                            omega, w, s, dt, reltol, abstol, maxiters, kwargs)
+                            omega, w, s, dt, reltol, abstol, maxiters, high_order_prob, kwargs)
 end
 
 function SciMLBase.solve!(cache::BDFCache{iip, T}) where {iip, T}
@@ -158,7 +161,7 @@
     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}
+function BDF_triangolo(cache::BDFCache{iip, T}, nxi::P, nxf::P, j0) where {P <: Integer, iip, T}
     @unpack prob, mesh, problem_size, zn, Jfdefun, N, abstol, maxiters, s, w, omega, halpha, u0, m_alpha, m_alpha_factorial, p = cache
     for n = nxi:min(N, nxf)
         n1 = n+1
@@ -315,11 +318,12 @@
 Jf_vectorfield(t, y, Jfdefun) = Jfdefun(t, y)
 
 function ABM_starting_term(cache::BDFCache{iip, T}, t) where {iip, T}
-    @unpack u0, m_alpha, mesh, m_alpha_factorial = cache
+    @unpack 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)
     for k = 1:m_alpha
-        ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*u0
+        ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*u0[:, k]
     end
     return ys
 end
@@ -336,13 +340,17 @@
 
 function _convert_single_term_to_vectorized_prob!(prob::FODEProblem)
     if SciMLBase.isinplace(prob)
-        new_prob = remake(prob; u0=[prob.u0], order=[prob.order])
+        if isa(prob.u0, AbstractArray)
+            new_prob = remake(prob; order=[prob.order])
+        else
+            new_prob = remake(prob; u0=[prob.u0], order=[prob.order])
+        end
         return new_prob
     else
         function new_f(du, u, p, t)
             du[1] = prob.f(u[1], p, t)
         end
-        new_fun = ODEFunction(new_f)
+        new_fun = ODEFunction{true}(new_f) # make in-place
         new_prob = remake(prob; f=new_fun, u0=[prob.u0], order=[prob.order])
         return new_prob
     end

src/fode/explicit_pi.jl

@@ -9,7 +9,7 @@
     y
     fy
     p
-    problem_size    
+    problem_size
     zn
 
     r
@@ -247,13 +247,13 @@ end
 function  PIEX_system_starting_term(cache::PIEXCache{iip, T}, t) where {iip, T}
     @unpack mesh, u0, m_alpha, m_alpha_factorial = cache
     t0 = mesh[1]
-    ys = zeros(size(u0, 1), 1)
+    ys = zeros(length(u0))
     for k = 1 : maximum(m_alpha)
         if length(m_alpha) == 1
-            ys = ys + (t-t0)^(k-1)/m_alpha_factorial[k]*u0[:, k]
+            ys = ys .+ (t-t0)^(k-1)/m_alpha_factorial[k]*u0[k]
         else
             i_alpha = findall(x -> x>=k, m_alpha)
-            ys[i_alpha, 1] = ys[i_alpha, 1] + (t-t0)^(k-1)*u0[i_alpha, k]./m_alpha_factorial[i_alpha, k]
+            ys[i_alpha] = ys[i_alpha] + (t-t0)^(k-1)*u0[i_alpha, k]./m_alpha_factorial[i_alpha, k]
         end
     end
     return ys