Solving DAE - Index may be greater than one.
Afficher commentaires plus anciens
This is my program to solve DAEs:
function btp_25Jan4
syms cm(x) Pe(x) g(x) delpie(x) Rad(x) Re Sc de Rr C0 D A B kp t Peo delp Rm Rg Mw kad n L;
eqs = [(1/(2*(cm-1)))*(g-(g^2)/5)*diff(cm,1)-(1-3*g/10)*diff(g,1)==(8/(Re*Sc*de/L))*(8-Pe*g)/(g^2),...
Pe*cm*Rr*g + 8*(cm-1)==(de/(C0*D))*g*(A*cm/(1+B*cm))*kp*(exp(-kp*t)),...
Pe==Peo*(1-delpie/delp)/( 1+Rad/Rm),...
delpie==2*Rg*C0*Rr*cm/Mw,...
Rad==kad*((A*cm/(1+B*cm))*(1-exp(-kp*t)))^n];
vars = [cm(x); Pe(x); g(x); delpie(x); Rad(x)];
Mu = 0.001;
Re = 600;
D = 1.4*(10^-9);
Sc = Mu/D;
de = (4*0.08/pi)^0.5;
Rr = 0.72;
C0 = 4*(10^-3);
A = 0.26;
B = 130;
kp = 0.58;
t = 2*3600;%2 hours
Peo = 24;
delp = 400000;
Rm = 1/(0.001*1.6*10^-11);
Rg = 8.314;
Mw = 20000;
kad = 4*(10^13);
n = 0.15;
L = 0.145;
go = (128/(3*Re*Sc*de/L))^(1/3);
[eqs, vars, ~] = reduceDifferentialOrder(eqs, vars);
[DAEs,DAEvars] = reduceDAEIndex(eqs,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)
DAEs = subs(DAEs);
f = daeFunction(DAEs, DAEvars);
F = @(x, Y, YP) f(x, Y, YP);
y0est = [0; 1];% have used only 2 because after reduceRedundancies only 2 DAEvars and DAEs are left
yp0est = zeros(2,1);
opt = odeset('RelTol', 10.0^(-2), 'AbsTol' , 10.0^(-2));
[y0, yp0] = decic(F, 0.0005, y0est,[], yp0est,[], opt);
end
But I get this error
Error using decic>sls (line 170)
Index may be greater than one.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);
Error in btp_25Jan4 (line 48)
[y0, yp0] = decic(F, 0.0005, y0est,[], yp0est,[], opt);
Please help me to know where have I done the mistake. Thanks.
3 commentaires
akshay
le 12 Mar 2016
hey were you able to resolve the error
I am getting same error...don't know how to proceed
Ghaith Salah
le 6 Mai 2018
Im having the same problem, kindly let me know if you resolved the error
Ramanuja Srinivasan
le 31 Août 2020
Hi, have you sorted out this problem, because i am encountering the same error. I would really appreciate your help.
Réponses (1)
Szar
le 13 Oct 2017
1 vote
Try fixing some of the initial guesses.
Catégories
En savoir plus sur Mathematics dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!