# Max function in differential equation

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:

``````# 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
``````

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

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

``````    # 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.

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

# Callback list

# 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:

``````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]

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:

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