0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ What do you mean by $k_t(s^-1)$? Does it perhaps mean that the dimensions of $k_t$ are $s^{-1}$? $\endgroup$ – Fabio Somenzi Nov 28 '18 at 18:18
  • $\begingroup$ It's some kind of rate so 'per second' $\endgroup$ – maria1991 Dec 1 '18 at 17:03
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks so much! I don't have access to MATLAB at the minute to give this ago but as soon as I do I'll let you know how I get on! $\endgroup$ – maria1991 Dec 1 '18 at 17:04
  • $\begingroup$ @maria1991 If the preceding answer has given you full satisfaction, please check the answer. $\endgroup$ – Jean Marie Mar 24 '19 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.