How to define stepchanges, variable BC for different timeintervalls using pdepe?

I am trying to define stepchanges of relative humidity at the boundary of a slab varying from 0.9 to 0.3 every 12 hours without any progress. The code is working without any issues for only one single BC but not otherwise. Any idea how I can get it? Thanks!

function parabolicmoisture
global  tp a_v tend v_s delta exi  
L=0.15;     
k=0.028;     
rho=156;
cp=990;
mu=4.25;
D=25*10^-6;
delta=D/mu;
v_s=17.28e-3;
exi=5/0.65;
a_v=(delta*v_s)/exi;
m = 0;
tp=24*3600;
tend=2*tp; 
x = linspace(0,L,1000);
t = linspace(0,tend,2000);


sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
RH = sol(:,:,1);
% A solution profile can also be illuminating.
plot(x,RH(end,:),'r--','LineWidth',2)
hold on
plot(x,RH(3*end/4,:),'b.-','LineWidth',1)
hold on
plot(x,RH(2*end/4,:),'k','LineWidth',1)
hold on
plot(x,RH(1*end/4,:),'*')
title(strcat('Solution at t = ', num2str(tend)))
legend('48h','36h','24h','12h')
xlabel('Distance x')
ylabel('RH (%)')
 
% %Plot surface temperature vs. time
% figure, plot(t,RH(:,1))
% title('Surface RH')
% xlabel('Time (s)')
% ylabel('RH (%)')
% --------------------------------------------------------------
 
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global a_v v_s delta exi
c=1/a_v;
f = 1*DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 0;
% --------------------------------------------------------------

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global  tend 
% 
   if  t<(tend/4)
         pl = ul-0.9;
   elseif  t>=(tend/4) || t<tend/2
          pl = ul-0.3;
   elseif  t>=tend/2  || t<3*tend/4
         pl = ul-0.9;
   elseif t>=3*tend/4
       pl = ul-0.3;
   end
 
ql = 0;
pr = ur;
qr = 0;

[Torsten]

You don't have to use || , but &&.

[Ali Karim]

Thank you for yout comment. My bad of course, I guess I mixed it up while tring different ways. However this is not the issue. It only takes the first defined BC, i.e for the intervall "t<(tend/4)" and neglects the other BC. Any idea?

[Torsten]

Then restart pdepe at tend/4, tend/2 and 3*tend/4 with the new boundary condition.

[Ali Karim]

Thanks again. Could you give an example? I cannot hav eany if or for somewhere else in the code if that is what you mean?

[Torsten]

Take a look at the last comment under

https://de.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results

Answer 1

[Bill Greene]

The reason that your code didn't show a change in the BC with time is due to the programming error that Torsten pointed out. If you had simply made the corrections he suggested (as I did) pdepe would have failed at the time of the first BC change. The problem is the sharp step change in BC value at the transition times. If you simply allow the BC to change over a small, but non-zero, time interval, pdepe will be able to obtain a solution. In the code I show below, I implemented this strategy and arbitrarily picked 10 seconds as the BC transition time. I also used the matlab function interp1 to substantially simplify the coding.

function parabolicmoisture
global  tp a_v tend v_s delta exi  
L=0.15;     
k=0.028;     
rho=156;
cp=990;
mu=4.25;
D=25*10^-6;
delta=D/mu;
v_s=17.28e-3;
exi=5/0.65;
a_v=(delta*v_s)/exi;
m = 0;
tp=24*3600;
tend=2*tp;
d=10; % small time over which bc changes occur
bcTimes=[0 tend/4 tend/4+d tend/2 tend/2+d 3*tend/4 3*tend/4+d tend];
bcVals=[.9 .9 .3 .3 .9 .9 .3 .3];
x = linspace(0,L,100);
t = linspace(0,tend,200);
% test the bc table
%figure; plot(t, interp1(bcTimes, bcVals, t)); return;
bcFunc=@(xl,ul,xr,ur,t) pdex1bc(xl,ul,xr,ur,t,bcTimes,bcVals);
sol = pdepe(m,@pdex1pde,@pdex1ic,bcFunc,x,t);
RH = sol(:,:,1);

figure; plot(t/tend, RH(:,1)); grid;
title('Solution at left end as a function of normalized time.');

% A solution profile can also be illuminating.
plot(x,RH(end,:),'r--','LineWidth',2)
hold on
plot(x,RH(3*end/4,:),'b.-','LineWidth',1)
hold on
plot(x,RH(2*end/4,:),'k','LineWidth',1)
hold on
plot(x,RH(1*end/4,:),'*')
title(strcat('Solution at t = ', num2str(tend)))
legend('48h','36h','24h','12h')
xlabel('Distance x')
ylabel('RH (%)')
 
% %Plot surface temperature vs. time
% figure, plot(t,RH(:,1))
% title('Surface RH')
% xlabel('Time (s)')
% ylabel('RH (%)')
% --------------------------------------------------------------
end
 
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global a_v v_s delta exi
c=1/a_v;
f = 1*DuDx;
s = 0;
end
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 0;
end

% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t,bcTimes,bcVals)
pl = ul-interp1(bcTimes, bcVals, t);
ql = 0;
pr = ur;
qr = 0;
end

Answer 2

[Ali Karim]

Thank you very much. I see the problem with my code now. Thank you Torsten and Bill Greene.