CAT arguments dimensions are not consistent using mldivide

My code is

    for i=1:100;
        dt=2;

        A=[e11(1,1) e12(1,1) 0 0;
        e21(1,1) e22(1,1) 0 0;
        0 e32(1,1) e33(1,1) 0;
        0 e42(1,1) e43(1,1) e44(1,1)];
        B=[((qf-qs)*dt)+(e11(1,1)*Vwt(1,i))+(e12(1,1)*p(1,i));
        ((Q+qf*hf-qs*hs(1,i))*dt)+(e21(1,1)*Vwt(1,i))+(e22(1,1)*p(1,i));
        ((Q-ar(1,1)*hc(1,1)*qdc(1,1))*dt)+(e32(1,1)*p(1,i))+(e33(1,1)*ar(1,i));
        ((Ps*(Vsd0-Vsd(1,1))*dt)/Td)+(((hf-hw(1,1))*qf*dt)/hc(1,1))+(e42(1,1)*p(1,i))+(e43(1,1)*ar(1,i))+(e44(1,1)*Vsd(1,i))];
        x=mldivide(A,B);
        Vwt(1,i+1)=x(1,1);
        p(1,i+1)=x(2,1);
        ar(1,i+1)=x(3,1);
        Vsd(1,i+1)=x(4,1);
        Vst(1,i+1)=Vt-Vwt(1,i+1);
  end
The error shown is
??? Error using ==> vertcat
CAT arguments dimensions are not consistent.

