MATLAB: Ode15s (…) must return a column vector

MATLABodeode15sode45

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 1×100 vector,
  • q_av(t) when evaluated returns a 1×100 vector,
  • q_a(t) when evaluated returns a 1×100 vector,
  • tspan is a 1×100 vector,
  • lv0 is 1×1,
  • ao0 is 1×1,
The first ODE @(t,V_lv) returns a 100×1 column vector which is expected for both cycles. The second ODE, @(t, V_ao) returns a 100×1 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.

Best Answer

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.
Related Question