MATLAB: ODE15s-Index exceeds matrix dimensions

ode15s

Hello all,

I am trying to solve 10 differential equations using ODE15s. But,

I got this error.

??? Index exceeds matrix dimensions.

Can someone help me with this code and tell me the mistake?

The mfiles are below.

function f=sofc(t,x) 
global u 
global m 
f=zeros(10,1); 
m=fsolve(@algebraeq,[5.35379707129652,0.380690723666416,0.266214812909653]); 
Rload=u; 
current=m(1); 
ethaa=m(2); 
ethac=m(3); 
dA=10^-4; 
dC=0.5*10^-4; 
dE=1.8*10^-4; 
L=0.04; 
B=0.04; 
F=96485.3; 
HA=10^-3; 
HC=10^-3; 
Taircatin=950; 
Tfuelanin=950; 
vaircatin=78.98; 
vfuelanin=2.88; 
rhocp=10^6; 
deltaHR=-241830; 
betha1=3.34*10^4; 
betha2=1.03*10^4; 
alpha=25; 
Deltaanode=5*10^-5; 
Deltacathode=5*10^-5; 
epsiloncat=5*10^-1; 
thoucat=3; 
epsilonan=5*10^-1; 
thouan=3; 
deltaan=5*10^-7; 
deltacat=5*10^-7; 
nuH2=7.07; 
nuo2=16.6; 
nuN2=17.9; 
nuH2o=12.7; 
R=8.314; 
Mo2=31.994; 
MH2=2.02; 
MH2o=18.02; 
MN2=28.02; 
naircatin=3.8*10^-2; 
nfuelanin=1.39*10^-3; 
yo2in=0.2; 
yN2in=0.8; 
vaircat=vaircatin; 
vfuelan=vfuelanin; 
ao2=34.57; 
bo2=0.1078e-2; 
do2=-784586; 
ah2o=34.36; 
bh2o=0.0627e-2; 
ch2o=0.56012e-5; 
ah2=27.67; 
bh2=0.3386e-2; 
an2=27.17; 
bn2=0.4180e-2; 
yH2in=0.9; 
yH2oin=0.1; 
deltaGR=-175933; 
deltaSR=-57; 
Tref=1300; 
%physical parameters variations 
Do2N2=1.013*10^-7*(x(1))^1.75*(1/Mo2+1/MN2)^0.5/ 
((x(2)*x(1)*9.86923267*10^-6+x(5)*R*9.86923267*10^-6*(1-x(4)/ 
x(3)))*(nuo2^(1/3)+nuN2^(1/3))^2); 
Dko2=97*deltacat*abs(sqrt(x(1)/Mo2)); 
Do2=1/(1/Do2N2+1/Dko2); 
ko2cat=(epsiloncat/(thoucat*Deltacathode))*Do2; 
Ho2catin=ao2*(Taircatin-298.15)+bo2*(Taircatin^2/2-298.15^2/2)-do2*(1/ 
Taircatin-1/298.15); 
HN2catin=an2*(Taircatin-298.15)+bn2*(Taircatin^2/2-298.15^2/2); 
Haircatin=yo2in*Ho2catin+yN2in*HN2catin; 
Ho2cat=ao2*(x(5)/x(3)-298.15)+bo2*((x(5)/x(3))^2/2-298.15^2/2)-do2*(1/ 
(x(5)/x(3))-1/298.15); 
HN2cat=an2*(x(5)/x(3)-298.15)+bn2*((x(5)/x(3))^2/2-298.15^2/2); 
Haircat=x(4)/x(3)*Ho2cat+(1-x(4)/x(3))*HN2cat; 
DkH2=97*deltaan*abs(sqrt(x(1)/MH2)); 
DH2H2o=1.013*10^-7*x(1)^1.75*(1/MH2+1/MH2o)^0.5/ 
((x(6)*x(1)*9.86923267*10^-6+x(7)*x(1)*9.86923267*10^-6)*(nuH2^(1/3)+nuH2o^¬(1/3))^2); 
DH2=1/(1/DH2H2o+1/DkH2); 
kH2an=(epsilonan/(thouan*Deltaanode))*DH2; 
HH2anin=ah2*(Tfuelanin-298.15)+bh2*(Tfuelanin^2/2-298.15^2/2); 
HH2oanin=ah2o*(Tfuelanin-298.15)+bh2o*(Tfuelanin^2/2-298.15^2/2)+ch2o*(Tfue-lanin^3/3-298.15^3/3)-241814; 
Hfuelanin=yH2in*HH2anin+yH2oin*HH2oanin; 
HH2an=ah2*(x(10)/x(8)-298.15)+bh2*((x(10)/x(8))^2/2-298.15^2/2); 
HH2oan=ah2o*(x(10)/x(8)-298.15)+bh2o*((x(10)/ 
x(8))^2/2-298.15^2/2)+ch2o*((x(10)/x(8))^3/3-298.15^3/3)-241814; 
Hfuelan=x(9)/x(8)*HH2an+(1-x(9)/x(8))*HH2oan; 
DkH2o=97*deltaan*abs(sqrt(x(1)/MH2o)); 
DH2oH2=DH2H2o; 
DH2o=1/(1/DH2oH2+1/DkH2o); 
kH2oan=(epsilonan/(thouan*Deltaanode))*DH2o; 
%Variables definition 
No2=ko2cat*(x(4)-x(2)/(R)); 
NH2=kH2an*(x(9)-x(6)/(R)); 
NH2o=kH2oan*(x(7)/(R)-x(8)*(1-x(9)/x(8))); 
U0=-1/(2*F)*(deltaGR-deltaSR*(x(1)-Tref)+R*x(1)*log((x(7)*x(1))/ 
(x(6)*x(1)*((x(2)*x(1))/(1.013*10^5))^0.5))); 
rhoE=1/betha1*exp(betha2/x(1)); 
% algebraic equations 

