# Max function in differential equation

**URL:** <https://discourse.julialang.org/t/max-function-in-differential-equation/89176>\
**Category:** General Usage\
**Tags:** question, diffeq\
**Created:** [October 24, 2022, 5:02am UTC](https://discourse.julialang.org/t/max-function-in-differential-equation/89176 "2022-10-24T05:02:52Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ge0rges](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ge0rges/32/38357_2.png) [@Ge0rges](https://discourse.julialang.org/u/Ge0rges)\
**Post date:** [October 24, 2022, 5:02am UTC](https://discourse.julialang.org/t/max-function-in-differential-equation/89176/1 "2022-10-24T05:02:52Z")

</div>

Hello all,

I have a set of differential equations I’m trying to solve. One of them includes a `max()` function. I am wondering if that is the reason why I am sometimes unable to solve this equation. Can the solver “integrate” `max()`? Do I need to specifically handle that discontinuity somehow?

Here is the model code:

```julia
# Load state conditions
    pOC_content = u[1]
    dOC_content = u[2]
    inorganic_carbon_content = u[3]
    cell_count = u[4]

    ## CELL COUNT
    # Growth
    growth = mu_max * (dOC_content / (Ks + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

    # Organic carbon requirement
    required_dOC_per_cell = maintenance_per_cell
    required_dOC = required_dOC_per_cell * cell_count

    # Starvation Deaths
    dOC_missing = max(required_dOC - dOC_content, 0)
    starvation_deaths = dOC_missing == 0 ? 0 : dOC_missing / required_dOC_per_cell

    # Total Deaths
    deaths = starvation_deaths

    ## CARBON
    dOC_consumption = required_dOC_per_cell * (cell_count - deaths)
    fixed_carbon = inorganic_carbon_fixing_rate * inorganic_carbon_content

    # Particulate Organic carbon
    du[1] = pOC_input_rate - eea_rate*cell_count

    # Dissolved Organic carbon
    du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count + fixed_carbon

    # Inorganic carbon
    du[3] = inorganic_carbon_input_rate + inorganic_carbon_per_cell * (deaths - growth) + required_dOC - fixed_carbon

    # Net cell count change
    du[4] = growth - deaths

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 24, 2022, 12:27pm UTC](https://discourse.julialang.org/t/max-function-in-differential-equation/89176/2 "2022-10-24T12:27:07Z")

</div>

the discontinuity is a problem. check out the documentation for callbacks.

---

<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 24, 2022, 12:42pm UTC](https://discourse.julialang.org/t/max-function-in-differential-equation/89176/3 "2022-10-24T12:42:30Z")

</div>

> [@Handling Instability When Solving ODE Problems](https://discourse.julialang.org/t/handling-instability-when-solving-ode-problems/9019/5):
>
> Thanks for the example. I will share as little as possible but just a snippet in order to explain what’s going on to others. This is a good teaching opportunity. I was able to find out that the problem is that, in the statement k9\*X3/(k10\*k11/(k12\*(1.-min(1.,k13+k14\*X5))) + X3) when X3 and X5 are dual numbers, the value of k11/(k12\*(1.-min(1.,k13+k14\*X5))) is NaN for the derivative components when the condition makes it the constant 1.0 (since the derivative of the minimum of in the denominat…

---

<div class="post-metadata">

**Author:** ![Ge0rges](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ge0rges/32/38357_2.png) [@Ge0rges](https://discourse.julialang.org/u/Ge0rges)\
**Post date:** [October 24, 2022, 10:01pm UTC](https://discourse.julialang.org/t/max-function-in-differential-equation/89176/4 "2022-10-24T22:01:45Z")

</div>

Thanks for that. I started by trying a soft max approach as follows.

```julia
    # Define a soft max() function.
    function softmax(x, y, k)
        return log(exp(k*x) + exp(k*y))/k
    end
    dOC_missing = softmax(required_dOC - dOC_content, 0, 10) # max(required_dOC - dOC_content, 0)

```

However the yielded instability in previously stable solutions, even for larger or smaller `K` (1, 10, 100). I then opted to implement the callback solution, however this runs out of iterations only in this case.

```julia
    ...
    condition(u, t, integrator) = p[2] * u[4] - u[2]
    affect!(integrator) = nothing
    max_cb = ContinuousCallback(condition, affect!, save_positions=(false,false))

    # Callback list
    cbs = CallbackSet(carbon_add_cb, max_cb)

    # Build the ODE Problem and Solve
    sol = solve(prob, Rosenbrock23(), callback=cbs, tstops=stops, isoutofdomain=is_invalid_domain, maxiters=1e7)

```

Tried with `autodiff=False` and with `abstol=1e-10, reltol=1e-10` as well.

Given that my max is `max(X, 0)`, the condition is simply `X`. I believe the reason this is happening is because the values in `X` approach equality but not exactly, therefore `X` tends to 0 rather is equal to it. What is the right way to handle this? I have seen this throughout my model where some variables which are analytically equal to 0 in certain cases tend to 0 instead (e.g. 1e-324).

Reproducible callback failing:

```julia
function solve_model(p, u0)
    p = [2.767317774014066e-7, 0.7433220714328409, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0]
    u0 = [1.6382099999999998e13, 3.4061512877350002e10, 0.0, 100000.0, 1.461e7]

    # Callback for the max() function in the model - causes a discontinuity.
    # If max equation changes, this condition will have to change.
    condition(u, t, integrator) = p[2] * u[4] - u[2]
    affect!(integrator) = nothing
    max_cb = ContinuousCallback(condition, affect!, save_positions=(false,false))

    # Build the ODE Problem and Solve
    prob = ODEProblem(model, u0[1:end-1], (0.0, last(u0)), p)
    sol = solve(prob, Rosenbrock23(), callback=max_cb, tstops=stops, isoutofdomain=is_invalid_domain, maxiters=1e7)

    return sol
end

# Implements the differential equations that define the model
function model(du,u,p,t)
    # Paramaters
    mu_max = p[1]
    maintenance_per_cell = p[2]
    dOC_per_cell = p[3]

    carrying_capacity = p[4]
    pOC_input_rate = p[5]
    dOC_input_rate = p[6]
    inorganic_carbon_input_rate = p[7]
    inorganic_carbon_fixing_rate = p[8]
    inorganic_carbon_per_cell = p[9]
    eea_rate = p[10]
    Ks = p[11]

    # Load state conditions
    pOC_content = u[1]
    dOC_content = u[2]
    inorganic_carbon_content = u[3]
    cell_count = u[4]

    growth = mu_max * (dOC_content / (Ks + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

    required_dOC = maintenance_per_cell * cell_count

    dOC_missing = max(required_dOC - dOC_content, 0)
    starvation_deaths = dOC_missing == 0 ? 0 : dOC_missing / required_dOC_per_cell

    deaths = starvation_deaths

    dOC_consumption = maintenance_per_cell * (cell_count - deaths)
    
    eea_rate = eea_rate * (1 / (1 + exp(-pOC_content))) * 2 - eea_rate
    
    du[1] = pOC_input_rate - eea_rate*cell_count

    du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count

    du[4] = growth - deaths
end

```

Example of working parameters:

```julia
p = [0.06, 0.1806273414373398, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0]

```