From f9f525704ba79f3b530cf3257d9dcf55956453fd Mon Sep 17 00:00:00 2001 From: fchen121 Date: Wed, 25 Jun 2025 09:43:49 -0400 Subject: [PATCH 1/9] Created FDE to ODE for one equation with alpha < 1 --- src/systems/diffeqs/basic_transformations.jl | 40 ++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 9312966528..69912dd007 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -184,6 +184,46 @@ function change_of_variables( return new_sys end +function FDE_to_ODE(eqs, alpha, epsilon, T; initial=0) + @independent_variables t + @variables x(t) + D = Differential(t) + δ = (gamma(alpha+1) * epsilon)^(1/alpha) + + a = pi/2*(1-(1-alpha)/((2-alpha) * log(epsilon^-1))) + h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(alpha - 1))) + + x_sub = (gamma(2-alpha) * epsilon)^(1/(1-alpha)) + x_sup = -log(gamma(1-alpha) * epsilon) + M = floor(Int, log(x_sub / T) / h) + N = ceil(Int, log(x_sup / δ) / h) + + function c_i(index) + h * sin(pi * alpha) / pi * exp((1-alpha)*h*index) + end + + function γ_i(index) + exp(h * index) + end + + new_eqs = Equation[] + rhs = initial + def = Pair{Num, Int64}[] + + for index in range(M, N-1; step=1) + new_z = Symbol(:z, :_, index-M) + new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) + new_eq = D(new_z) ~ eqs - γ_i(index)*new_z + push!(new_eqs, new_eq) + push!(def, new_z=>0) + rhs += c_i(index)*new_z + end + eq = x ~ rhs + push!(new_eqs, eq) + @mtkcompile sys = System(new_eqs, t; defaults=def) + return sys +end + """ change_independent_variable( sys::System, iv, eqs = []; From 6fd7d31bd7322cca2f56d963ba1e75ec09b979eb Mon Sep 17 00:00:00 2001 From: fchen121 Date: Wed, 25 Jun 2025 21:21:06 -0400 Subject: [PATCH 2/9] Implement FDE to ODE for single equation with alpha > 1 --- src/systems/diffeqs/basic_transformations.jl | 58 +++++++++++++++----- 1 file changed, 43 insertions(+), 15 deletions(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 69912dd007..d997ca549a 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -188,18 +188,29 @@ function FDE_to_ODE(eqs, alpha, epsilon, T; initial=0) @independent_variables t @variables x(t) D = Differential(t) - δ = (gamma(alpha+1) * epsilon)^(1/alpha) - a = pi/2*(1-(1-alpha)/((2-alpha) * log(epsilon^-1))) - h = 2*pi*a / log(1 + (2/epsilon * (cos(a))^(alpha - 1))) + alpha_0 = alpha + if (alpha > 1) + coeff = 1/(alpha - 1) + m = 2 + while (alpha - m > 0) + coeff /= alpha - m + m += 1 + end + alpha_0 = alpha - m + 1 + end - x_sub = (gamma(2-alpha) * epsilon)^(1/(1-alpha)) - x_sup = -log(gamma(1-alpha) * epsilon) + δ = (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) / pi * exp((1-alpha)*h*index) + h * sin(pi * alpha_0) / pi * exp((1-alpha_0)*h*index) end function γ_i(index) @@ -210,18 +221,35 @@ function FDE_to_ODE(eqs, alpha, epsilon, T; initial=0) rhs = initial def = Pair{Num, Int64}[] - for index in range(M, N-1; step=1) - new_z = Symbol(:z, :_, index-M) - new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) - new_eq = D(new_z) ~ eqs - γ_i(index)*new_z - push!(new_eqs, new_eq) - push!(def, new_z=>0) - rhs += c_i(index)*new_z + if (alpha < 1) + for index in range(M, N-1; step=1) + new_z = Symbol(:z, :_, index-M) + new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) + new_eq = D(new_z) ~ eqs - γ_i(index)*new_z + push!(new_eqs, new_eq) + push!(def, new_z=>0) + rhs += c_i(index)*new_z + end + else + for index in range(M, N-1; step=1) + new_z = Symbol(:z, :_, index-M) + start = eqs + γ = γ_i(index) + for k in range(1, m; step=1) + new_z = Symbol(:z, :_, index-M, :_, k) + new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) + new_eq = D(new_z) ~ start - γ*new_z + start = k * new_z + push!(new_eqs, new_eq) + push!(def, new_z=>0) + end + rhs += coeff*c_i(index)*new_z + end end eq = x ~ rhs push!(new_eqs, eq) - @mtkcompile sys = System(new_eqs, t; defaults=def) - return sys + @named sys = System(new_eqs, t; defaults=def) + return mtkcompile(sys) end """ From 8a1101d9fd8668d6e3484900e80660361b32e6a6 Mon Sep 17 00:00:00 2001 From: fchen121 Date: Wed, 25 Jun 2025 21:38:47 -0400 Subject: [PATCH 3/9] Update FDE_to_ODE to handle multiple FDEs --- src/systems/diffeqs/basic_transformations.jl | 112 ++++++++++--------- 1 file changed, 62 insertions(+), 50 deletions(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index d997ca549a..7344125c38 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -184,71 +184,83 @@ function change_of_variables( return new_sys end -function FDE_to_ODE(eqs, alpha, epsilon, T; initial=0) +function FDE_to_ODE(eqs, alphas, epsilon, T; initials=0) @independent_variables t @variables x(t) D = Differential(t) - - alpha_0 = alpha - if (alpha > 1) - coeff = 1/(alpha - 1) - m = 2 - while (alpha - m > 0) - coeff /= alpha - m - m += 1 + i = 0 + all_eqs = Equation[] + all_def = Pair{Num, Int64}[] + + function FDE_helper(sub_eq, α; initial=0) + alpha_0 = α + if (α > 1) + coeff = 1/(α - 1) + m = 2 + while (α - m > 0) + coeff /= α - m + m += 1 + end + alpha_0 = α - m + 1 end - alpha_0 = alpha - 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) + δ = (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))) - function c_i(index) - h * sin(pi * alpha_0) / pi * exp((1-alpha_0)*h*index) - end + 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 γ_i(index) - exp(h * index) - end + function c_i(index) + h * sin(pi * alpha_0) / pi * exp((1-alpha_0)*h*index) + end - new_eqs = Equation[] - rhs = initial - def = Pair{Num, Int64}[] - - if (alpha < 1) - for index in range(M, N-1; step=1) - new_z = Symbol(:z, :_, index-M) - new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) - new_eq = D(new_z) ~ eqs - γ_i(index)*new_z - push!(new_eqs, new_eq) - push!(def, new_z=>0) - rhs += c_i(index)*new_z + function γ_i(index) + exp(h * index) end - else - for index in range(M, N-1; step=1) - new_z = Symbol(:z, :_, index-M) - start = eqs - γ = γ_i(index) - for k in range(1, m; step=1) - new_z = Symbol(:z, :_, index-M, :_, k) + + new_eqs = Equation[] + rhs = initial + def = Pair{Num, Int64}[] + + if (α < 1) + 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) ~ start - γ*new_z - start = 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 + for index in range(M, N-1; step=1) + new_z = Symbol(:z, :_, i) + i += 1 + γ = γ_i(index) + for k in range(1, m; step=1) + new_z = Symbol(:z, :_, index-M, :_, k) + new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) + new_eq = D(new_z) ~ sub_eq - γ*new_z + sub_eq = 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 + push!(new_eqs, x ~ rhs) + return (new_eqs, def) + end + + for (eq, alpha) in zip(eqs, alphas) + (new_eqs, def) = FDE_helper(eq, alpha) + append!(all_eqs, new_eqs) + append!(all_def, def) end - eq = x ~ rhs - push!(new_eqs, eq) - @named sys = System(new_eqs, t; defaults=def) + @named sys = System(all_eqs, t; defaults=all_def) return mtkcompile(sys) end From 8d2123417724a42fa595c555a1d21385e2cb953f Mon Sep 17 00:00:00 2001 From: fchen121 Date: Thu, 26 Jun 2025 15:38:31 -0400 Subject: [PATCH 4/9] Fix bug in FDE to ODE for alpha > 1 --- src/systems/diffeqs/basic_transformations.jl | 24 ++++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 7344125c38..887b73f99d 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -184,15 +184,14 @@ function change_of_variables( return new_sys end -function FDE_to_ODE(eqs, alphas, epsilon, T; initials=0) +function FDE_to_ODE(eqs, variables, alphas, epsilon, T; initials=0) @independent_variables t - @variables x(t) D = Differential(t) i = 0 all_eqs = Equation[] all_def = Pair{Num, Int64}[] - function FDE_helper(sub_eq, α; initial=0) + function FDE_helper(sub_eq, sub_var, α; initial=0) alpha_0 = α if (α > 1) coeff = 1/(α - 1) @@ -222,10 +221,10 @@ function FDE_to_ODE(eqs, alphas, epsilon, T; initials=0) end new_eqs = Equation[] - rhs = initial def = Pair{Num, Int64}[] - if (α < 1) + if (α < 1) + rhs = initial for index in range(M, N-1; step=1) new_z = Symbol(:z, :_, i) i += 1 @@ -236,27 +235,32 @@ function FDE_to_ODE(eqs, alphas, epsilon, T; initials=0) rhs += c_i(index)*new_z end else + rhs = 0 + for (index, value) in enumerate(initial) + rhs += value * t^(index - 1) / gamma(index) + end for index in range(M, N-1; step=1) new_z = Symbol(:z, :_, i) i += 1 γ = γ_i(index) + base = sub_eq for k in range(1, m; step=1) new_z = Symbol(:z, :_, index-M, :_, k) new_z = ModelingToolkit.unwrap(only(@variables $new_z(t))) - new_eq = D(new_z) ~ sub_eq - γ*new_z - sub_eq = k * new_z + 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 end - push!(new_eqs, x ~ rhs) + push!(new_eqs, sub_var ~ rhs) return (new_eqs, def) end - for (eq, alpha) in zip(eqs, alphas) - (new_eqs, def) = FDE_helper(eq, alpha) + for (eq, cur_var, alpha, init) in zip(eqs, variables, alphas, initials) + (new_eqs, def) = FDE_helper(eq, cur_var, alpha; initial=init) append!(all_eqs, new_eqs) append!(all_def, def) end From b936c8b64af0cdeaf98b2f8a94d693de542e71dc Mon Sep 17 00:00:00 2001 From: fchen121 Date: Thu, 26 Jun 2025 21:31:13 -0400 Subject: [PATCH 5/9] Create test case for FDE to ODE --- src/ModelingToolkit.jl | 4 +- src/systems/diffeqs/basic_transformations.jl | 6 +- test/fractional_to_ordinary.jl | 68 ++++++++++++++++++++ 3 files changed, 73 insertions(+), 5 deletions(-) create mode 100644 test/fractional_to_ordinary.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index d23d173b9a..4e9be46555 100644 --- a/src/ModelingToolkit.jl +++ b/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 diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 887b73f99d..648b3d77aa 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/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 diff --git a/test/fractional_to_ordinary.jl b/test/fractional_to_ordinary.jl new file mode 100644 index 0000000000..b2f65a1d8b --- /dev/null +++ b/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) \ No newline at end of file From 90534fab01a517d65134b5b157a6f12e919c8d97 Mon Sep 17 00:00:00 2001 From: fchen121 Date: Thu, 26 Jun 2025 21:45:11 -0400 Subject: [PATCH 6/9] Add docstring to fractional_to_ordinary --- src/systems/diffeqs/basic_transformations.jl | 26 ++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 648b3d77aa..390c605d82 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -184,6 +184,32 @@ function change_of_variables( return new_sys end +""" +Generates the system of ODEs to find solution to FDEs. + +Example: + +```julia +@independent_variables t +@variables x(t) +D = Differential(t) +tspan = (0., 1.) + +α = 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((3/2*time^(α/2) - time^4)^2, sol(time, idxs=x), atol=1e-3) + time += 0.1 +end +``` +""" function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0) @independent_variables t D = Differential(t) From 17f978be492aacd195b0dd3c3190d08faad4c48e Mon Sep 17 00:00:00 2001 From: fchen121 Date: Thu, 3 Jul 2025 11:08:04 -0400 Subject: [PATCH 7/9] Add new function to deal with multiple term problem --- Project.toml | 2 + src/ModelingToolkit.jl | 3 +- src/systems/diffeqs/basic_transformations.jl | 69 ++++++++++++++++++++ 3 files changed, 73 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 863277711a..4063a1f216 100644 --- a/Project.toml +++ b/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" diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 4e9be46555..81ed7d4272 100644 --- a/src/ModelingToolkit.jl +++ b/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 diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 390c605d82..7f004b799a 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/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 = []; From 4efd06f79ec152b65d538ea62f92099d626db378 Mon Sep 17 00:00:00 2001 From: fchen121 Date: Tue, 8 Jul 2025 22:19:02 -0400 Subject: [PATCH 8/9] Improve test cases --- test/fractional_to_ordinary.jl | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/test/fractional_to_ordinary.jl b/test/fractional_to_ordinary.jl index b2f65a1d8b..23a0e39595 100644 --- a/test/fractional_to_ordinary.jl +++ b/test/fractional_to_ordinary.jl @@ -7,6 +7,7 @@ using Test @variables x(t) D = Differential(t) tspan = (0., 1.) +timepoint = [0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.] function expect(t, α) return (3/2*t^(α/2) - t^4)^2 @@ -18,7 +19,7 @@ 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) +sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10) time = 0 while(time <= 1) @@ -32,11 +33,11 @@ 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) +sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10) time = 0 while(time <= 1) - @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-3) + @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-7) time += 0.1 end @@ -46,11 +47,11 @@ 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) +sol = solve(prob, radau5(), saveat=timepoint, abstol = 1e-10, reltol = 1e-10) time = 0 while(time <= 1) - @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-3) + @test isapprox(expect(time, α), sol(time, idxs=x), atol=1e-7) time += 0.1 end @@ -60,9 +61,25 @@ end 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]) +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]) 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) \ No newline at end of file +@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.) + +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) + +@test isapprox(expect(5000), sol(5000, idxs=x_0), atol=1e-6) \ No newline at end of file From 4f0753cc10f3d42a15f79f568e6c7c0de1c2cfa3 Mon Sep 17 00:00:00 2001 From: fchen121 Date: Thu, 10 Jul 2025 10:58:28 -0400 Subject: [PATCH 9/9] Update docstring --- src/systems/diffeqs/basic_transformations.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 7f004b799a..10ac80299e 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -294,6 +294,26 @@ function fractional_to_ordinary(eqs, variables, alphas, epsilon, T; initials = 0 return mtkcompile(sys) end +""" +Generates the system of ODEs to find solution to FDEs. + +Example: + +```julia +@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 linear_fractional_to_ordinary(degrees, coeffs, rhs, epsilon, T; initials = 0) @independent_variables t @variables x_0(t)