Help solving very complex first order ODEs using ode45 - MATLAB - movement of water

Question

I am trying to solve this complex set of first order ODEs using ode45 on MATLAB. They describe the change in volume of water in the inner layer ($X_I$) and outer layer ($X_O$) of the Venus Flytrap's upper leaf (trap). This movement of water is what allows the trap to snap shut. I am looking to solve the ODEs for time for different parts of the prey capture process (which will have different conditions - eg. for the closing process, where the trap moves from an open to semi-closed state, we have $u_a$ and $u_c$ equal to zero.) These u terms are water transport rates driven by various gradients, $\alpha$ is the water supply rate from the roots and $\mu$ is the water consumption rate.

$dX_O/dt = \frac{\alpha X_O^2}{X_O^2 + X_I^2} + u_h + u_a - u_c -\mu X_O$

$dX_I/dt = \frac{\alpha X_I^2}{X_O^2 + X_I^2} - u_h - u_a + u_c -\mu X_I $

Where:

$X_O + X_I = 1$

$k_t(s^-1) = 10$

$k_f(s^-1) = 2$

$k_d(s^-1) = 0.000045$

$\alpha (h^-1) = 1$

$\mu (h^-1) = 1$

and

$u_h = k_t max\{X_I(t) - X_O(t),0\}\delta _t(t)$

$u_a =\left\{ \begin{array}{@{}ll@{}} k_f X_I(t) \delta _t(t), & \text{if}\ C \geq 14 \\ 0, & \text{if}\ C<14 \end{array}\right.$

where $\delta_ t(t) = \left\{ \begin{array}{@{}ll@{}} 1, & \text{if}\ 0 \leq t \leq 0.3 \\ 0, & \text{if}\ t>0.3 \end{array}\right.$

$u_c(t) = k_d \delta _c(t)$

where $\delta _c(t) =\left\{ \begin{array}{@{}ll@{}} 1, & \text{if}\ T_{start} \leq t \leq T_{start}+T_D \ (T_D=15 \ hours) \\ 0, & \text{otherwise} \end{array}\right.$

So far I have attempted to create a script file which defines the variables:

a = 3600;
mu = 3600;
Xo = 1 - Xi;

if (0 <= t1) && (t1 <= 0.3)
    Dt1 = 1;
else
    Dt1 = 0;
end

Uh = 10*max((Xo - Xi),0)*d;

if C >= 14
Ua = 2*Xi*Dt1; else Ua = 0; end

if (0 <= t2) && (t2 <= 54000)
    Dc = 1;
else 
    Dc = 0;
end

Uc = 0.000045*Dc;  

where a represents $\alpha$ and Dt and Dc represent $\delta _t$ and $\delta _c$ respectively. Note that I converted the values of $\alpha$ and $\mu$ into seconds. I created two times, t1 and t2 so that I could put $T_{start}$ as $0$ for $\delta _c(t)$.

If anyone could help me I would be extremely grateful, I've spent about 15 hours on this now and am still getting nowhere, despite my numerous searches online for hints.

Answer 1

If $X_I + X_O = 1$, then $\frac{d}{dt}(X_I+X_O) = 0$, which implies $\alpha = \mu$. We can also simplify the problem by letting ode45 compute, say, $X_O(t)$ and then computing $X_I(t) = 1 - X_O(t)$.

MATLAB's ode45 and its siblings are passed a handle to a function of $t$ and $X_O$ that computes $\frac{dX_O}{dt}$. That function may be built in two stages. First we write a function with all the parameters as inputs (besides $t$ and $X_O$). It may look like this:

function xdot = dxodt(t, x, kt, kf, kd, alpha, mu, C, Tstart, TD)
  c = alpha * x^2 / (2*x^2 - 2*x + 1);
  deltat = 1.0 * (0 <= t && t <= 0.3);
  uh = kt*max(1-2*x,0)*deltat;
  if C >= 14
    ua = kf * (1-x) * deltat;
  else
    ua = 0;
  end
  deltac = 1.0 * (Tstart <= t && t <= Tstart + TD);
  uc = kd * deltac;
  xdot = c + uh + ua - uc - mu*x;
end

We then fix the values of the parameters and invoke the solver by some code like this:

function venusflytrap
  %VENUSFLYTRAP simulate Venus Flytrap

  kt = 10;       % 1/s
  kf = 2;        % 1/s
  kd = 0.000045; % 1/s
  alpha = 3600;  % 1/s
  mu = alpha;
  Tstart = 0;
  TD = 15*3600;
  % These values are uneducated guesses.
  C = 1;
  xinit = 0;
  Tfin = 1;

  options = odeset('stats','on');
  fh = @(t,x) dxodt(t, x, kt, kf, kd, alpha, mu, C, Tstart, TD);
  sol = ode45(fh, [0, Tfin], xinit, options);
  plot(sol.x,sol.y)
end

If these two functions are placed in the same file, venusflytrap should come first, and the file should be named venusflytrap.m.