# Solving different parts of state vector with different schemes/solvers

**URL:** <https://discourse.julialang.org/t/solving-different-parts-of-state-vector-with-different-schemes-solvers/99305>\
**Category:** Modelling & Simulations\
**Created:** [May 23, 2023, 9:33pm UTC](https://discourse.julialang.org/t/solving-different-parts-of-state-vector-with-different-schemes-solvers/99305 "2023-05-23T21:33:13Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![p\_f](https://avatars.discourse-cdn.com/v4/letter/p/45deac/32.png) [@p\_f](https://discourse.julialang.org/u/p_f)\
**Post date:** [May 23, 2023, 9:33pm UTC](https://discourse.julialang.org/t/solving-different-parts-of-state-vector-with-different-schemes-solvers/99305/1 "2023-05-23T21:33:13Z")

</div>

Hi all, I hope this question finds you well.

I have two state vectors `X1` and `X2` that comprise my total state `X = [X1; X2]`.

X1 evolves according to `dX1 = f(X1,X2) ` and X2 evolves according to `dX2 = g(X2) - h(X1,X2,dX1)` (note dependence on `dX1`). I have a specific timestepping scheme I’d like to use for `X2` where I treat `g` explicitly and `h` implicitly, and `h` is linear in X2 so I have a function for X2(t+dt) ready to go in terms of X2(t).

But for X1, the function `f(X1,X2)` is fairly complex so I like using the Julia solvers to adaptively timestep, do local error analysis etc.

I’m currently solving the system by passing `dX1 = f(X1,X2)` and `dX2 = 0` to a julia solver, and then using a DiscreteCallback to timestep X2 manually after every timestep. However, I assume the Julia solvers will find this approach weird (since f is returning 0 for dX2) and will result in poor solving.

I’m running a test with `dtmax = 1e-5` to try and ensure that X2 is being solved well, but now I’m losing most of the benefits of the lovely prewritten adaptive solvers that are so good for X1. That might work just fine, but the current ETA is about 7 hours so I thought I’d see if there’s an obvious better way of handling this situation.

Before trying this, I was just passing the rhs function `dX = [dX1; dX2]` to julia solvers and using a stiff solver, but I found myself running into trouble when moving certain parameters too high. E.g. julia solvers needing increasingly small (\< 10^-16) timesteps, I believe for X2. I would see X2 occasionally go negative which it shouldn’t, but the IMEX scheme I’ve implemented for X2(t+1) is positive by definition (is this iterative implicit solvers vs rearranged for X2(t+1)?).

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [May 24, 2023, 12:06am UTC](https://discourse.julialang.org/t/solving-different-parts-of-state-vector-with-different-schemes-solvers/99305/2 "2023-05-24T00:06:20Z")

</div>

Can you use a [`SplitODEProblem`](https://docs.sciml.ai/DiffEqDocs/latest/solvers/split_ode_solve/)?

For maintaining a positivity constraint, I have used `isoutofdomain` (under MIscelaneous [here](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/)) to reduce the timestep and try again when part of the state goes out of bounds.

You might also try [Manifold Projection](https://docs.sciml.ai/DiffEqDocs/latest/examples/kepler_problem/#Manifold-Projection) with a residual of `min(0, positive_state)`, although I have not tried that.