Add new function to deal with multiple term problem

Author: fchen121

Project.toml

@@ -45,6 +45,7 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -142,6 +143,7 @@ OrdinaryDiffEq = "6.82.0"
 OrdinaryDiffEqCore = "1.15.0"
 OrdinaryDiffEqDefault = "1.2"
 OrdinaryDiffEqNonlinearSolve = "1.5.0"
+Plots = "1.40.15"
 PrecompileTools = "1"
 Pyomo = "0.1.0"
 REPL = "1"

src/ModelingToolkit.jl

@@ -299,7 +299,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        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, fractional_to_ordinary
+       noise_to_brownians, Girsanov_transform, change_of_variables,
+       fractional_to_ordinary, linear_fractional_to_ordinary
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation

src/systems/diffeqs/basic_transformations.jl

@@ -294,6 +294,75 @@ function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0
     return mtkcompile(sys)
 end
 
+function linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initials = 0)
+    @independent_variables t
+    @variables x_0(t)
+    D = Differential(t)
+    i = 0
+    all_eqs = Equation[]
+    all_def = Pair{Num, Int64}[]
+
+    function fto_helper(sub_eq, α)
+        δ = (gamma(α+1) * epsilon)^(1/α)
+        a = pi/2*(1-(1-α)/((2-α) * log(epsilon^-1)))
+        h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(α - 1)))
+
+        x_sub = (gamma(2-α) * epsilon)^(1/(1-α))
+        x_sup = -log(gamma(1-α) * epsilon)
+        M = floor(Int, log(x_sub / T) / h)
+        N = ceil(Int, log(x_sup / δ) / h)
+
+        function c_i(index)
+            h * sin(pi * α) / pi * exp((1-α)*h*index)
+        end
+
+        function γ_i(index)
+            exp(h * index)
+        end
+
+        new_eqs = Equation[]
+        def = Pair{Num, Int64}[]
+        sum = 0
+        for index in range(M, N-1; step=1)
+            new_z = Symbol(:z, :_, i)
+            i += 1
+            new_z = ModelingToolkit.unwrap(only(@variables $new_z(t)))
+            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
+        return (new_eqs, def, sum)
+    end
+
+    previous = x_0
+    for i in range(1, ceil(Int, degrees[1]); step=1)
+        new_x = Symbol(:x, :_, i)
+        new_x = ModelingToolkit.unwrap(only(@variables $new_x(t)))
+        push!(all_eqs, D(previous) ~ new_x)
+        push!(all_def, previous => initials[i])
+        previous = new_x
+    end
+
+    new_rhs = -rhs
+    for (degree, coeff) in zip(degrees, coeffs)
+        rounded = ceil(Int, degree)
+        new_x = Symbol(:x, :_, rounded)
+        new_x = ModelingToolkit.unwrap(only(@variables $new_x(t)))
+        if isinteger(degree)
+            new_rhs += coeff * new_x
+        else
+            (new_eqs, def, sum) = fto_helper(new_x, rounded - degree)
+            append!(all_eqs, new_eqs)
+            append!(all_def, def)
+            new_rhs += coeff * sum
+        end
+    end
+    push!(all_eqs, 0 ~ new_rhs)
+    @named sys = System(all_eqs, t; defaults=all_def)
+    return mtkcompile(sys)
+end
+
 """
     change_independent_variable(
         sys::System, iv, eqs = [];