Added matrix form for linear FDE to ODE and new test case

Author: fchen121

src/systems/diffeqs/basic_transformations.jl

@@ -342,14 +342,14 @@ prob = ODEProblem(sys, [], tspan)
 sol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
 ```
 """
-function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initials = 0, symbol = :x, iv = only(@independent_variables t))
+function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initials = 0, symbol = :x, iv = only(@independent_variables t), matrix=false)
     previous = Symbol(symbol, :_, 0)
     previous = ModelingToolkit.unwrap(only(@variables $previous(iv)))
     @variables x_0(iv)
     D = Differential(iv)
     i = 0
     all_eqs = Equation[]
-    all_def = Pair{Num, Int64}[]
+    all_def = Pair[]
 
     function fto_helper(sub_eq, α)
         δ = (gamma(α+1) * epsilon)^(1/α)
@@ -370,16 +370,29 @@ function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initial
         end
 
         new_eqs = Equation[]
-        def = Pair{Num, Int64}[]
-        sum = 0
-        for index in range(M, N-1; step=1)
+        def = Pair[]
+        if matrix
             new_z = Symbol(:ʐ, :_, i)
             i += 1
-            new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
-            new_eq = D(new_z) ~ sub_eq - γ_i(index)*new_z
+            γs = diagm([γ_i(index) for index in M:N-1])
+            cs = [c_i(index) for index in M:N-1]
+
+            new_z = only(@variables $new_z(iv)[1:N-M])
+            new_eq = D(new_z) ~ -γs*new_z .+ sub_eq
+            sum = dot(cs, new_z)
+            push!(def, new_z=>zeros(N-M))
             push!(new_eqs, new_eq)
-            push!(def, new_z=>0)
-            sum += c_i(index)*new_z
+        else
+            sum = 0
+            for index in range(M, N-1; step=1)
+                new_z = Symbol(:ʐ, :_, i)
+                i += 1
+                new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
+                new_eq = D(new_z) ~ sub_eq - γ_i(index)*new_z
+                push!(new_eqs, new_eq)
+                push!(def, new_z=>0)
+                sum += c_i(index)*new_z
+            end
         end
         return (new_eqs, def, sum)
     end

test/fractional_to_ordinary.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, ODEInterfaceDiffEq, SpecialFunctions
+using ModelingToolkit, OrdinaryDiffEq, ODEInterfaceDiffEq, SpecialFunctions, LinearAlgebra
 using Test
 
 # Testing for α < 1
@@ -55,31 +55,37 @@ while(time <= 1)
     time += 0.1
 end
 
-# # Testing for example 2 of Section 7
-# @independent_variables t
-# @variables x(t) y(t)
-# D = Differential(t)
-# tspan = (0., 220.)
+# Testing for example 2 of Section 7
+@independent_variables t
+@variables x(t) y(t)
+D = Differential(t)
+tspan = (0., 220.)
 
-# sys = fractional_to_ordinary([1 - 4*x + x^2 * y, 3*x - x^2 * y], [x, y], [1.3, 0.8], 10^-8, 220; initials=[[1.2, 1], 2.8]; matrix=true)
-# prob = ODEProblem(sys, [], tspan)
-# sol = solve(prob, radau5(), abstol = 1e-8, reltol = 1e-8)
+sys = fractional_to_ordinary([1 - 4*x + x^2 * y, 3*x - x^2 * y], [x, y], [1.3, 0.8], 10^-8, 220; initials=[[1.2, 1], 2.8]; matrix=true)
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-8, reltol = 1e-8)
 
-# @test isapprox(1.0097684171, sol(220, idxs=x), atol=1e-5)
-# @test isapprox(2.1581264031, sol(220, idxs=y), atol=1e-5)
+@test isapprox(1.0097684171, sol(220, idxs=x), atol=1e-5)
+@test isapprox(2.1581264031, sol(220, idxs=y), atol=1e-5)
 
-# #Testing for example 3 of Section 7
-# @independent_variables t
-# @variables x_0(t)
-# D = Differential(t)
-# tspan = (0., 5000.)
+#Testing for example 3 of Section 7
+@independent_variables t
+@variables x_0(t)
+D = Differential(t)
+tspan = (0., 5000.)
+
+function expect(t)
+    return sqrt(2) * sin(t + pi/4)
+end
+
+sys = linear_fractional_to_ordinary([3, 2.5, 2, 1, .5, 0], [1, 1, 1, 4, 1, 4], 6*cos(t), 10^-5, 5000; initials=[1, 1, -1])
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
 
-# function expect(t)
-#     return sqrt(2) * sin(t + pi/4)
-# end
+@test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-5)
 
-# sys = linear_fractional_to_ordinary([3, 2.5, 2, 1, .5, 0], [1, 1, 1, 4, 1, 4], 6*cos(t), 10^-5, 5000; initials=[1, 1, -1])
-# prob = ODEProblem(sys, [], tspan)
-# sol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
+msys = linear_fractional_to_ordinary([3, 2.5, 2, 1, .5, 0], [1, 1, 1, 4, 1, 4], 6*cos(t), 10^-5, 5000; initials=[1, 1, -1], matrix=true)
+mprob = ODEProblem(sys, [], tspan)
+msol = solve(prob, radau5(), abstol = 1e-5, reltol = 1e-5)
 
-# @test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-5)
+@test isapprox(expect(5000), msol(5000, idxs=x_0), atol=1e-5)