How fix this error: "Error CASSENS2 returns a vector of length 3, but the length of initial conditions vector is 4."."

function xBL = calculaxBL(T)
global P A B C Psat a1 r qsi rohfo ropoliol phi1

P = 0.921;
A = 6.820820;
B= 1550.590847;
C = 17.08326013;
Psat =  10^(A-B/(T+C))/101.235;
a1 = P/Psat;  
r = 3.0975; 
qsi= 1.588; 
rohfo=1380;
ropoliol=1100;

   [X,FVAL,EXITFLAG] = fzero(@(x) log(1-x)+x*(1-1/r)+qsi*x^2-log(a1), [0.4 0.999]); % Flory-Huggins resolution to obtain phi2
  if (EXITFLAG==1)
       phi2 = X;
  else
      print "erro na fzero EXITFLAG =" EXITFLAG 
  end 
  phi1 = 1 - phi2;
  xBL = 1/(1+(1-phi2)*ropoliol/(phi2*rohfo*r));

[Birdman]

Where is CASSENS2 ?

[brbecas]

    function dydt = cassens2(t,y)

    global c1 c2 c3 c4 c5 c6 d1 d2  lambda W mhfo MMhfo MMpoliol R VT mpoliol ropoliol k1 k2 k3 k4 k5 AOH EOH OH0 NCO0 rnco rw XOH AW EW rohfo rop L XW W0 U Text  T k6 k7 k8 xBL NCO OH
    T    = y(1);
    XOH = y(2);
    XW  = y(3);
    L   = y(4);
    c1=1800; 
    c2=1; 
    c3=836.6; 
    c4=0;  
    c5=870; 
    L=0.242;
    c6=1210; 
    d1=252400; 
    d2=86000; 
    rop=1100;
    rohfo=1380; 
    ropoliol=1100; 
    lambda=0.168;
    OH0=2419700;
    W0=708.89; 
    W=0.0000001; 
    NCO0=162711; 
    EW=44100; 
    AW=25.8; 
    AOH=2.6; 
    rnco=NCO0/OH0;
    rw=W0/OH0;
    XOH=0.00001;
    XW=0.00001;  
    MMhfo=164; 
    MMpoliol=508; 
    mhfo=0.05;
    mpoliol=0.1457; 
    VT=0.0002; 
    Text = 298; 
    EOH=43800;
    R=8.31; 
    U = 0.0003;
    k1 = c1+c2*c3+c4*c5+L*c6;
    k2 = -d1*OH0/rop;
    k3 = -d2*W0/rop;
    OH = OH0*(1-XOH);
    NCO = (rnco-2*rw*XW-XOH)*(1/(1+(L*rop/rohfo)));
    rnco = NCO0/OH0;
    rw = W0/OH0;
    k4 = AOH*OH*NCO;
    k5 = AW*(1-XW)*(1/(1+(L*rop/rohfo)));
    k6= (VT-((mpoliol/ropoliol)+(mhfo/rohfo)));
    k7 = mpoliol/MMpoliol;
    k8 = MMhfo/MMpoliol;

    % calculating dxBL/dT using finite differenced with deltaT = 0.5 C
    dxBLdT = calculaxBL(T+0.5)-calculaxBL(T-0.5);
    dXOHdt = k4*exp(-EOH/(R*y(1)));
    dXWdt = k5*exp(-EW/(R*y(1)));
    dTdt = (k2/k1)*dXOHdt+(k3/k4)*dXWdt - U*(T-Text);
    dLdt = k8*(1./(1-xBL)^2)*dxBLdT*dTdt;

    dydt =[  dTdt ; dXOHdt ; dXWdt ; dLdt  ] ;
    end

Answer 1

[Walter Roberson]

One of your global variables is probably not initialized, so it is probably empty.

You should avoid using global variables. You should https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[brbecas]

This helped. But I keep finding some errors...maybe you can understand what is the problem.

Error using fzero (line 274) The function values at the interval endpoints must differ in sign.

Error in calculaxBL (line 35) [x fval exitflag] = fzero(fun,x0,options);

function xBL = calculaxBL(T)
global P A1 A2 A3 A4 Pc Tc Psat a1 r qsi rohfo ropoliol  

P = 0.921; 
A1=-7.93972; 
A2=1.8547; 
A3=-2.93418; 
A4=-5.48416; 
Pc = 2901; 
Tc = 444.42;
Psat = Pc*exp(Tc/T*(A1*(1-T/Tc)+A2*(1-T/Tc)^1.5+A3*(1-T/Tc)^2.5+A4*(1-T/Tc)^5)); 
a1 = P/Psat;  
r = 3.0975; 
qsi= 1.588;  
rohfo=1380;
ropoliol=1100;

fun = @(x) log(1-x)+x*(1-1/r)+qsi*x^2-log(a1); 
x0 = [0.4 0.999]; 
options = optimset('Display','iter'); 
[x fval exitflag] = fzero(fun,x0,options);

phi2 = x;
%phi1 = 1 - phi2;

xBL = 1/(1+(1-phi2)*ropoliol/(phi2*rohfo*r));
end

[Walter Roberson]

Your code is suggestive that you are using ode45() or similar, which in turn would tend to suggest that T might start at 0. However, your function divides by T so Psat would have a divide by 0 there. The limit of fun() as T approaches 0 is undefined, being +inf on one side and -inf on the other. We therefore need to know what the valid range is for T to be able to solve this.

If you take your code and put in a symbolic T and symbolic x, then you can see that fun has an exp() in T with no x in it, and you can see that there are some terms in x and a term log(1-x) . For positive T < 22221/50 (about 444) the exp() term is real-valued and so cannot balance out the complex component of log(1-x) for x>1 . There just might be T values that exactly balance out the complex component for x>1 but if so they are very very narrow -- not zero crossings if they exist but rather just touching 0 with a jump of more than 1E41 for going even a hair over (arctan with a large multiplier changes quadrants)

For 0 < x < 1 then solutions can be found for up to about T = 444.4199999999522 with the top end seeming to be about x = .9999670374 . More typical solutions are for T roughly 225