ode15s step size warning with ODE system using anonymous functions

4 views (last 30 days)
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!

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!