# DifferentialEquations: Enforce that variables are inside a domain (isoutofdomain, PositiveDomain)

**URL:** <https://discourse.julialang.org/t/differentialequations-enforce-that-variables-are-inside-a-domain-isoutofdomain-positivedomain/30452>\
**Category:** General Usage\
**Tags:** diffeq, sde\
**Created:** [October 29, 2019, 4:11pm UTC](https://discourse.julialang.org/t/differentialequations-enforce-that-variables-are-inside-a-domain-isoutofdomain-positivedomain/30452 "2019-10-29T16:11:10Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![LolianSh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/loliansh/32/10834_2.png) [@LolianSh](https://discourse.julialang.org/u/LolianSh)\
**Post date:** [October 29, 2019, 4:11pm UTC](https://discourse.julialang.org/t/differentialequations-enforce-that-variables-are-inside-a-domain-isoutofdomain-positivedomain/30452/1 "2019-10-29T16:11:11Z")

</div>

Hi,  
I am using DifferentialEquations.jl to solve a system of SDEs. This is a model of an electrochemical system, where I have reaction-diffusion equations.

My variables are the following (\theta, CO, OH, S\_t). All of these are extended variables, so they are arrays, which I concatenate to form the varibles array `u`.

So far I am using the `EM()` solver with a fixed time-step (`dt = dt_stable`).

During the integration I want to ensure that my variables do not go out of some physical bound (0\<\theta, CO, OH \< 1, S\_{t, max} \< S\_t \< S\_{t, min}).

I saw on the online documentation that I could use this solution:

`isoutofdomain = (u, p, t) -> (any(x -> x < 0, p.view(u, \theta)) || any(x -> x > 1, p.view(u, \theta)) || any(x -> x < p.St_min, p.view(u, St)))`

or using the `GeneralDomain` callback (whihc I still did not implement).  
(`p.view` is a `NamedTupleShape` object that I use to unpack `u` )

Reading on the shortcomings of these methods I saw that they simply reject the proposed integration step and adjust the time stepping in order to not go out of bounds in the first place.  
This is certainly nice but I am not sure whether the extrapolation of the solution makes sense when dealing with noise in my system.  
Additionally, I want to avoid the overhead of having to check my solution at each time step and repeat the integration a number of times to ensure it.

I was thinking of simply clamping the array `u` after the integration step, in order to ensure my bounds. This way I would also avoid defining the differential equations outside the domain of my variables. Do you think this would be a reasonable solution ?

So I would like to ask how I can manually change the state of `integrator.u` right after it has been updated, before going to the next integration step and before saving the result to the output object (either the output of the `solve()` function or by means of a `SavingCallback`).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 29, 2019, 4:20pm UTC](https://discourse.julialang.org/t/differentialequations-enforce-that-variables-are-inside-a-domain-isoutofdomain-positivedomain/30452/2 "2019-10-29T16:20:43Z")

</div>

> [@LolianSh](#):
>
> I was thinking of simply clamping the array `u` after the integration step, in order to ensure my bounds. This way I would also avoid defining the differential equations outside the domain of my variables. Do you think this would be a reasonable solution ?

For a fixed time step that is reasonable. All of the other strategies are made for adaptive integration.

> [@LolianSh](#):
>
> So I would like to ask how I can manually change the state of `integrator.u` right after it has been updated, before going to the next integration step and before saving the result to the output object

`DiscreteCallback`. Make `condition(u,t,integrator) = true` so it happens after every step, and then `affect!(integ) = (integ.u[...] = ...)`. You may then want to set `save_positions = (false,false)` in the construction of the callback so it doesn’t add extra save points.