MATLAB: Error using odearguments: nputs must be floats, namely single or double.

ode45odearguments

I'm writing up a code to solve a system of mass balances. The balances are ordinary differential equations, so I used ode45 to solve the system. Additionally, inside the function that has the differential equations, I used vpasolve because I have another system of equations for the parameters of one of the differential equation. When I run up the code I got the following error:
Error using odearguments (line 113)
Inputs must be floats, namely single or double.
Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in TrayBioreactor_1 (line 100)
    [t,M]=ode45(@TrayFunction_1,ti,Mi,[],par);
– Part of the code:
% Initial conditions 
      % Initial water activity
      awi=(1-exp(-0.0004373*(par.T+57.926)*par.phii^(1.2626)));
      % Saturation vapor pressure of water (Antoine equation) $Pa$
      par.Psat=133.322*exp(18.3036-(3816.44/(par.T+227.02)));
      % Initial moisture content in the interparticle space air $\frac{kg_{H_2O}}{kg_{da}}$
      Hi=((1/0.0288)*(par.Psat/par.P)*awi)/(55.5556-55.5556*awi*(par.Psat/par.P));
      % Initial total amount of water
      Wi=(1-par.epsilon)*par.rho_s*par.phii+par.epsilon*par.rho_a*Hi;
      double(Wi)
      % Initial oxygen concentration
      Oi=1;
      % Initial carbon dioxide concentration
      Ci=1;
      Mi=[Oi Ci Wi];
      ti=(0:1:500);
      [t,M]=ode45(@TrayFunction_1,ti,Mi,[],par);
Function:
function dM=TrayFunction_1(t,M,par)
% Components
      O=M(1);       % Oxygen concentration $\frac{molO_2}{m^3_{a}}$
      C=M(2);       % Carbon dioxide concentration $\frac{molCO_2}{m^3_{a}}$
      W=M(3);       % Total amount of water in a control volume $\frac{kg_{H_2O}}{m^3_{b}}$
% Oxygen balance
      dO=(1/(par.epsilon*par.At*par.dz))*(par.Ky_O*par.At*(par.O_h-O)-(1-par.epsilon)*par.rho_s*par.At*par.dz*par.r_O);
% Carbon dioxide balance
      dC=(1/(par.epsilon*par.At*par.dz))*((1-par.epsilon)*par.rho_s*par.At*par.dz*par.r_C-par.Ky_C*par.At*(C-par.C_h));
% Complementary equations
      syms aw H phi
      % Water activity
      eqn1=aw-(1-exp(-0.0004373*(par.T+57.926)*phi^(1.2626)))==0;
      % Moisture content of the interparticle space air 
      eqn2=H-((1/0.0288)*(par.Psat/par.P)*aw)/(55.5556-55.5556*aw*(par.Psat/par.P))==0;
      % Moisture content in the solids
      eqn3=(1-par.epsilon)*par.rho_s*phi+par.epsilon*par.rho_a*H==W;
      [aw, H, phi]=vpasolve([eqn1, eqn2, eqn3], [aw, H, phi]);
      double(aw); double(H); double(phi);
% Water balance 
      dW=(1/(par.At*par.dz))*(-par.Ky_eW*par.At*(aw-par.aw_eq)-par.Ky_dW*par.At*(H-par.H_h));
      dM=[dO dC dW]';
return
I really appreciate all your comments.

Best Answer

    double(aw); double(H); double(phi);
These commands don't change the variables aw, H, or phi in any way. They convert those symbolic variables into double precision arrays ... and then throws those double precision arrays away. If you want those variables to become double arrays, call the double function with an output.
awD = double(aw);
HD = double(H);
phiD = double(phi);
Now use those new variables to compute the variable your ODE function will return.
Related Question