Skip to content

DE Transformation (Change of Variables) #3695

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 24 commits into from
Jun 12, 2025
Merged

Conversation

fchen121
Copy link
Contributor

@fchen121 fchen121 commented Jun 4, 2025

Created change of variable functions for ODE and SDE (#141) using #1073 as a basis.
I commented out the "Linear transformation to diagonal system" test as I couldn't figure out how to fix it.

@test isapprox(sol[x][end], new_sol[x][end], rtol=1e-4)


# # Linear transformation to diagonal system
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

u' = Au
A = P^{-1}DP is the eigendecomposition (D is the diagonal matrix of eigenvalues and P is the matrix of eigenvectors)
z = Pu
u' = P^{-1}DPu
z' = Dz

Comment on lines 101 to 103
if !iscomplete(sys)
sys = mtkcompile(sys)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

probably best to apply this before simplification so you already have the brownians?

Comment on lines 113 to 119
if neqs !== nothing
isSDE = true
neqs = [neqs[i,:] for i in 1:size(neqs,1)]

brownvars = map([Symbol(:B, :_, i) for i in 1:length(neqs[1])]) do name
unwrap(only(@brownian $name))
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AayushSabharwal if the system isn't simplified, what's a quick function to find the brownians?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's the other way around - any system with noise_eqs won't have brownians. If a system isn't simplified and has brownian terms in the equations, brownians(sys) contains the list of brownian variables.

Copy link
Contributor Author

@fchen121 fchen121 Jun 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a way to isolate the diffusion coefficients of the brownians before simplifying?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Manually, with Symbolics.linear_expansion. See

Is = Int[]
Js = Int[]
vals = Num[]
new_eqs = copy(eqs)
dvar2eq = Dict{Any, Int}()
for (v, dv) in enumerate(var_to_diff)
dv === nothing && continue
deqs = 𝑑neighbors(graph, dv)
if length(deqs) != 1
error("$(eqs[deqs]) is not handled.")
end
dvar2eq[fullvars[dv]] = only(deqs)
end
for (j, bj) in enumerate(brown_vars), i in 𝑑neighbors(graph, bj)
push!(Is, i)
push!(Js, j)
eq = new_eqs[i]
brown = fullvars[bj]
(coeff, residual, islinear) = Symbolics.linear_expansion(eq, brown)
islinear || error("$brown isn't linear in $eq")
new_eqs[i] = 0 ~ residual
push!(vals, coeff)
end
g = Matrix(sparse(Is, Js, vals))
sys = state.sys
@set! sys.eqs = new_eqs
@set! sys.unknowns = [v
for (i, v) in enumerate(fullvars)
if !iszero(new_idxs[i]) &&
invview(var_to_diff)[i] === nothing]
ode_sys = mtkcompile(
sys; simplify, inputs, outputs, disturbance_inputs, kwargs...)
eqs = equations(ode_sys)
sorted_g_rows = zeros(Num, length(eqs), size(g, 2))
for (i, eq) in enumerate(eqs)
dvar = eq.lhs
# differential equations always precede algebraic equations
_iszero(dvar) && break
g_row = get(dvar2eq, dvar, 0)
iszero(g_row) && error("$dvar isn't handled.")
g_row > size(g, 1) && continue
@views copyto!(sorted_g_rows[i, :], g[g_row, :])
end
# Fix for https://github.com/SciML/ModelingToolkit.jl/issues/2490
if sorted_g_rows isa AbstractMatrix && size(sorted_g_rows, 2) == 1
# If there's only one brownian variable referenced across all the equations,
# we get a Nx1 matrix of noise equations, which is a special case known as scalar noise
noise_eqs = reshape(sorted_g_rows[:, 1], (:, 1))
is_scalar_noise = true
elseif __num_isdiag_noise(sorted_g_rows)
# If each column of the noise matrix has either 0 or 1 non-zero entry, then this is "diagonal noise".
# In this case, the solver just takes a vector column of equations and it interprets that to
# mean that each noise process is independent
noise_eqs = __get_num_diag_noise(sorted_g_rows)
is_scalar_noise = false
else
noise_eqs = sorted_g_rows
is_scalar_noise = false
end
.

@ChrisRackauckas ChrisRackauckas merged commit 9d0b7ef into SciML:master Jun 12, 2025
39 of 49 checks passed
@ChrisRackauckas
Copy link
Member

@AayushSabharwal what's the right place in the new docs to add this?

@AayushSabharwal
Copy link
Member

The docstring should be added to the ## Model transformations section of model_building.md. Beyond that, an example? We could have a page demonstrating all the different model transformations.

@ChrisRackauckas
Copy link
Member

@fchen121 can you follow up with that PR? Then today we can talk about going towards Fractional Differential Equations

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants