How to convert an ODESystem to NonlinearSystem?

Specific Domains Modelling & Simulations
,
1

I know there is a “convert_system” function which can convert a “NonlinearSystem” to an “ODESystem”. Can we do it backwards? So that we can have both steady-state and dynamic simulation with a single code base.

2 Likes
2

Good point. Open an issue. Right now, if you have a NonlinearSystem you can use the SteadyStateDiffEq.jl methods to solve the nonlinear system via ODEs, but it’s not at the symbolic level.

1 Like
3

It turns out that a similar issue has been opened 4 weeks ago.
link to issue

4

Okay cool. Yeah we’ll get to it: it’s not too difficult. Also the NonlinearSolve.jl docs are what’s on my table right now, and there’s some simple fixes there that would make this straight forward. So check back there in a little bit, and remind me if I drop it.

5

I am also interested in converting an ODESystem to a NonlinearSystem. I wonder if this feature has been implemented or if there is a good workaround @ChrisRackauckas. I am aware of doing something like this for an ODESystem but then we don’t have the nice solution interface anymore

dfg = NonlinearProblem(odeprob)
nsol = solve(dfg, SciMLNLSolve.NLSolveJL())
1 Like
6

The solution interface is the same? MWE?

7

Here is an MWE showcasing what would be convenient but doesn’t work currently. It boils down to that for an ODEProblem, you can access the solution with sol[u(t,x,y)] while for a NonlinearProblem created from that ODEProblem you cannot access it via nsol[u(t,x,y)] or nsol[u(x,y)]. I could rewrite the equation to a nonlinear problem which would allow to access nsol[u(x,y)], but it would be nice if there was a convenience function that would do that conversion. I believe that is also what the original post was looking for.

using ModelingToolkit, OrdinaryDiffEq, MethodOfLines, Plots, SciMLNLSolve

@parameters x y t
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt  = Differential(t)

# Boundary conditions
bcs = [
	u(0,x,y) ~ 0.0,
	u(t,0,y) ~ 0.0, u(t,1,y) ~ 0.0,
	u(t,x,0) ~ 0.0, u(t,x,1) ~ 0.0
]

# Space and time domains
domains = [
	t ∈ IntervalDomain(0.0,1.0),
	x ∈ IntervalDomain(0.0,1.0),
	y ∈ IntervalDomain(0.0,1.0)
]

# Discretization
dx = 0.04
nd = round(Int, 1 / dx) + 1
discretization = MOLFiniteDifference([x=>dx,y=>dx],t)

eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + sin(pi*x)*sin(pi*y)
@named pde_system = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

prob = discretize(pde_system,discretization)

sol = solve(prob,TRBDF2(),saveat=0.1)

# solution interface accessible via u(t,x,y)
sol[u(t,x,y)]

dfg = NonlinearProblem(prob)
nsol = solve(dfg, SciMLNLSolve.NLSolveJL())

# solution interface not accessible via u(t,x,y): ArgumentError: u(t, x, y) is neither an observed nor a state variable.
nsol[u(t,x,y)] # does not work

# instead, we have to access the solution via the solution vector - this is not convenient for larger models
M_ss = zeros(nd,nd)
M_ss[2:end-1, 2:end-1] .= reshape(nsol.u, nd-2, nd-2)

MWE adapted from Control of PDE systems · JuliaSimControl (juliahub.com)

8

@ChrisRackauckas - you can find the MWE in the post above. And btw, big thanks to all your great work with SciML and MethodOfLines.jl in particular :clap:

9

Yeah that’s Unable to Index into a steadystate solution using symbols · Issue #338 · SciML/MethodOfLines.jl · GitHub, it’s on our mind and needs to get fixed.

1 Like