Csdrhrt

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.

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.