DE Transformation (Change of Variables)

[fchen121]

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.

[ChrisRackauckas]

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

[ChrisRackauckas]

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

[ChrisRackauckas]

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

[AayushSabharwal]

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.

[fchen121]

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

[AayushSabharwal]

Manually, with Symbolics.linear_expansion. See

ModelingToolkit.jl/src/systems/systems.jl

Lines 98 to 156 in faf322f

	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

.