v=U0-ethaa-ethac-rhoE*dE*current/(L*B); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%physical properties of the solid 
cph2solid=ah2+bh2*x(1); 
cpo2solid=ao2+bo2*x(1)+do2*x(1)^(-2); 
cpo2cat=ao2+bo2*x(1)+do2*x(1)^(-2); 
cpn2cat=an2+bn2*x(5)/x(3); 
cpaircat=x(4)/x(3).*cpo2cat+(1-x(4)/x(3)).*cpn2cat; 
cph2an=ah2+bh2*x(10)/x(8); 
cph2oan=ah2o+bh2o*x(10)/x(8)+ch2o*(x(10)/x(8))^2; 
cpfuelan=x(9)/x(8)*cph2an+(1-x(9)/x(8))*cph2oan; 
% cpfuelan=31.4; 
%Differential equations 
f(1)=1/(rhocp*(dA+dC+dE))*((-deltaHR/(2*F)-v)*current/(L*B)+(alpha 
+cph2solid/(2*F)*current/(L*B))*(x(10)/x(8)-x(1))+(alpha+cpo2solid/ 
(4*F)*current/(L*B))*(x(5)/x(3)-x(1))); 
f(2)=R/Deltacathode*(No2-current/(L*B*4*F)); 
f(3)=1/(L*B*HC)*(naircatin-vaircat*x(3)*B*HC-No2*L*B); 
f(4)=1/(L*B*HC)*(naircatin*yo2in-vaircat*x(4)*B*HC-No2*L*B); 
f(5)=1/(L*B*HC*cpaircat)*(naircatin*Haircatin-vaircat*x(3)*Haircat*B*HC 
+L*B*alpha*(x(1)-x(5)/x(3))-No2*Ho2cat*L*B); 
f(6)=R/Deltaanode*(NH2-(1/(L*B))*(current/(2*F))); 
f(7)=R/Deltaanode*(-NH2o+(1/(L*B))*(current/(2*F))); 
f(8)=1/(L*B*HA)*(nfuelanin-vfuelan*x(8)*B*HA+L*B*(NH2o-NH2)); 
f(9)=1/(L*B*HA)*(nfuelanin*yH2in-vfuelan*x(9)*B*HA-NH2*L*B); 
f(10)=1/(L*B*HA*cpfuelan)*(nfuelanin*Hfuelanin- 
vfuelan*x(8)*Hfuelan*B*HA+alpha*(x(1)-x(10)/x(8))*L*B+L*B*(HH2oan*NH2o- 
HH2an*NH2)); 
function g=algebraeq(m) 
global u 
global x 
Rload=u; 
dE=1.8*10^-4; 
L=0.04; 
B=0.04; 
F=96485.3; 
gammaA=5.7*10^7; 
gammaC=7*10^9; 
thetaaC=1.4; 
thetacC=0.6; 
thetaaA=2; 
thetacA=1; 
betha1=3.34*10^4; 
betha2=1.03*10^4; 
EA=140000; 
EC=160000; 
deltaGR=-175933; 
deltaSR=-57; 
Tref=1300; 
U0=-1/(2*F).*(deltaGR-deltaSR.*(x(1)-Tref)+R.*x(1).*log((x(7).*x
(1))./ 
(x(6).*x(1).*((x(2).*x(1))./(1.013*10^5))^0.5))); 
rhoE=1/betha1*exp(betha2/x(1)); 
% algebraic equations 
v=U0-m(2)-m(3)-rhoE*dE*m(1)/(L*B); 
g(1)=m(1)/(L*B)-gammaA*((x(6)*x(1))/(1.013*10^5))*((x(7)*x(1))/ 
(1.013*10^5))*exp(-EA/(R*x(1)))*(exp(thetaaA*F/(R*x(1))*m(2))-exp(- 
thetacA*F/(R*x(1))*m(2))); 
g(2)=m(1)/(L*B)-gammaC*((x(2)*x(1))/(1.013*10^5)).^0.25*exp(-EC/ 
(R*x(1)))*(exp(thetaaC*F*m(3)/(R*x(1)))-exp(-thetacC*F*m(3)/ 
(R*x(1)))); 
g(3)=v-Rload*m(1); 
and then, I have this script to integrate. 
global u 
u=5; 
options=odeset('RelTol',1e-10,'AbsTol',1e-10); 
tspan=[0 100]; 
x0=[1053.96206726841,19.7782866165184,12.0239706084605,2.40128132003310,
11462.7463172666,88.1877747305364,12.2062806809450,12.0659722222222
,
10.6185407334819,12057.5693431174,5.35379707129652,0.380690723666416,0.26621481290965¬3]; 
[t x]=ode15s(@sofc,tspan,x0,options);

