ode15s step size warning with ODE system using anonymous functions

Hello all, I have a system of (1D) equations I would like to solve:

Note that the first equation is actually only evaluated at the top of the grid (indicated in the code by subscript 'fs'). From a previous suggestion I went with the method of lines, discretising in space first and then solving the set of ODEs using ode15s. However, I am running into a problem with respect to the time step size. 'Unable to meet integration tolerances without reducing the step size below the smallest value allowed'. For reference, this is the ode function file:

function dydt = ode_system(t,y,nx,T_fs,T_w,rho_fr_bc1,rho_fr_bc2,rhoa_fun,drhoadT_fun,dQdT_fun,ddQdTT_fun,k_fr_fun,cp_fr_fun,D_eff,h_m,W_a,W_fs,L_sub,k_ice)
x_fs = y(1); dx = x_fs/nx; T_fr = y(2:nx-1); rho_fr = y((nx-1)+1:end);
%% Use anonymous functions to calculate k_fr & rho_a with the given initial conditions for rho_fr, T_fr
k_fr = k_fr_fun(rho_fr); rho_a = rhoa_fun(T_fr);
dTfsdx = gradient([T_fr ; T_fs],dx); dTfsdx = dTfsdx(end);
%% First ODE, evaluated at top of grid
dxfsdt = 1./(rho_fr_bc2)*(h_m*(W_a - W_fs) - D_eff*rhoa_fun(T_fs)*dQdT_fun(T_fs)*dTfsdx);
%% Allocating space
dTdx = zeros(size(T_fr)); ddTdxx = dTdx; drhodx = dTdx; dkdx = dTdx; drhoadx = dTdx; dQdx = dTdx; % Q = rho_v/rho_a ddQdxx = dTdx;
%% Derivatives w.r.t. T
drhoadT = drhoadT_fun(T_fr); dQdT = dQdT_fun(T_fr); ddQdTT = ddQdTT_fun(T_fr);
%% Equations making use of boundary conditions
dTdx(1) = (T_fr(3) - T_w)/(2*dx); dTdx(nx-2) = (T_fs - T_fr(nx-3))/(2*dx); ddTdxx(1) = (T_fr(3) - 2*T_fr(2) + T_w)/(dx^2); ddTdxx(nx-2) = (T_fs - 2*T_fr(nx-2) + T_fr(nx-3))/(dx^2); drhodx(1) = (rho_fr(3) - rho_fr_bc1)/(2*dx); drhodx(nx-2) = (rho_fr_bc2 - rho_fr(nx-3))/(2*dx); dkdx(1) = (k_fr(3) - k_ice)/(2*dx); dkdx(nx-2) = (k_fr_fun(rho_fr_bc2) - k_fr(nx-3))/(2*dx); drhoadx(1) = drhoadT(1)*dTdx(1); drhoadx(nx-2) = drhoadT(nx-2)*dTdx(nx-2); dQdx(1) = dQdT(1)*dTdx(1); dQdx(nx-2) = dQdT(nx-2)*dTdx(nx-2); ddQdxx(1) = ddQdTT(1)*ddTdxx(1); ddQdxx(nx-2) = ddQdTT(nx-2)*ddTdxx(nx-2);
%% System of equations only making use of neighbouring internal points
for i = 3:nx-3 dTdx(i-1) = (T_fr(i+1) - T_fr(i-1))/(2*dx); ddTdxx(i-1) = (T_fr(i+1) - 2*T_fr(i) + T_fr(i-1))/(dx^2); drhodx(i-1) = (rho_fr(i+1) - rho_fr(i-1))/(2*dx); dkdx(i-1) = (k_fr(i+1) - k_fr(i-1))/(2*dx); drhoadx(i-1) = drhoadT(i-1)*( dTdx(i-1) ); dQdx(i-1) = dQdT(i-1)*dTdx(i-1); ddQdxx(i-1) = ddQdTT(i-1)*ddTdxx(i-1); end
%% Create ODEs
dTdt = ( 1./(cp_fr_fun(T_fr).*rho_fr) ).*( dkdx.*dTdx + k_fr.*ddTdxx + L_sub.*drhodx ); drhodt = -D_eff.*( drhoadx.*dQdx + rho_a.*ddQdxx);
%% Collect ODEs as output
dydt = [dxfsdt ; dTdt ; drhodt];
end

Now I am not sure what exactly the problem is. I strongly suspect it has to do with the anonymous functions used by the ode function, which are these:

syms rhoa_fun(T_fun) rhov_fun(T_fun) Q_fun(T_fun) k_fr_fun(rho_fun) 
rhoa_fun = @(T_fun) ((p_0 - (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15)))))./(R_s.*T_fun))+((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun)); 
rhov_fun = @(T_fun) (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun); 
k_fr_fun = @(rho_fun) (rho_fun - rho_ice)./(rho_air - rho_ice)*(k_air - k_ice) + k_ice; 
cp_fr_fun = @(rho_fun) (rho_fun - rho_ice)./(rho_air - rho_ice)*(c_p_air - c_p_ice) + c_p_ice; 
drhoadT_fun_sym(T_fun) = diff(rhoa_fun,T_fun); 
drhoadT_fun = matlabFunction(drhoadT_fun_sym); 
Q_fun = @(T_fun) ((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun))./(((p_0 - (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15)))))./(R_s.*T_fun))+((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun))); 
dQdT_fun_sym(T_fun) = diff(Q_fun,T_fun); 
dQdT_fun = matlabFunction(dQdT_fun_sym); ddQdTT_fun_sym(T_fun) = diff(Q_fun,T_fun,2); 
ddQdTT_fun = matlabFunction(ddQdTT_fun_sym);

Although it might also relate to my usage of the gradient function? I use this as a quick fix as I lack information beyond the grid. After ignoring dxfsdt in the ode system, however, the step size warning still pops up.

If anyone has experience with this problem, I would appreciate the help!

Answer 1

[Jasvin]

Hi,

For resolving these issues with tolerances, most likely the ODE solver chosen, in this case, “ode15s” is not well-suited for this problem and instead experimentation should be done to find the right match.

The following guide details the steps involved in finding the right ODE solver,

https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html

Hope this helps!