Piecewise function in a DAE
22 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello, please can you help me with the following question:
I've been working in a code (m. file) , a DAEs problem , the code is long, but summarizing , I have a picewise function and I tried to include the picewise according to: https://la.mathworks.com/matlabcentral/answers/408447-how-to-include-a-piecewise
But i got the following error:
Index in position 1 exceeds array bounds (must not exceed 7).
Error in DAEs3 (line 8)
DUst = in2(8,:);
Please anyone can tell me what Im doing wrong, Thank you!!!
Here is the code:
syms t Us(t) Rs(t) HComb(t) Ms(t) MCO2(t) MCaO(t) Mg(t) Tg(t) Ts(t) Tw(t) Tm(t) g(t) yee(t) dAGSdz1(t) dAGWdz1(t) dAWSdz1(t) Dh1(t) AG1(t)
%Code ........
emi_CO2=piecewise(0<=Tg<1300, 0.2, 1300<=Tg<=4300, 0.3);
emi_H2O=piecewise(0<=Tg<1300, 0.15, 1300<=Tg<=4300, 0.12);
emiG=emi_CO2+emi_H2O-emi_CO2.*emi_H2O ;
K=(1-emiG).*(1-emiW).*(fiWS.*(1-emiG).*(1-emiS)+(1-fiWS));
emiWS=(emiW.*emiS.*(1-emiG))./(1-K);
emiGS=(emiG.*emiS.*(1+fiWS.*(1-emiW).*(1-emiG)))./(1-K);
emiGW=(emiG.*emiW.*(1+fiWS.*(1-emiS).*(1-emiG)))./(1-K);
%more code
eqs=[diff(Us(t),t)==((((pgi-PCO2g).*kCaCO3.*6.*((1-Us(t)).^(2./3)))./(xCO2ks.*rhoks1.*RCO2.*Ts(t).*dp.*wp)))
diff(Rs(t),t)==diff(Us(t),t).*(leyCaCO3/100).*Mks1.*xCO2ks.*deltahCO2
diff(HComb(t),t)==Fbr1.*huBr.*(diff((1-exp(-a.*(t-L1).^b)),1))
%more equations]
vars=[Us(t), Rs(t), HComb(t), Tg(t), Ts(t), Tw(t), Tm(t), Ms(t), MCO2(t), MCaO(t), Mg(t)];
origVars = length(vars)
origeqns = length(eqs)
M = incidenceMatrix(eqs, vars)
[eqns, vars] = reduceDifferentialOrder(eqs, vars)
isLowIndexDAE(eqs,vars)
[DAEs,DAEvars] = reduceDAEIndex(eqs,vars)
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)
pDAEs = symvar(DAEs)
pDAEvars = symvar(DAEvars)
extraParams = setdiff(pDAEs, pDAEvars)
ff = daeFunction(DAEs, DAEvars, g(t), yee(t), dAGSdz1(t), dAGWdz1(t), dAWSdz1(t), Dh1(t), AG1(t),'file','DAEs3.m');
g=@(t) interp1(ft,f,t);
yee=@(t) interp1(ft,ye,t);
AG1=@(t) interp1(ft,AG,t);
dAGSdz1=@(t) interp1(ft,dAGSdz,t);
dAGWdz1=@(t) interp1(ft,dAGWdz,t);
dAWSdz1=@(t) interp1(ft,dAWSdz,t);
Dh1=@(t) interp1(ft,Dh,t);
F=@(t,Y,YP) ff(t,Y,YP,g(t),yee(t),dAGSdz1(t), dAGWdz1(t), dAWSdz1(t), Dh1(t), AG1(t));%,'file','DAEs3.m');
t0 = 0;
y0est = [0.25;0.25*Mks1.*xCO2ks.*deltahCO2;Fbr1.*huBr;Tgas+273.15;1110;Tsol+273.15;370];
yp0est= [0;0;0;0;0;0;0];
[tSol, ySol]=ode15i(F, [t0:1:L1], y0est, yp0est)
0 Kommentare
Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!