MATLAB: Different input per time step ODE45

different inputode45timesteps

Hi
I'm trying to describe the waterflow in soil. The change in soil water content (dx/dt) is calculated by an ODE45 for three different depths. With only one input u (for example constant irrigation) it works. But irrigation is not constant, so u is different for every t. This is given in a vector u. t = 0:1:150; u is a vector of 150*1. but this keeps giving me the error:
Error using odearguments (line 92) @(TIME,X)FFLOW(TIME,X,U,P) returns a vector of length 152, but the length of initial conditions vector is 3. The vector returned by @(TIME,X)FFLOW(TIME,X,U,P) and the initial conditions vector must have the same number of elements. Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin); Error in Main_program (line 37) [time, xt] = ode45(fh, time, x0);
So my question is: How can I give a different input in my ODE function for every different t?
My main code is this:
p = [25;25;25];           %(constant) depth of soil layers [cm]
% If time is not imported
time = (0:1:150)';     % time vector
x0 = [0.058;0.042;0.046];       % vector for initial water content values   
% Calculate water flow
fh = @(time,x)fflow(time, x, u, p);
[time, xt] = ode45(fh, time, x0);
fflow, where the u is eventually used looks like this(K and diff are just different constant values):
qin = u;  % Volume flux density as input (=irr-ET)
q12 = -K(1)*((-abs(diff1)/p(1))+1);       % Volume flux density layer 1 to 2 [cm3 cm-2 h-1]
q23 = -K(2)*((-abs(diff2)/p(2))+1);       % Volume flux density layer 2 to 3 [cm3 cm-2 h-1]
q34 = -K(3)*((-abs(diff3)/p(3))+1);       % Volume flux density layer 3 to 4 [cm3 cm-2 h-1]
% Compute change in water content theta (x)
dx1dt = (qin-q12)/p(1);       
dx2dt = (q12/p(1)-q23/p(2));
dx3dt = (q23/p(2)-q34/p(3));
dxdt = [dx1dt; dx2dt; dx3dt];
If more information is needed, please ask. I hope you can help.

Best Answer

function main
p = [25;25;25];           %(constant) depth of soil layers [cm]
% If time is not imported
time = (0:1:150)';     % time vector
x0 = [0.058;0.042;0.046];       % vector for initial water content values   
% Calculate water flow
fh = @(t,x)fflow(t, x, time, u, p);
[T,X] = ode45(fh, time, x0);
function dxdt = fflow(t,x,time,u,p)
qin = interp1(time,u,t);                  % Volume flux density as input (=irr-ET)
q12 = -K(1)*((-abs(diff1)/p(1))+1);       % Volume flux density layer 1 to 2 [cm3 cm-2 h-1]
q23 = -K(2)*((-abs(diff2)/p(2))+1);       % Volume flux density layer 2 to 3 [cm3 cm-2 h-1]
q34 = -K(3)*((-abs(diff3)/p(3))+1);       % Volume flux density layer 3 to 4 [cm3 cm-2 h-1]
% Compute change in water content theta (x)
dx1dt = (qin-q12)/p(1);       
dx2dt = (q12/p(1)-q23/p(2));
dx3dt = (q23/p(2)-q34/p(3));
dxdt = [dx1dt; dx2dt; dx3dt];
Best wishes
Torsten.
Related Question