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

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.

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

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.

• 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! – maria1991 Dec 1 '18 at 17:04
• @maria1991 If the preceding answer has given you full satisfaction, please check the answer. – Jean Marie Mar 24 '19 at 20:41