Best Answer

@Mona, if you are writing files this long you should learn how to use the debugger (try this tutorial). Your error message should give you a file and line number for where the error occurred, e.g.,

Error in ==> algebraeq at 43

You should be able to just click on the underlined location to go to that line. Then click beside the line number for that line (you should see a red dot). Finally, run the script again. It should stop at the line number and you can look at the values of the variables by hovering your mouse over them.

Related Solutions

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

Maybe you mean

function main
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
F_FEED = 100;           %TOTAL FLOWRATE OF THE FEED (mol/s)
T0 = 650;          %INLET TEMPERATURE OF THE REACTOR (K)
P0 = 303975;    %INLET PRESSURE OF THE REACTOR (Pa)
WEIGHT_cat = 100;       %CATALYST WEIGHT (kg)
R = ToBeSpecified;       
Ua = ToBeSpecified;
Feed = zeros(8,1);
Feed(1) = 0.014*F_FEED;           %INLET FLOWRATE OF G (mol/s)
Feed(2) = 0.634*F_FEED;           %INLET FLOWRATE OF H2O (mol/s)
Feed(3) = 0.352*F_FEED;           %INLET FLOWRATE OF N2 (mol/s)
Feed(8) = T0; 
F_G0 = Feed(1);                 
F_H2O0 = Feed(2);              
F_N20 = Feed(3);                
F_A0 = Feed(4);               
F_AC0 = Feed(5);                
F_AA0 = Feed(6);                
F_FA0= Feed(7);                
F_T0 = F_FEED;
[W,F] = ode45(@Reactors , [0 WEIGHT_cat], Feed);
end
function dFdW = Reactors(t,F)
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
dFdW = zeros(8,1);
%dFdW = [r_G]
%       [r_H2O]
%       [r_N2]
%       [r_O2]
%       [r_A]
%       [r_AC]
%       [r_AA]
%       [r_FA]
%       [T]     Temperature
%reaction
Hf_G = -577.9;                    %Heat of formation: G (kJ/mol)
Hf_H2O = -241.818;                %Heat of formation: H2O (kJ/mol)
Hf_N2 = 0;                        %Heat of formation: N2 (kJ/mol)
Hf_A = -65;                       %Heat of formation: A (kJ/mol)
Hf_AC = -370.06;                  %Heat of formation: AC (kJ/mol)
Hf_AA = -166.1;                   %Heat of formation: AA (kJ/mol)
Hf_FA = -108.6;                   %Heat of formation: FA (kJ/mol)
Cp_G = -35.3267443 + 0.645106658*F(8) - 0.00061741*F(8)^2 + 0.000000239229*F(8)^3;               %provide integrated form
Cp_H2O = 32.95266198 - 0.00391977*F(8) + 0.0000222314*F(8)^2 - 0.0000000104953*F(8)^3;
Cp_N2 = 28.883 -0.00157*F(8)+ 0.00000808*F(8)^2 -0.000000002871*F(8)^3;
Cp_A = -14.9617445 + 0.373240705*F(8) - 0.00034559*F(8)^2 + 0.000000129748*F(8)^3;
Cp_AC = -11.7479865 + 0.402868987*F(8) - 0.00031782*F(8)^2 + 0.000000109142*F(8)^3;
Cp_AA = -5.82005966 + 0.23787576*F(8) - 0.00017362*F(8)^2 + 0.0000000543415*F(8)^3;
Cp_FA = 19.10519375 + 0.048680586*F(8) + 0.00000878405*F(8)^2 - 0.0000000150776*F(8)^3;
integraldeltaCp_rxn_1 = 2*(32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)) + (-14.9617445*(F(8)-298) + 0.5*0.373240705*(F(8)^2-298^2) - (1/3)*0.00034559*(F(8)^3-298^3) + (1/4)*0.000000129748*(F(8)^4-298^4)) - ((-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4))); 
integraldeltaCp_rxn_2 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (-11.7479865*(F(8)-298) + 0.5*0.402868987*(F(8)^2-298^2) - (1/3)*0.00031782*(F(8)^3-298^3) + (1/4)*0.000000109142*(F(8)^4-298^4)) - (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)));
integraldeltaCp_rxn_3 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
integraldeltaCp_rxn_4 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
Hrxn_1 =  2*Hf_H2O + Hf_A - Hf_G + (integraldeltaCp_rxn_1);               %Heat of reaction: reaction 1
Hrxn_2 =  2*Hf_H2O + Hf_AC - Hf_G + (integraldeltaCp_rxn_2);              %Heat of reaction: reaction 2
Hrxn_3 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_3);        %Heat of reaction: reaction 3
Hrxn_4 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_4);        %Heat of reaction: reaction 4
%INITIAL CONCENTRATIONS  
VOLUMETRICFLOWRATEin = F_FEED*R*T0;                         %Initial volumetric flowrate into the reactor (m^3/s) (assuming ideal gas in)   NEGLECTING P/P0 to start...
C_G0 =      F_G0/VOLUMETRICFLOWRATEin;                      %Inital concentration of G
C_H2O0 =    F_H2O0/VOLUMETRICFLOWRATEin;                    %Inital concentration of H2O
C_N20 =     F_N20/VOLUMETRICFLOWRATEin;                     %Inital concentration of N2
C_A0 =  0;                                                  %Outlet concentration of A

