清空筛选条件
清空筛选条件

ODE15s-Index exceeds matrix dimensions

1 次查看（过去 30 天）
Mona
Mona 2011-4-12
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);

请先登录，再回答此问题。

采纳的回答

Andrew Newell
Andrew Newell 2011-4-12
@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.

更多回答（1 个）

Jan
Jan 2011-4-12
Instead of posting large sections of unformatted code (please read the "Markup" link on this page), the complete error message would be more helpful, because it contains the number of the line, which causes the problem. In addition post this line only.
In general "dbstop if error" helps to inspect the problems: The debugger stops, when the error occurs. Then you can observethe current variables in the Command window, Workspace Browser, etc.
  4 个评论
显示 2更早的评论
Jan
Jan 2011-4-12
Look on top of the "Add an Answer" field. Direct link: http://www.mathworks.com/matlabcentral/answers/help/markup

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by