Unable to perform assignment because the left and right sides have a different number of elements.

Hi all, I don't know why does my code give an error message:

Unable to perform assignment because the left and right sides have a different number of elements.

    Error in eulerian>odesfun (line 115)
  dydt(1) = (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2 ;

Error in eulerian (line 53)
dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in,  NSO2,  NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);

The code is given below, and on the attachment:

function eulerian()

global S_total C_total Ca_total CSO2_bulkslurry CSO3 ...
    CCO2_bulkslurry CCaCO3 CCaSO3 CH CCO2_in...
    E NSO2 E_CO2 NCO2  RCaSO3 RCaCO3 SO2_g CO2_g 

global V_Headspace F CSO2_in R Temp HSO2 kga kLa DCa2 ...
    DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 ...
    MWCaCO3 Kad  KCO2 KHCO3 KSO2 KHSO3 Kw  CH0 

% Constant Parameters
KSO2        = 6.24;
KHSO3       = 5.68e-5;
DCa2        = 1.39e-9;
DSO2        = 2.89e-9;
R           = 8.314;
Temp        = 323.15;
HSO2        = 149;
kga         = -4.14e-6;
kLa         = -7e-2;
V_Headspace = 1.5e-6;
F           = 1.67e-4;
CSO2_in     = 6.51e-2;
KCO2        = 1.7e-3;
KHCO3       = 6.55e-8;
DCO2        = 3.53e-9;
kLa_CO2     = 9.598e-4;
HCO2        = 5.15e+3;
CCO2_in     = 0;
KSPCaSO3    = 1.07e-7;
BETCaSO3    = 10;
ktot        = 8.825e-3;
BETCaCO3    = 12.54;
MWCaCO3     = 100.0869;
Kad         = 0.84;
CH0         = 7.414e-6;
Kw          = 5.3e-8;

tspan = [600;1200;1800;2400;3000;10200;17400;24600;31800;39000;...
       46200;53400;60600;67800;75000;82200;89400;96600;103800;111000;...
       118200;125400;132600;139800;147000;154200;161400;168600;175800;183000];

y0 = [1.72E-02; 6.97E-08; 3.34E-01; 4.88E-06; 4.88E-06; 4.86E+01; 0];
y0 = y0';

delta = 0.01;
nloop = 30;
t0 = 600.;
t = t0;

for i = 1:nloop

dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in,  NSO2,  NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);

y = y0+ dydt.'*delta;

%Variables  SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3 are passed as y(1)...y(7)        
 SO2_g       = y(1);
 CO2_g       = y(2);
 S_total     = y(3);
 C_total     = y(4);
 Ca_total    = y(5);
 CCaCO3      = y(6);
 CCaSO3      = y(7);

% CH, CSO2_bulkslurry , E, NSO2, CCO2_bulkslurry, E_CO2, NCO2, CSO3, RCaSO3 and RCaCO3 are passed as variable Parameters

CH              = fsolve(@(CH) CH + 2.* Ca_total - ((S_total.* KSO2.*CH)/(CH^2 + KSO2.*CH + KSO2.*KHSO3))- 2.*((S_total.*KSO2.*KHSO3)/(CH^2 ...
                + KSO2.*CH + KSO2.*KHSO3))-((C_total.*KCO2.*CH)/(CH^2 + KCO2.*CH + KCO2.*KHCO3))- 2.*((C_total.*KCO2.*KHCO3)/(CH^2 ...
                + KCO2.*CH + KCO2.*KHCO3))- Kw/CH, CH0);

CSO2_bulkslurry = (S_total.* CH.^2)/(CH.^2 + KSO2.* CH + KSO2.* KHSO3);

E               = 1 + DCa2.* Ca_total / (DSO2.* S_total);

NSO2            =(SO2_g.* R.* Temp - HSO2.* CSO2_bulkslurry) / ( (1/kga) + HSO2/ (E.* kLa) ) ;

CCO2_bulkslurry = (C_total.* CH.^2)/(CH.^2 + KCO2.* CH + KCO2.* KHCO3);

E_CO2           = 1 + DCa2  * Ca_total / (DCO2 .* C_total);