C_AC0 = 0;                                                  %Outlet concentration of AC

C_AA0 = 0;                                                  %Outlet concentration of AA

C_FA0 = 0;                                                  %Outlet concentration of FA

C_T0 =      F_FEED/VOLUMETRICFLOWRATEin;                    %Initial total concentration
%FINAL Concentrations
F_T = F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7);       
C_G =   C_T0*(F(1)/F_T)*(T0/F(8));    %Outlet concentration of G              
C_H2O = C_T0*(F(2)/F_T)*(T0/F(8));      %Outlet concentration of H2O
C_N2 =  C_T0*(F(3)/F_T)*(T0/F(8));      %Outlet concentration of N2              NEGLECTING P/P0 to start...
C_A =   C_T0*(F(4)/F_T)*(T0/F(8));    %Outlet concentration of A
C_AC =  C_T0*(F(5)/F_T)*(T0/F(8));      %Outlet concentration of AC
C_AA =  C_T0*(F(6)/F_T)*(T0/F(8));      %Outlet concentration of AA
C_FA =  C_T0*(F(7)/F_T)*(T0/F(8));      %Outlet concentration of FA
%k-VALUES 
k1 = 23.3244*exp(-45007.1/(R*F(8)));           %Reaction rate constant k1 (m^3/kgcat*s)
k2 = 2528.8*exp(-79244.89/(R*F(8)));           %Reaction rate constant k2 (m^3/kgcat*s)
k3 = 2150.595*exp(-73702.78/(R*F(8)));         %Reaction rate constant k3 (m^3/kgcat*s)
k4 = 0.146768*exp(-24596.97/(R*F(8)));         %Reaction rate constant k4 (m^3/kgcat*s)
%reaction rates
r_1 = -k1*C_G;                                     %Rate of reaction 1
r_2 = -k2*C_G;                                     %Rate of reaction 2
r_3 = -k3*C_G;                                     %Rate of reaction 3
r_4 = -k4*C_G;                                     %Rate of reaction 4
r_G = -(k1+k2+k3+k4)*C_G;                  %Rate of formation of G
r_N2 = 0;                                  %Rate of formation of N2
r_H2O = k1*(C_G)^2+(k3+k4)*G_G;            %Rate of formation of H2O
r_A = k1*C_G;                              %Rate of formation of A
r_AC = k2*C_G;                             %Rate of formation of AC
r_AA = (k3+k4)*C_G;                        %Rate of formation of AA
r_FA = (k3+k4)*C_G;                        %Rate of formation of FA
%Differential Equations to Solve
dFdW(1) = (r_G);
dFdW(2) = (r_H2O);
dFdW(3) = (r_N2);
dFdW(4) = (r_A);
dFdW(5) = (r_AC);
dFdW(6) = (r_AA);
dFdW(7) = (r_FA);
dFdW(8) = (Ua(Ta-T0)+(-r_1)*(-Hrxn_1)+(-r_2)*(-Hrxn_2)+(-r_3)*(-Hrxn_3)+(-r_4)*(-Hrxn_4))/(F(1)*Cp_G + F(2)*Cp_H2O + F(3)*Cp_N2 + F(4)*Cp_A + F(5)*Cp_AC + F(6)*Cp_AA + F(7)*Cp_FA); %dTdW from pg 545 in Fogler;   
end

MATLAB: Simplify expression after collect function

call vpa() on the expression.

Best Answer

Related Solutions

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

MATLAB: Simplify expression after collect function

Related Question