Step Function as Input Boundary Condition to MATLAB PDEPE solver

1 view (last 30 days)
Siamack on 19 Mar 2018
Commented: Bill Greene on 22 Mar 2018
Dear Matlab Users I have the following code to find the transient solution to the 1D heat conduction equation. The code below intends to use a rectangular heat pulse as the initial temperature but i receive an error ( Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t). Can anyone tell me what is issue with the boundary condition as I set it and if there is a better way to deal with the matter.
function out = parabolic_silicon(~)
global rho cp k T_i T_o
L = 10; % [cm]
k = 1.3; % [W/cmK]
rho = 2.33; % [g/cm^3]
cp = 0.7; % [J/gK]
T_i = 295; % [K]
T_o = T_i + 5.2e-3; % [K]
t_end= 1; % [s]
m = 0;
x = linspace(0,L,500);
t = linspace(0,t_end,500);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
Temperature = sol(:,:,1);
L = strcat('t = ',num2str(t(:)),' Seconds');
plot(x,sol);
% legend(L,'location','northeast');
% % --------------------------------------------------------------
% OUTPUT
out = {x,t,sol};
% % --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global rho cp k
c = rho*cp;
f = k*DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 295;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global T_i T_o
pl = ul - (T_i + (T_o*(heaviside(t) - heaviside(t - 0.5))));
ql = 0;
pr = ur - T_i;
qr = 0;

Bill Greene on 22 Mar 2018
I suggest you smooth the sharp corners in your boundary condition profile using a Heaviside function similar to the code below. Note, you can change the sharpness of the corners by varying the r variable.
function k=heavisideSmooth(t)
k1 = 0; k2 = 1; c = 0;
r = 1e3;
ecr = exp(c*r);
erx = exp(r*t);
k = (k1*ecr + k2*erx)./(ecr+erx);
end
Bill Greene on 22 Mar 2018
The error message indicates the ODE solver wasn't able to take even a tiny, first step. So I assumed that was due to the sharp discontinuity in the BC at t=0.