Index in position 1 exceeds array bounds (must not exceed 1)
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi everyone, when I run the code below (or in the attachment), I get the error message:
"
Index in position 1 exceeds array bounds (must not exceed 1).
Error in sym/subsref (line 859)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8,param9,param10,param11,param12,param13)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]
Error in samplecode>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H,I,J,K,LL,MM)
Error in sym>funchandle2ref (line 1291)
S = x(S{:});
Error in sym>tomupad (line 1204)
x = funchandle2ref(x);
Error in sym (line 211)
S.s = tomupad(x);
Error in checkDAEInput (line 25)
eqs = sym(eqs);
Error in sym/decic (line 144)
[eqs, vars] = checkDAEInput(eqs, vars);
Error in samplecode (line 59)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
"
when I debbug, it point me here:
"
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
"
However I don't see anything wrong with this entry. What can I do to fix this error?
The complete code is:
clear all
close all
clc
syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H I J K LL MM
eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;
eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;
eqn3 = FF == (d(x) * b(x)) / a(x);
eqn4 = G == (d(x) * c(x)) / b(x);
eqn5 = H == d(x) * e(x);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];
vars = [a(x); b(x); c(x); d(x); e(x)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H, I, J, K, LL, MM);
A = 2.89e-9;
B = 2.35e-9;
C = 1.69e-9;
D = 1.6e-8;
E = 9.25e-9;
FF = 6.24;
G = 5.68e-5;
H = 5.3e-8;
I = 4.01e-5;
J = 7.21e-05 ;
K = 149;
LL = 84.9;
MM = 3.47;
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H, I, J, K, LL, MM);
vars;
a0 = K * LL;
b0 = (FF * K * LL)/MM;
c0 = (FF * G * K * LL)/(MM)^2;
d0 = sym(3.47);
e0 = H/3.47;
y0est = [a0; b0; c0; d0; e0];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10.0^(-2), 'AbsTol' , 10.0^(-2));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')
clear all
0 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Integration dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!