NCO2            = kLa_CO2 .* E_CO2 .* (CO2_g .* R .* Temp / HCO2 - CCO2_bulkslurry);

CSO3            = (S_total.* KSO2 .* KHSO3)/(CH.^2 + KSO2 .* CH + KSO2 .* KHSO3);

RCaSO3          =  0.001 .* 1.62e-3 .* exp(-5152/ Temp) .* ( (Ca_total .* CSO3 / KSPCaSO3) - 1).^3 / BETCaSO3;

RCaCO3          = - ktot .* BETCaCO3 .* MWCaCO3 .* CCaCO3 .* CH .* (1 - (Kad .* CH)/(1+ Kad .* CH));

y0 = y;
t = t + delta;
CH0 = CH;

[t,y] = ode45(@odesfun,tspan,y0)

figure (1)

plot (t,y)

legend('Mod-SO_{2-gas}', 'Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}',...
    'Mod-CaSO_{3}.^{1}/_{2}H_{2}O','Mod-H^+','Mod-pH', 'Location','best')
end

end

function dydt = odesfun(t,y,V_Headspace, F, CSO2_in, CCO2_in,  NSO2,  NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g)

%Variables  SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3 are passed as y(1)...y(7)        

dydt = zeros(7,1);

% SO2 in gas phase
dydt(1) = (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2 ;

% CO2 in gas phase.
dydt(2) = (1 ./V_Headspace).* (F .* CCO2_in - F .* CO2_g) - NCO2 ;

% Total sulfur balance
dydt(3) = NSO2 -  RCaSO3;    

% Total sulfur balance
dydt(4) = RCaCO3 - NCO2;

% Total calcium balance
 dydt(5) =  RCaCO3 - RCaSO3;

% Limestone dissolution submodel
 dydt(6) = RCaCO3;

% Hannebachite crystalization submodel    
 dydt(7) = RCaSO3;

end

Please help.

[Kevin Chng]

Try to change this

 dydt(1) = (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2 ;

to

size( (1 ./V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2)

To see what is the size of output, it might not a single variable

[Walter Roberson]

Unable to perform assignment because the left and right sides have a different number of elements.

I see that and the first thing I think of is "Uninitialized global"

global S_total C_total [...]

And look at that, a bunch of globals.

I would not bother going any further until the code had been rewritten to avoid using global.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Dursman Mchabe]

Thanks a lot for the comments, Kevin and Walter. When use

dydt = zeros(7,1);

% SO2 in gas phase
size((1 /V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2);

% CO2 in gas phase.
size((1 /V_Headspace).* (F .* CCO2_in - F .* CO2_g) - NCO2);

% Total sulfur balance
size(NSO2 -  RCaSO3);    

% Total sulfur balance
size(RCaCO3 - NCO2);

% Total calcium balance
 size(RCaCO3 - RCaSO3);

% Limestone dissolution submodel
 size(RCaCO3);

% Hannebachite crystalization submodel    
 size(RCaSO3);

I get the error message:

Error using  - 
Matrix dimensions must agree.

Error in eulerian>odesfun (line 118)
size((1 /V_Headspace).* (F .* CCO2_in - F .* CO2_g) - NCO2);

Error in eulerian (line 53)
dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in,  NSO2,  NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);

Again, when I comment out the globals, the variable parameters become undefined, I will get errors like:

Undefined function or variable 'NSO2'.

Error in eulerian (line 53)
dydt = odesfun(t,V_Headspace, F, CSO2_in, CCO2_in,  NSO2,  NCO2, RCaSO3, RCaCO3,SO2_g, CO2_g,y0);

[Walter Roberson]

You would not just comment out the globals: you would add each of them (or a struct or array of them) as parameters to the function, and use parameterization http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Dursman Mchabe]

Thanks a lot for the comment Walter. I have re-written the code, and I have added each parameter at the function for odes. To simplify further, I have introduced CH as a set of experimental values, instead of fsolve equation. The code now looks like this:

function eulerian1()

% Constant Parameters
KSO2        = 6.24;
KHSO3       = 5.68e-5;
DCa2        = 1.39e-9;
DSO2        = 2.89e-9;
R           = 8.314;
Temp        = 323.15;
HSO2        = 149;
kga         = 4.14e-6;
kLa         = 7e-2;
V_Headspace = 1.5e-6;
F           = 1.67e-4;
CSO2_in     = 6.51e-2;
KCO2        = 1.7e-3;
KHCO3       = 6.55e-8;
DCO2        = 3.53e-9;
kLa_CO2     = 9.598e-4;
HCO2        = 5.15e+3;
CCO2_in     = 0;
KSPCaSO3    = 1.07e-7;
BETCaSO3    = 10;
ktot        = 8.825e-3;
BETCaCO3    = 12.54;
MWCaCO3     = 100.0869;
Kad         = 0.84;
CH0         = 7.414e-6;
Kw          = 5.3e-8;

tspan = [600;1200;1800;2400;3000;10200;17400;24600;31800;39000;...
       46200;53400;60600;67800;75000;82200;89400;96600;103800;111000;...
       118200;125400;132600;139800;147000;154200;161400;168600;175800;183000];

y0 = [1.72E-02; 6.97E-08; 3.34E-01; 4.88E-06; 4.88E-06; 4.86E+01; 0];
y0 = y0';

delta = 0.01;
nloop = 30;
t0 = 600.;
t = t0;

for i = 1:nloop

y = y0 + odesfun(t,y,Parameter).'*delta;

%Variables  SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3 are passed as y(1)...y(7)        
 SO2_g       = y(1);
 CO2_g       = y(2);
 S_total     = y(3);
 C_total     = y(4);
 Ca_total    = y(5);
 CCaCO3      = y(6);
 CCaSO3      = y(7);

% CH, CSO2_bulkslurry , E, NSO2, CCO2_bulkslurry, E_CO2, NCO2, CSO3, RCaSO3 and RCaCO3 are passed as variable Parameters

CH              = [7.41E-06;9.33E-06;1.20E-05;1.05E-05;1.74E-05;3.72E-05;3.55E-05;1.00E-04;4.07E-02;2.45E-01;...
                  6.17E-01;1.32E+00;2.29E+00;2.34E+00;2.40E+00;1.82E+00;1.38E+00;2.09E+00;1.82E+00;1.58E+00;...
                  2.29E+00;1.62E+00;1.12E+00;8.91E-01;8.51E-01;7.59E-01;8.71E-01;1.12E+00;1.00E+00;8.51E-01];

CSO2_bulkslurry = (S_total.* CH.^2)/(CH.^2 + KSO2.* CH + KSO2.* KHSO3);

E               = 1 + DCa2.* Ca_total / (DSO2.* S_total);

NSO2            =(SO2_g.* R.* Temp - HSO2.* CSO2_bulkslurry) / ( (1/kga) + HSO2/ (E.* kLa) );

CCO2_bulkslurry = (C_total.* CH.^2)/(CH.^2 + KCO2.* CH + KCO2.* KHCO3);

E_CO2           = 1 + DCa2  * Ca_total / (DCO2 .* C_total);

NCO2            = kLa_CO2 .* E_CO2 .* (CO2_g .* R .* Temp / HCO2 - CCO2_bulkslurry);

CSO3            = (S_total.* KSO2 .* KHSO3)/(CH.^2 + KSO2 .* CH + KSO2 .* KHSO3);

RCaSO3          =  0.001 .* 1.62e-3 .* exp(-5152 / Temp) .* ( (Ca_total .* CSO3 / KSPCaSO3) - 1).^3 / BETCaSO3;

RCaCO3          = - ktot .* BETCaCO3 .* MWCaCO3 .* CCaCO3 .* CH .* (1 - (Kad .* CH)/(1+ Kad .* CH));

y0 = y;
t = t + delta;
CH0 = CH;

[t,y] = ode45(@odesfun,tspan,[V_Headspace,F,SO2_in,CCO2_in,NSO2,NCO2,RCaSO3,RCaCO3,SO2_g,CO2_g],y0)

figure (1)

plot (t,y)

legend('Mod-SO_{2-gas}', 'Mod-CO_{2-gas}','Mod-S_{total}','Mod-C_{total}','Mod-Ca^{2+}_{total}','Mod-CaCO_{3}',...
    'Mod-CaSO_{3}.^{1}/_{2}H_{2}O', 'Location','best')
end

end

function dydt = odesfun(t,y,Parameter)
%Parameters
V_Headspace     = Parameter(1); 
F               = Parameter(2);
CSO2_in         = Parameter(3);
CCO2_in         = Parameter(4);
NSO2            = Parameter(5);
NCO2            = Parameter(6);
RCaSO3          = Parameter(7);
RCaCO3          = Parameter(8);
SO2_g           = Parameter(9);
CO2_g           = Parameter(10);
CH              = Parameter(11);
CSO2_bulkslurry = Parameter(12);
E_CO2           = Parameter(13);
NCO2            = Parameter(14);
CSO3            = Parameter(15);
RCaSO3          = Parameter(16);
RCaCO3          = Parameter(17);

dydt = zeros(7,1);

% SO2 in gas phase
dydt(1) =   (1 /V_Headspace).* (F.* CSO2_in - F.* SO2_g ) - NSO2;

% CO2 in gas phase.
dydt(2) =   (1 /V_Headspace).* (F .* CCO2_in - F .* CO2_g) - NCO2;

% Total sulfur balance
dydt(3) =   NSO2 -  RCaSO3;    

% Total sulfur balance
dydt(4) =   RCaCO3 - NCO2;

% Total calcium balance
dydt(5) =   RCaCO3 - RCaSO3;

% Limestone dissolution submodel
dydt(6) =   RCaCO3;

% Hannebachite crystalization submodel    
dydt(7) =   RCaSO3;

end

When I run the code, I now get an error message, that says

Undefined function or variable 'y'.

Error in eulerian1 (line 45)
y = y0 + odesfun(t,y,Parameter).'*delta;

Indicating that I should define y and Parameter somewhere in the eulerian1 function. Where can I define it and how should I define it? Am I in the right track? Regards Dursman

[Torsten]

Instead of correcting technical errors in your code, it might be better to first explain what you are trying to do. From your code above, I can't tell that.

[Dursman Mchabe]

Thanks a lot your comment Torsten. What I'm trying to do is to solve 7 ODEs,

dydt(1) ... dydt(7).

These ODEs have variable parameters:

CH, CSO2_bulkslurry , E, NSO2, CCO2_bulkslurry, E_CO2, NCO2, CSO3, RCaSO3 and RCaCO3.

In turn, these variable parameters consist of constant parameters and variables. Of which the constant parameters are:

V_Headspace F CSO2_in R Temp HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3 BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad  KCO2 KHCO3 KSO2 KHSO3 Kw  CH0

The 7 variables are:

SO2_g, CO2_g, S_total, C_total, Ca_total, CCaCO3 and CCaSO3.

All I have to do is to solve the 7 ODEs over a period of 183000 seconds.

The equation to evaluate one of the variable parameters, CH ,is rather complicated, and requires fsolve.

CH = fsolve(@(CH) CH + 2.* Ca_total - ((S_total.* KSO2.*CH)/(CH^2 + KSO2.*CH + KSO2.*KHSO3))- 2.*((S_total.*KSO2.*KHSO3)/(CH^2 
                + KSO2.*CH + KSO2.*KHSO3))-((C_total.*KCO2.*CH)/(CH^2 + KCO2.*CH + KCO2.*KHCO3))- 2.*((C_total.*KCO2.*KHCO3)/(CH^2 + KCO2.*CH + KCO2.*KHCO3))- Kw/CH, CH0);

As a result, I thought of just using a vector of experimental values.

CH = [7.41E-06;9.33E-06;1.20E-05;1.05E-05;1.74E-05;3.72E-05;3.55E-05;1.00E-04;4.07E-02;2.45E-01; 6.17E-01;1.32E+00;2.29E+00;2.34E+00;2.40E+00;1.82E+00;1.38E+00;2.09E+00;1.82E+00;1.58E+00; 2.29E+00;1.62E+00;1.12E+00;8.91E-01;8.51E-01;7.59E-01;8.71E-01;1.12E+00;1.00E+00;8.51E-01];

I hope I unpacked it sufficiently. In a nutshell, I am trying to use ode45 to solve 7 ODEs.