tanh gives value larger than 1

I'm trying to solve a set of ode with ode45. Parameter alpha_a and alpha_m are defined through a tanh function. Plotting this function gives values way larger than 1 (increasing ad infinitum), while the idea is to set a two-phase parameter that is between 0 and 1 (to estimate growth factors).

The main code is the following:

Y0 = [600; 10; 10; 0.01; 0; 0; 0];
[ts, ys] = ode45(@Pinto1, tspan, Y0); % solves the set of ode's with a general ode solver

function dBdt = Pinto1(t,B)
S = B(1); % substrate concentration
xa = B(2); % anodophilic / electrogens concentration
xm = B(3); % methanogens concentration
Mox = B(4); % oxidised mediator concentration
Q_hold = B(5); % space declaration for Q
alpha_a_hold = B(6);
alpha_m_hold= B(7);

% here are the formulas for S, xa, xm, Mox, Q

alpha_a = (1+tanh(Kx*(xa+xm-Xma)))/2; % biofilm retention fraction for electrogens
alpha_m = (1+tanh(Kx*(xa+xm-Xmm)))/2; % biofilm retention fraction for methanogens
dxadt = mua*xa - Kda*xa - alpha_a*D*xa; % electrogens concentration difference over time
dxmdt = mum*xm - Kdm*xm - alpha_m*D*xm; % methanogens concentration difference over time

% mua, Kda, mum, Kdm, and D are constants or variables defined unrelated to alpha_a and alpha_m.
dBdt = [dSdt; dxadt; dxmdt; dMoxdt; Q; alpha_a; alpha_m];
end

The calculation of S, Mox, and Q are unrelated to alpha_a and alpha_m.

My output vector contains alpha_a and alpha_m values of well over 3 within the part that I'm plotting.

As far as I know, tanh should have output values between -1 and 1.

Am I declaring my input and output wrong? Or is tanh not suited for use in an ode system? I cannot recreate the issue with just constants (leaving out xa and xm)

Answer 1

[Bjorn Gustavsson]

Maybe I'm getting things completely wrong, but the way I interpret your question you've plotted the last 2 components out of B returned bythe ODE-solver, and they grow in ways that "troubles" you?

If that's the case, then you've defined ODE-function such that the last two components of dBdt are alpha_a; alpha_m - meaning that the ode-solver returns the solution of the (components of the coupled ODE's) for the equations:

dF_1dt = alpha_a;
dF_2dt = alpha_m;

With alpha_a and alpha_m defined to be between 0 and 1 - this guarantees that F_1 and F_2 will grow for all times - maybe what you're interested in is the gradients of the F_1 and F_2 (last and second last component in your B) with respect to time - those would be your alpha_a and alpha_m.

HTH

[Bas de Wolf]

Hi Bjorn,

Thank you for your reply!

I plotted the gradient of my alpha values and plotted these: looks better now! I also messed this up for Q.

I think I misunderstood what the ode45 function does. I think you're saying that the output of the d/dt functions of the ode system is not the same as the 'input' for the next time step. So this would mean I'm getting a derivative of the actual values that I'm looking at, if I understand correctly.

Thank you for clarifying! I might have to take a look at my old uni books for ode systems...

[Bjorn Gustavsson]

Sounds about right. The ode45 et al functions integrates your ODE-function. That is they return the y(t) that solves:, and the function f(t,y) has to return ...

Have a look at your uni-litterature - and also the ODE-function documentation and examples...

All the best,