Error of function (eps: class must be double) trying to solve a symbolic DAE system with ODE15i
8 views (last 30 days)
Show older comments
Many greetings to anyone who stumbles upon this thread. I'm currently trying to solve thisalgebraic-differential equations system for a chemical reaction rate regression analysis. [OH] is a variable x(t) and [DCF] is my variable y(t). All the other terms are simple constants.
And to do so I stumbled upon this guide for solving DAEs symbolically: https://www.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html
For those that have tried this approach I'd like to ask for help. During Step 4 where I'm instructed to find optimal initial conditions for my variables I get the following errors:
Error using eps
Class must be 'single' or 'double'.
Error in odenumjac (line 88)
br = eps(classF) ^ (0.875);
Error in ode15ipdupdate (line 25)
[dfdy,dfdy_options.fac,NF] = odenumjac(odefun,{t,y,yp,extras{:}},f,dfdy_options);
Error in ode15ipdinit (line 76)
ode15ipdupdate([],odefun,t0,y0,yp0,f0,dfdy,dfdyp,dfdy_options,dfdyp_options,extras);
Error in decic (line 70)
ode15ipdinit(odefun,t0,y0,yp0,res,options,varargin);
Error in parametrows (line 63)
[y0,yp0] = decic(F,0,y0est,[],yp0est,[])
My code is pretty much a copy and paste of that guide, except with my own equations. I honestly don't know how to debug or fix it. I get the same error of the eps function even if I don't try to find optimal conditions and rather jump straight to ode15i. Any idea why? Here's my code:
datos %loading of data
error_residual %loading of other data
%here I begin
syms x(t) y(t) FeOOH OZ ka kb kc kd ke kf kg kh ki kj kk kl k_o k_oh k_nine
eqn1 = ka*x(t)^3+x(t)^2*(kb*y(t)+kc*FeOOH+kd*FeOOH/y(t)-ke*OZ)+x(t)*(kf*OZ*y(t)+kg*(OZ^2)-kh*OZ*FeOOH+ki*(OZ*FeOOH/y(t)))+(-kj*OZ^2*y(t)-kk*OZ^2*FeOOH+kl*(OZ^2*FeOOH/y(t))) == 0;
eqn2 = diff(y(t),1) == -k_o*OZ*y(t)-k_oh*y(t)*x(t)-k_nine*FeOOH*OZ;
eqns = [eqn1 eqn2];
vars = [x(t);y(t)];
origVars = length(vars);
M = incidenceMatrix(eqns,vars);
[DAEs, DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs,pDAEvars);
f = daeFunction(DAEs,DAEvars,FeOOH,OZ,k_o,k_nine,k_oh,ka,kb,kc,kd,ke,kf,kg,kh,ki,kj,kk,kl);
k_o = double(x_parameters(1)); %x_parameters is the array with the data for the constants
k_oh = double(x_parameters(2));
k_9 = double(x_parameters(3));
ka =double(x_parameters(4));
kb = double(x_parameters(5));
kc =double(x_parameters(6));
kd =double(x_parameters(7));
ke = double(x_parameters(8));
kf =double(x_parameters(9));
kg = double(x_parameters(10));
kh = double(x_parameters(11));
ki = double(x_parameters(12));
kj = double(x_parameters(13));
kk = double(x_parameters(14));
kl = double(x_parameters(15));
FeOOH = double(g(1)); %likewise with x_parameters, g is data array
OZ = double(ox(1));
F = @(t,Y,YP) f(t,Y,YP,FeOOH,OZ,k_o,k_nine,k_oh,ka,kb,kc,kd,ke,kf,kg,kh,ki,kj,kk,kl);
y0est = [0;0.054768223996110]; %this numbers are my initial numerical guesses for the initial conditions
yp0est = [0;-0.012990851139745];
[y0,yp0] = decic(F,0,y0est,[],yp0est,[])
I literally followed that guide step by step, as well as the one for solving DAEs using Mass matrix: https://www.mathworks.com/help/symbolic/solve-daes-using-mass-matrix-solvers.html. I get the same error with this one.
3 Comments
Answers (0)
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!