ode15s (...) must return a column vector

I am attaching a snippet of my code. I do not understand why I get this error and appreciate any help.

clear all
clc

%Define constants
HR = 90/60;
T = 1/HR;
alpha_s = 0.25;
alpha_r = 0.6;
T_s = T*alpha_s;
T_r = T*alpha_r;

stepsize = 100;
%tspan = linspace(0,T,stepsize);
tspan = linspace(0,T,stepsize);
norm = @(t) (1+floor(length(stepsize)*(t - min(tspan))/(max(tspan)-min(tspan))));

V_tot = 2500;
V_lvSys = 130;
V_lvDia = 50;

V_lvMax = V_lvSys;
V_lvMin = V_lvDia;

V_lv = V_lvSys;
V_lv0 = 0.2*V_lvMin;

V_stroke = V_lvMax - V_lvMin;
CO = HR*V_stroke;

p_lvSys = 120;
p_lvDia = 5;
p_lvMax = p_lvSys;
p_lvMin = p_lvDia;

%solve LV equations...
E_m = mean(p_lvDia)/(mean(V_lvDia)-V_lv0);
E_M = mean(p_lvSys)/(mean(V_lvSys)-V_lv0);

%newest version of E_lv piecewise
E_lv = @(t) ((((E_M - E_m)/2)*(1-cos(pi.*t./T_s))+E_m).*(t<T_s)) +...
            ((((E_M-E_m)./2)*(1+cos(pi.*(t-T_s)./(T_r-T_s)))+E_m).*((T_s<=t)&(t<T_r))) +...
            (E_m.*((T_r<=t)&(t<T)));

p_lv = @(t) E_lv(t).*(V_lv - V_lv0);

p_ao = @(t) 0.99.*p_lv(t);
p_vc = @(t) 1.1.*p_lv(t);

R_av = @(t) (mean(p_lvMax)-p_ao(t))./CO;
R_mv = @(t) (p_vc(t)-mean(p_lvMin))./CO;

q_av = @(t) (((p_lv(t)-p_ao(t))/R_av(t)).*(p_lv(t)>p_ao(t)));
q_mv = @(t) ((p_vc(t)-p_lv(t))/R_mv(t)).*(p_vc(t)>p_lv(t));

n=0;
lv0 = 130;
ao0 = (0.25*V_tot);
sa0 = 0.2*V_tot;
sv0 = 0.7*V_tot;
vc0 = 0.075*V_tot;

V_lv = zeros(100,1);
V_ao = zeros(100,1);
%LV=zeros(100,1);
LV=[];
AO=[];
SA=[];
SV=[];
VC=[];
time=[];

while n <=19 %loop number of heartbeats
   n   
        %LV ODE solution
        %lv0 = 130;
        [t,V_lv] = ode15s(@(t,V_lv) q_mv(t) - q_av(t), tspan, lv0);

          LV = [LV;V_lv];
          %lv0 = LV((n+1)*100,1);
          %lv0 = V_lv(end)
          lv0 = V_lv(end); %redefine initial condition as last value of latest solution

          E_ao = p_lvMax/(0.025*V_tot);

          %solve AO equations...

          p_sa = @(t) 0.99*p_ao(t);

          p_aoMean = (mean(p_ao(t)));
          p_saMean = (mean(p_sa(t)));

          R_a = ((p_aoMean - p_saMean)/CO);

          %define q_a
          q_a = @(t) ((p_ao(t) - p_sa(t))./(R_a));   

          %AO ODE solution
          [t,V_ao] = ode15s(@(t,V_ao) q_av(t) - q_a(t), tspan, ao0);

          AO = [AO;V_ao];
          ao0 = V_ao(end);
          ao0 = ao0(:);
          p_ao = @(t) E_ao*transpose(V_ao).*(t/t);

          %solve SA equations
          p_vc = @(t) 1.1*p_lv(t);
          p_sv = @(t) 1.1*p_vc(t);

          R_s = @(t) (p_saMean - p_sv(t))/CO;

          %define q_s
          q_s = @(t) (p_sa(t) - p_sv(t))./R_s(t);
          %SA ODE solution
          [t,V_sa] = ode15s(@(t,V_sa) q_a(t)-q_s(t),tspan, sa0);

          SA = [SA;V_sa];
          sa0 = V_sa(end);

          %Solve SV equations
          R_v = @(t) (p_sv(t)-p_vc(t))/CO;        

          %define q_v
          q_v = @(t) (p_sv(t) - p_vc(t))/R_v(t);

          %SV ODE solution
          [t,V_sv] = ode15s(@(t,V_sv) q_s(t) - q_v(t), tspan, sv0);

          SV = [SV;V_sv];
          sv0 = V_sv(end);

          %Solve VC equations

          %VC ODE solution
          [t,V_vc] = ode15s(@(t,V_vc) q_v(t) - q_mv(t),tspan, vc0);

          VC = [VC;V_vc];
          vc0 = V_vc(end);

          n
          time = [time;(n+1)*t];
     n=n+1;
  end

