Create test case for FDE to ODE

Author: fchen121

src/ModelingToolkit.jl

@@ -298,8 +298,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        tunable_parameters, isirreducible, getdescription, hasdescription,
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
-export liouville_transform, change_independent_variable, substitute_component,
-       add_accumulations, noise_to_brownians, Girsanov_transform, change_of_variables
+export liouville_transform, change_independent_variable, substitute_component, add_accumulations,
+       noise_to_brownians, Girsanov_transform, change_of_variables, fractional_to_ordinary
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation

src/systems/diffeqs/basic_transformations.jl

@@ -184,14 +184,14 @@ function change_of_variables(
     return new_sys
 end
 
-function FDE_to_ODE(eqs, variables, alphas, epsilon, T; initials=0)
+function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0)
     @independent_variables t
     D = Differential(t)
     i = 0
     all_eqs = Equation[]
     all_def = Pair{Num, Int64}[]
 
-    function FDE_helper(sub_eq, sub_var, α; initial=0)
+    function fto_helper(sub_eq, sub_var, α; initial=0)
         alpha_0 = α
         if (α > 1)
             coeff = 1/(α - 1)
@@ -260,7 +260,7 @@ function FDE_to_ODE(eqs, variables, alphas, epsilon, T; initials=0)
     end
 
     for (eq, cur_var, alpha, init) in zip(eqs, variables, alphas, initials)
-        (new_eqs, def) = FDE_helper(eq, cur_var, alpha; initial=init)
+        (new_eqs, def) = fto_helper(eq, cur_var, alpha; initial=init)
         append!(all_eqs, new_eqs)
         append!(all_def, def)
     end

test/fractional_to_ordinary.jl

@@ -0,0 +1,68 @@
+using ModelingToolkit, OrdinaryDiffEq, ODEInterfaceDiffEq, SpecialFunctions
+using Test
+
+# Testing for α < 1
+# Uses example 1 from Section 7 of https://arxiv.org/pdf/2506.04188
+@independent_variables t
+@variables x(t)
+D = Differential(t)
+tspan = (0., 1.)
+
+function expect(t, α)
+    return (3/2*t^(α/2) - t^4)^2
+end
+
+α = 0.5
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-3)
+    time += 0.1
+end
+
+α = 0.3
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-3)
+    time += 0.1
+end
+
+α = 0.9
+eqs = (9*gamma(1 + α)/4) - (3*t^(4 - α/2)*gamma(5 + α/2)/gamma(5 - α/2))
+eqs += (gamma(9)*t^(8 - α)/gamma(9 - α)) + (3/2*t^(α/2)-t^4)^3 - x^(3/2)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1)
+
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-10, reltol = 1e-10)
+
+time = 0
+while(time <= 1)
+    @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-3)
+    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.)
+
+sys = fractional_to_ordinary([1 - 4*x + x^2 * y, 3*x - x^2 * y], [x, y], [1.3, 0.8], 10^-6, 220; initials=[[1.2, 1], 2.8])
+prob = ODEProblem(sys, [], tspan)
+sol = solve(prob, radau5(), abstol = 1e-8, reltol = 1e-8)
+
+@test isapprox(1.0097684171, sol(220, idxs=x), atol=1e-3)
+@test isapprox(2.1581264031, sol(220, idxs=y), atol=1e-3)