Error in ==> boiler at 124 B=[((qf-qs)*dt)+(e11(1,1)*Vwt(1,i))+(e12(1,1)*p(1,i)); Please help

Answer 1

[per isakson]

Without values of e12, e21, e22, etc. we cannot find out why this error is thrown.

Here are some links on debugging in Matlab

• Debugging in MATLAB
• Debugging session
• Debugging Approaches

.

Add the following lines before B=[((qf- ...

    disp( size((qf-qs)*dt)+(e11(1,1)*Vwt(1,i))+(e12(1,1)*p(1,i)));
    disp( size((Q+qf*hf-qs*hs(1,i))*dt)+(e21(1,1)*Vwt(1,i))+(e22(1,1)*p(1,i)));
    disp( size((Q-ar(1,1)*hc(1,1)*qdc(1,1))*dt)+(e32(1,1)*p(1,i))+(e33(1,1)*ar(1,i)));
    disp( size((Ps*(Vsd0-Vsd(1,1))*dt)/Td)+(((hf-hw(1,1))*qf*dt)/hc(1,1))+(e42(1,1)*p(1,i))+(e43(1,1)*ar(1,i))+(e44(1,1)*Vsd(1,i)));

[BITS]

e12,e21 etc are defined above

[BITS]

The problem lies in (1,i)

[per isakson]

I cannot see the definitions. copy&paste and run throws

    Undefined function 'e11' for input arguments of type 'double'.

Probably one of those values is different from what you expect.

[BITS]

the full code is

Vd=40; %volume of drum m^3
Vr=37; %volume of riser m^3
Vdc=11; %volume of downcomer m^3
Vt=Vd+Vr+Vdc;%total volume
Ad=20; %drum area at normal operating level m^2
md=95000;%mass of drum in kg
mt=300000; %total metal mass kg
mr=160000; %total riser mass kg
k=25;%friction coefficent in downcomer riser loop
B=0.3;%parameter
Td=12; % residence time of steam min drum s
Cp=0.420*10^-3;%specific heat of steel j/kgdegC
n=100;
t=zeros(1,n);
t(1,1)=0;
qf=50; %feed flow rate kg/s(to be specified by user)
p=zeros(1,n);
p(1,1)=8.5; %steam pressue MPa(to be specified by user)
Vwt=zeros(1,n);
Vwt(1,1)=57.2;%initial total volume of water in m^3
Vst=zeros(1,n);
Vst(1,1)=Vt-Vwt(1,1);
pf=9.142235;%pressure at which feed is entering Mpa(to be specified by user)
hf=1035.1*10^-3;%-0.0003*pf^6+0.019*pf^5-0.47*pf^4+6.1*pf^3-46*pf^2+2.4e+002*pf+5.7e+002;%enthalpy of saturated feed j/kg
hs=zeros(1,n);
hs(1,1)=(0.059*p(1,1)^3-2.3*p(1,1)^2+10*p(1,1)+2.8e+003)*10^-3; %specific enthalpy of steam (J/kg)
hw=zeros(1,n);
hw(1,1)=(0.24*p(1,1)^3-7.5*p(1,1)^2+1.2e+002*p(1,1)+7.1e+002)*10^-3;%specific enthalpy of water (J/kg)
Ps=zeros(1,n);
Ps(1,1)=(0.032*p(1,1)^3-0.53*p(1,1)^2+8.2*p(1,1)-5.9); % density of steam(kg/m^3)
Pw=zeros(1,n);
Pw(1,1)=(-0.056*p(1,1)^3+1.6*p(1,1)^2-32*p(1,1)+9e+002); %density of water(kg/m^3)
U=zeros(1,n);
U(1,1)=((0.032*3)*p(1,1)^2-(0.53*2)*p(1,1)+8.2); %dPs/dp
V=zeros(1,n);
V(1,1)=-(0.056*3)*p(1,1)^2+(1.6*2)*p(1,1)-32; %dPw/dp
W=zeros(1,n);
W(1,1)=((0.059*3)*p(1,1)^2-(2.3*2)*p(1,1)+10)*10^-3; %dhs/dp
X=zeros(1,n);
X(1,1)=((0.24*3)*p(1,1)^2-15*p(1,1)+002)*10^-3; %dhw/dp
ts=zeros(1,n);
ts(1,1)=(0.06*p(1,1)^3-1.9*p(1,1)^2+27*p(1,1)+1.7e+002)+273.15; %temperature of steam(K)
O=zeros(1,n);
O(1,1)=(0.06*3)*p(1,1)^2-(1.9*2)*p(1,1)+27; %dts/dp
g=9.81; %acceleration due to gravity m/s^2
%first system is at equilibrium
qs=qf;
Q=qs*hs(1,1)-qf*hf;%in kW
Vsd0=7.7;
hc=zeros(1,n);
hc(1,1)=hs(1,1)-hw(1,1);%enthalpy of condensation in kJ/kg
Vsd=zeros(1,n);
Vsd(1,1)=Vsd0-((Td*(hw(1,1)-hf)*qf)/(Ps(1,1)*hc(1,1)));%volume of steam in the drum below water level
Lr=28.9;%length of riser in m
Ldc=28.9;%length of downcomer in m
Adc=Vdc/Ldc;% area of downcomer in m^2
Ar=Vr/Lr;% area of riser in m^2
%ii=((2*Pw*Adc*(Pw-Ps)*g*Vr)/k)^0.5;
%jj=(Q/hc);
%tt=(Pw/(Pw-Ps));
%uu=(Ps/(Pw-Ps));
%guess=[0 1];
%result=fsolve(@eqns,guess);
%x0=[0 0];
%fsolve(@myfun,x0)
ar=zeros(1,n);
ar(1,1)=0.051;
aar=zeros(1,n);
aar(1,1)=(Pw(1,1)/(Pw(1,1)-Ps(1,1)))*(1-(Ps(1,1)/((Pw(1,1)-Ps(1,1))*ar(1,1)))*log(1+((Pw(1,1)-Ps(1,1))*ar(1,1))/Ps(1,1)));
neta=zeros(1,n);
neta(1,1)=ar(1,1)*(Pw(1,1)-Ps(1,1))/Ps(1,1);
daarp=zeros(1,n);
daarp(1,1)=(1/(Pw(1,1)-Ps(1,1))^2)*(Pw(1,1)*U(1,1)-Ps(1,1)*V(1,1))*((1+(Pw(1,1)/(Ps(1,1)*(1+neta(1,1)))))-...
    (((Ps(1,1)+Pw(1,1))*log(1+neta(1,1)))/(neta(1,1)*Ps(1,1))));
daarar=zeros(1,n);
daarar(1,1)=((Pw(1,1)/(Ps(1,1)*neta(1,1))))*((log(1+neta(1,1))/neta(1,1))-(1/(1+neta(1,1))));
Vwd=zeros(1,n);
Vwd(1,1)=Vwt(1,1)-Vdc-((1-aar(1,1))*Vr);
lw=zeros(1,n);
lw(1,1)=Vwd(1,1)/Ad;
ls=zeros(1,n);
ls(1,1)=Vsd(1,1)/Ad;
l=zeros(1,n);
l(1,1)=lw(1,1)+ls(1,1);
qsd=zeros(1,n);
qsd(1,1)=(Ps(1,1)*Vsd0/Td);
qdc=zeros(1,n);
qdc(1,1)=((2*Pw(1,1)*Adc*(Pw(1,1)-Ps(1,1))*g*aar(1,1)*Vr)/k)^0.5;
qct=zeros(1,n);
qct(1,1)=(hw(1,1)-hf(1,1))*qf/hc(1,1);
qr=zeros(1,n);
qr(1,1)=qdc(1,1);
e11=zeros(1,n);
e11(1,1)=Pw(1,1)-Ps(1,1);
e12=zeros(1,n);
e12(1,1)=(Vwt(1,1)*V(1,1))+(Vst(1,1)*U(1,1));
e21=zeros(1,n);
e21(1,1)=(Pw(1,1)*hw(1,1))-(Ps(1,1)*hs(1,1));
e22=zeros(1,n);
e22(1,1)=(Vwt(1,1)*(hw(1,1)*V(1,1)+Pw(1,1)*X(1,1)))+(Vst(1,1)*(hs(1,1)*U(1,1)+Ps(1,1)*W(1,1)))-Vt+(mt*Cp*O(1,1));
e32=zeros(1,n);
e32(1,1)=((Pw(1,1)*X(1,1)-ar(1,1)*hc(1,1)*V(1,1))*(1-aar(1,1))*Vr)+...
    ((((1-ar(1,1))*hc(1,1)*U(1,1))+Ps(1,1)*W(1,1))*aar(1,1)*Vr)+((Ps(1,1)+(Pw(1,1)-Ps(1,1))*ar(1,1))*hc(1,1)*Vr*daarp(1,1))-Vr+(mr*Cp*O(1,1));
e33=zeros(1,n);
e33(1,1)=((1-ar(1,1))*Ps(1,1)+ar(1,1)*Pw(1,1))*hc(1,1)*Vr*daarar(1,1);
e42=zeros(1,n);
e42(1,1)=(Vsd(1,1)*U(1,1))+((1/hc(1,1))*(Ps(1,1)*Vsd(1,1)*W(1,1)+Pw(1,1)*Vwd(1,1)*X(1,1)-Vsd(1,1)-...
    Vwd(1,1)+md*Cp*O(1,1)))+(ar(1,1)*(1+B)*Vr*(aar(1,1)*U(1,1)+(1-aar(1,1))*V(1,1)+(Ps(1,1)-Pw(1,1))*daarp(1,1)));
e43=zeros(1,n);
e43(1,1)=ar(1,1)*(1+B)*(Ps(1,1)-Pw(1,1))*Vr*daarar(1,1);
e44=zeros(1,n);
e44(1,1)=Ps(1,1);
%give step input
qs=qs+10;
%Q=Q+10^6;
%solving the dynamic equations
for i=1:100;
    dt=2;

      A=[e11(1,i) e12(1,i) 0 0;
      e21(1,i) e22(1,i) 0 0;
      0 e32(1,i) e33(1,i) 0;
      0 e42(1,i) e43(1,i) e44(1,i)];
      B=[((qf-qs)*dt)+(e11(1,i)*Vwt(1,i))+(e12(1,i)*p(1,i));
      ((Q+qf*hf-qs*hs(1,i))*dt)+(e21(1,i)*Vwt(1,i))+(e22(1,i)*p(1,i));
      ((Q-ar(1,i)*hc(1,i)*qdc(1,i))*dt)+(e32(1,i)*p(1,i))+(e33(1,i)*ar(1,i));
      ((Ps*(Vsd0-Vsd(1,i))*dt)/Td)+(((hf-hw(1,i))*qf*dt)/hc(1,i))+(e42(1,i)*p(1,i))+(e43(1,i)*ar(1,i))+(e44(1,i)*Vsd(1,i))];
      x=mldivide(A,B);
      Vwt(1,i+1)=x(1,1);
      p(1,i+1)=x(2,1);
      ar(1,i+1)=x(3,1);
      Vsd(1,i+1)=x(4,1);
      Vst(1,i+1)=Vt-Vwt(1,i+1);
      hs(1,i+1)=(0.059*p(1,i+1)^3-2.3*p(1,i+1)^2+10*p(1,i+1)+2.8e+003)*10^-3; %specific enthalpy of steam (J/kg)
      hw(1,i+1)=(0.24*p(1,i+1)^3-7.5*p(1,i+1)^2+1.2e+002*p(1,i+1)+7.1e+002)*10^-3;%specific enthalpy of water (J/kg)
      Ps(1,i+1)=0.032*p(1,i+1)^3-0.53*p(1,i+1)^2+8.2*p(1,i+1)-5.9; % density of steam(kg/m^3)
      Pw(1,i+1)=-0.056*p(1,i+1)^3+1.6*p(1,i+1)^2-32*p(1,i+1)+9e+002; %density of water(kg/m^3)
      U(1,i+1)=(0.032*3)*p(1,i+1)^2-(0.53*2)*p(1,i+1)+8.2; %dPs/dp
      V(1,i+1)=-(0.056*3)*p(1,i+1)^2+(1.6*2)*p(1,i+1)-32; %dPw/dp
      W(1,i+1)=((0.059*3)*p(1,i+1)^2-(2.3*2)*p(1,i+1)+10)*10^-3; %dhs/dp
      X(1,i+1)=((0.24*3)*p(1,i+1)^2-15*p(1,i+1)+002)*10^-3; %dhw/dp
      ts(1,i+1)=0.06*p(1,i+1)^3-1.9*p(1,i+1)^2+27*p(1,i+1)+1.7e+002+273.15; %temperature of steam(K)
      O(1,i+1)=(0.06*3)*p(1,i+1)^2-(1.9*2)*p(1,i+1)+27;
      hc(1,i+1)=hs(1,i+1)-hw(1,i+1);%enthalpy of condensation in J/kg
      Vsd(1,i+1)=Vsd0-((Td*(hw(1,i+1)-hf)*qf)/(Ps(1,i+1)*hc(1,i+1)));
      aar(1,i+1)=(Pw(1,i+1)/(Pw(1,i+1)-Ps(1,i+1)))*(1-(Ps(1,i+1)/((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1)))...
          *log(1+((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))/Ps(1,i+1)));
      neta(1,i+1)=ar(1,i+1)*(Pw(1,i+1)-Ps(1,i+1))/Ps(1,i+1);
      daarp(1,i+1)=(1/(Pw(1,i+1)-Ps(1,i+1))^2)*(Pw(1,i+1)*U(1,i+1)-Ps(1,i+1)*V(1,i+1))*((1+(Pw(1,i+1)/(Ps(1,i+1)*(1+neta(1,i+1)))))-...
      (((Ps(1,i+1)+Pw(1,i+1))*log(1+neta(1,i+1)))/(neta(1,i+1)*Ps(1,i+1))));
      daarar(1,i+1)=((Pw(1,i+1)/(Ps(1,i+1)*neta(1,i+1))))*((log(1+neta(1,i+1))/neta(1,i+1))-(1/(1+neta(1,i+1))));
      Vwd(1,i+1)=Vwt(1,i+1)-Vdc-((1-aar(1,i+1))*Vr);
      lw(1,i+1)=Vwd(1,i+1)/Ad;
      ls(1,i+1)=Vsd(1,i+1)/Ad;
      l(1,i+1)=lw(1,i+1)+ls(1,i+1);
      qsd(1,i+1)=(Ps(1,i+1)*Vsd0/Td);
      qdc(1,i+1)=((2*Pw(1,i+1)*Adc*(Pw(1,i+1)-Ps(1,i+1))*g*aar(1,i+1)*Vr)/k)^0.5;
      qct(1,i+1)=((hw(1,i+1)-hf)*qf/hc(1,i+1))+((1/hc(1,i+1)*dt)*(Ps(1,i+1)*Vst(1,i+1)*W(1,i+1)+...
          Pw(1,i+1)*Vwt(1,i+1)*X(1,i+1)-Vt+mt*Cp*O(1,i+1))*(p(1,i+1)-p(1,i)));
      qr(1,i+1)=qdc(1,i+1)-((Vr/dt)*((aar(1,i+1)*U(1,i+1)+(1-aar(1,i+1))* V(1,i+1)+...
          (Pw(1,i+1)-Ps(1,i+1))*daarp(1,i+1)))*(p(1,i+1)-p(1,i)))+...
          ((((Pw(1,i+1)-Ps(1,i+1))*Vr*daarar(1,i+1)*(ar(1,i+1)-ar(1,i+1)))/dt));
      e11(1,i+1)=Pw(1,i+1)-Ps(1,i+1);
      e12(1,i+1)=(Vwt(1,i+1)*V(1,i+1))+(Vst(1,i+1)*U(1,i+1));
      e21(1,i+1)=(Pw(1,i+1)*hw(1,i+1))-(Ps(1,i+1)*hs(1,i+1));
      e22(1,i+1)=(Vwt(1,i+1)*(hw(1,i+1)*V(1,i+1)+Pw(1,i+1)*X(1,i+1)))+...
          (Vst(1,i+1)*(hs(1,i+1)*U(1,i+1)+Ps(1,i+1)*W(1,i+1)))-Vt+(mt*Cp*O(1,i+1));
      e32(1,i+1)=((Pw(1,i+1)*X(1,i+1)-ar(1,i+1)*hc(1,i+1)*V(1,i+1))*(1-aar(1,i+1))*Vr)+...
          ((((1-ar(1,i+1))*hc(1,i+1)*U(1,i+1))+Ps(1,i+1)*W(1,i+1))*aar(1,i+1)*Vr)+...
          ((Ps(1,i+1)+(Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))*hc(1,i+1)*Vr*daarp(1,i+1))-Vr+(mr*Cp*O(1,i+1));
      e33(1,i+1)=((1-ar(1,i+1))*Ps(1,i+1)+ar(1,i+1)*Pw(1,i+1))*hc(1,i+1)*Vr*daarar(1,i+1);
      e42(1,i+1)=(Vsd(1,i+1)*U(1,i+1))+((1/hc(1,i+1))*(Ps(1,i+1)*Vsd(1,i+1)*W(1,i+1)+...
          Pw(1,i+1)*Vwd(1,i+1)*X(1,i+1)-Vsd(1,i+1)-Vwd(1,i+1)+md*Cp*O(1,i+1)))+...
          (ar(1,i+1)*(1+B)*Vr*(aar(1,i+1)*U(1,i+1)+(1-aar(1,i+1))*V(1,i+1)+(Ps(1,i+1)-Pw(1,i+1))*daarp(1,i+1)));
      e43(1,i+1)=ar(1,i+1)*(1+B)*(Ps(1,i+1)-Pw(1,i+1))*Vr*daarar(1,i+1);
      e44(1,i+1)=Ps(1,i+1);
  end

Now run it

[per isakson]

It is possible to attach files. Please use it next time.

I added the four lines disp(size( .. as I propose in the answer above. There was cut&paste errors, which I have fixed. And I put the code in a function, cssm1.

All the four expressions, which you try to concatenate vertically, B=[((qf ... return values of different size.

Running cssm1 on R2013a returned

    >> cssm1
       1.0e+04 *
        3.1353    3.1353
       1.0e+04 *
        2.2456    2.2456
         1     1
         1   100
    Error using vertcat
    Dimensions of matrices being concatenated are not consistent.
    Error in cssm1 (line 131)
          B=[((qf-qs)*dt)+(e11(1,i)*Vwt(1,i))+(e12(1,i)*p(1,i));