I am working on this code that attempts to model a cardiovascular system: changes in compartment volume (veins, heart, etc), pressures, elastances, etc. Each cycle of the while loop is a heart beat.

During the first cycle, everything evaluates properly; all ODEs generate solutions.

During the the second cycle, I get a solution for V_lv, the first ODE, but get an error for V_ao, the second ODE. Quickly describing those 2 ODEs:

• q_mv(t) when evaluated returns a 1x100 vector,
• q_av(t) when evaluated returns a 1x100 vector,
• q_a(t) when evaluated returns a 1x100 vector,
• tspan is a 1x100 vector,
• lv0 is 1x1,
• ao0 is 1x1,

The first ODE @(t,V_lv) returns a 100x1 column vector which is expected for both cycles. The second ODE, @(t, V_ao) returns a 100x1 column vector for the first cycle, then throws the error: must return a column vector, for the 2nd cycle through.

The only thing that changes during each cycle is defining p_ao(t). This function defines a pressure and is itself a function of change in Volume in the ao compartment. So, the solution to V_ao is used to define a new p_ao(t) for the subsequent cycle.

I am still fairly amateur when it comes to MATLAB so hopefully it is a simple thing I am doing wrong.

Thank you.

Answer 1

[Walter Roberson]

You have in a loop,

          q_a = @(t) ((p_ao(t) - p_sa(t))./(R_a));   

          %AO ODE solution
          [t,V_ao] = ode15s(@(t,V_ao) q_av(t) - q_a(t), tspan, ao0);

          AO = [AO;V_ao];
          ao0 = V_ao(end);
          ao0 = ao0(:);
          p_ao = @(t) E_ao*transpose(V_ao).*(t/t);

Each time the q_a is assigned to in the loop, the current value of p_ao and p_sa and R_a will be usedl

The first time through, p_ao is from the early line

p_ao = @(t) 0.99.*p_lv(t);

and the first time that the ode15s call is executed everything is okay. But then a few lines after that, you redefine p_ao in terms of V_ao which is the vector of values (one per timepoint) returned by the ode15s call. So after the first call to ode15s, p_ao is redefined to return a vector that happens to be a row vector. Then when q_a is defined again, it uses that new definition of p_ao, so your ode objective returns a row vector of values, which is not permitted. If you had happened to use V_ao instead of V_ao.' then you would have received an error about the lengths not matching.

[William Banks]

Don't I want p_ao(t) to evaluate as a row vector just as p_lv(t) does (for t = tspan)?

I redefine p_ao(t) to be in terms of the vector V_ao. I take its transpose because initially it is a column vector, but I need it to be a row vector. Is this incorrect?

p_ao = @(t) 0.99.*p_lv(t); %1x100 for t = tspan
p_ao = @(t) E_ao*V_ao.*(t/t); %1x100 for t = tspan

Why do they result in different behaviors in ode15s?

edit: I took

p_aoMean = (mean(p_ao(t)));
p_sa = @(t) 0.99*p_ao(t); 
p_saMean = (mean(p_sa(t)));
R_a = ((p_aoMean - p_saMean)/CO);
q_a = @(t) ((p_ao(t) - p_sa(t))./(R_a));

out of the loop ad that allows me to cycle though the loop error free. However, I don't understand why that fixes the problem.

[Walter Roberson]

"I redefine p_ao(t) to be in terms of the vector V_ao"

And that is making p_ao(t) be vector valued.

Perhaps you wanted to do something like interpolate the V_ao values with respect to t.

[William Banks]

Thank you! I think this conversation provided the insight I needed to solve this problem. I understand now. I need to index into V_ao to only use the value corresponding to a given time step in tspan.

I greatly appreciate your help!