Combined matrix ver. with regular

Author: fchen121

src/systems/diffeqs/basic_transformations.jl

@@ -211,12 +211,15 @@ while(time <= 1)
 end
 ```
 """
-function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0, additional_eqs = [], iv = only(@independent_variables t))
+function fractional_to_ordinary(
+        eqs, variables, alphas, epsilon, T;
+        initials = 0, additional_eqs = [], iv = only(@independent_variables t), matrix=false
+)
     D = Differential(iv)
     i = 0
     all_eqs = Equation[]
-    all_def = Pair{Num, Int64}[]
-
+    all_def = Pair[]
+    
     function fto_helper(sub_eq, sub_var, α; initial=0)
         alpha_0 = α
 
@@ -248,38 +251,61 @@ function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0
         end
 
         new_eqs = Equation[]
-        def = Pair{Num, Int64}[]
+        def = Pair[]
 
-        if (α < 1)  
-            rhs = initial
-            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)
-                rhs += c_i(index)*new_z
+        if matrix
+            new_z = Symbol(:ʐ, :_, i)
+            i += 1
+            γs = diagm([γ_i(index) for index in M:N-1])
+            cs = [c_i(index) for index in M:N-1]
+
+            if (α < 1)
+                new_z = only(@variables $new_z(iv)[1:N-M])
+                new_eq = D(new_z) ~ -γs*new_z .+ sub_eq
+                rhs = dot(cs, new_z) + initial
+                push!(def, new_z=>zeros(N-M))
+            else
+                new_z = only(@variables $new_z(iv)[1:N-M, 1:m])
+                new_eq = D(new_z) ~ -γs*new_z + hcat(fill(sub_eq, N-M, 1), collect(new_z[:, 1:m-1]*diagm(1:m-1)))
+                rhs = coeff*sum(cs[i]*new_z[i, m] for i in 1:N-M)
+                for (index, value) in enumerate(initial)
+                    rhs += value * iv^(index - 1) / gamma(index)
+                end
+                push!(def, new_z=>zeros(N-M, m))
             end
+            push!(new_eqs, new_eq)
         else
-            rhs = 0
-            for (index, value) in enumerate(initial)
-                rhs += value * iv^(index - 1) / gamma(index)
-            end
-            for index in range(M, N-1; step=1)
-                new_z = Symbol(:ʐ, :_, i)
-                i += 1
-                γ = γ_i(index)
-                base = sub_eq
-                for k in range(1, m; step=1)
-                    new_z = Symbol(:ʐ, :_, index-M, :_, k)
+            if (α < 1)  
+                rhs = initial
+                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) ~ base - γ*new_z
-                    base = k * new_z
+                    new_eq = D(new_z) ~ sub_eq - γ_i(index)*new_z
                     push!(new_eqs, new_eq)
                     push!(def, new_z=>0)
+                    rhs += c_i(index)*new_z
+                end
+            else
+                rhs = 0
+                for (index, value) in enumerate(initial)
+                    rhs += value * iv^(index - 1) / gamma(index)
+                end
+                for index in range(M, N-1; step=1)
+                    new_z = Symbol(:ʐ, :_, i)
+                    i += 1
+                    γ = γ_i(index)
+                    base = sub_eq
+                    for k in range(1, m; step=1)
+                        new_z = Symbol(:ʐ, :_, index-M, :_, k)
+                        new_z = ModelingToolkit.unwrap(only(@variables $new_z(iv)))
+                        new_eq = D(new_z) ~ base - γ*new_z
+                        base = k * new_z
+                        push!(new_eqs, new_eq)
+                        push!(def, new_z=>0)
+                    end
+                    rhs += coeff*c_i(index)*new_z
                 end
-                rhs += coeff*c_i(index)*new_z
             end
         end
         push!(new_eqs, sub_var ~ rhs)
@@ -385,78 +411,6 @@ function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initial
     return mtkcompile(sys)
 end
 
-function fto_matrix(eqs, variables, alphas, epsilon, T; initials = 0, additional_eqs = [], iv = only(@independent_variables t))
-    D = Differential(iv)
-    i = 0
-    all_eqs = Equation[]
-    all_def = Pair[]
-
-    function fto_helper(sub_eq, sub_var, α; initial=0)
-        alpha_0 = α
-
-        if (α > 1)
-            coeff = 1/(α - 1)
-            m = 2
-            while (α - m > 0)
-                coeff /= α - m
-                m += 1
-            end
-            alpha_0 = α - m + 1
-        end
-
-        δ = (gamma(alpha_0+1) * epsilon)^(1/alpha_0)
-        a = pi/2*(1-(1-alpha_0)/((2-alpha_0) * log(epsilon^-1)))
-        h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(alpha_0 - 1)))
-
-        x_sub = (gamma(2-alpha_0) * epsilon)^(1/(1-alpha_0))
-        x_sup = -log(gamma(1-alpha_0) * epsilon)
-        M = floor(Int, log(x_sub / T) / h)
-        N = ceil(Int, log(x_sup / δ) / h)
-
-        function c_i(index)
-            h * sin(pi * alpha_0) / pi * exp((1-alpha_0)*h*index)
-        end
-
-        function γ_i(index)
-            exp(h * index)
-        end
-
-        new_eqs = Equation[]
-        def = Pair[]
-        new_z = Symbol(:ʐ, :_, i)
-        i += 1
-        γs = diagm([γ_i(index) for index in M:N-1])
-        cs = [c_i(index) for index in M:N-1]
-
-        if (α < 1)
-            new_z = only(@variables $new_z(iv)[1:N-M])
-            new_eq = D(new_z) ~ -γs*new_z .+ sub_eq
-            rhs = dot(cs, new_z) + initial
-            push!(def, new_z=>zeros(N-M))
-        else
-            new_z = only(@variables $new_z(iv)[1:N-M, 1:m])
-            new_eq = D(new_z) ~ -γs*new_z + hcat(fill(sub_eq, N-M, 1), collect(new_z[:, 1:m-1]*diagm(1:m-1)))
-            rhs = coeff*sum(cs[i]*new_z[i, m] for i in 1:N-M)
-            for (index, value) in enumerate(initial)
-                rhs += value * iv^(index - 1) / gamma(index)
-            end
-            push!(def, new_z=>zeros(N-M, m))
-        end
-        push!(new_eqs, new_eq)
-        push!(new_eqs, sub_var ~ rhs)
-        return (new_eqs, def)
-    end
-
-    for (eq, cur_var, alpha, init) in zip(eqs, variables, alphas, initials)
-        (new_eqs, def) = fto_helper(eq, cur_var, alpha; initial=init)
-        append!(all_eqs, new_eqs)
-        append!(all_def, def)
-    end
-    append!(all_eqs, additional_eqs)
-    @named sys = System(all_eqs, iv; defaults=all_def)
-    return mtkcompile(sys)
-end
-
 """
     change_independent_variable(
         sys::System, iv, eqs = [];

test/fractional_to_ordinary.jl

@@ -30,7 +30,7 @@ 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 = fto_matrix(eqs, x, α, 10^-7, 1)
+sys = fractional_to_ordinary(eqs, x, α, 10^-7, 1; matrix=true)
 
 prob = ODEProblem(sys, [], tspan)
 sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10)
@@ -55,31 +55,31 @@ 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 = fto_matrix([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])
-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
+# 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)
+# 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)
 
-@test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-5)
+# @test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-5)