PDEPE boundary condition outputting odd results.

Scott Lines (view profile)

on 22 Mar 2019
Latest activity Answered by Scott Lines

on 22 Mar 2019

Torsten (view profile)

For the boundary condition I've set the left (pl) to a constant value of 280 with the idea being that this value will always be 280, however when i output the results the value is instead 4e10 and appears to be changing. What am I doing incorrectly so that my boundary condition is a) not constant and b) way above the value I've stated.
I've gone over the pdepe page a few times trying to figure it out but unfortauntely haven't made much progress.
Thank you for any help provided.
m = 0;
tspan = linspace(0,365*20*86400,21);
xmesh = linspace(0,20,201);
sol = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
u = sol(:,:,1);
figure, plot(u(2,:),xmesh)
ylabel('Depth (m)')
xlabel('Oxygen Concentration (g/m3)');
set(gca, 'XAxisLocation', 'top')
set(gca, 'YDir','reverse')
function[c,f,s]=pdefun(x,t,u,DuDx)
c = 0.18;
f = 1e-8*DuDx;
s = -1e-9*u;
end
function u=icfun(x)
u = 0;
end
function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
pl = 280;
ql = 1;
pr = ur;
qr = 0;
end

Show 1 older comment
Scott Lines

Scott Lines (view profile)

on 22 Mar 2019 Torsten

on 22 Mar 2019
ok, then
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
Scott Lines

Scott Lines (view profile)

on 22 Mar 2019
hmm maybe I have my initial conditions incorrect, by setting u0 = 280 the entire column is set to 280 instead just the initial u value, instead I'm trying to just have the initial u value at 280 and the rest starting at 0.

R2018b

Torsten (view profile)

on 22 Mar 2019
Edited by Torsten

Torsten (view profile)

on 22 Mar 2019

function main
m = 0;
tspan = linspace(0,365*20*86400,21);
xmesh = linspace(0,20,201);
sol = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
u = sol(:,:,1);
figure, plot(u(2,:),xmesh)
ylabel('Depth (m)')
xlabel('Oxygen Concentration (g/m3)');
set(gca, 'XAxisLocation', 'top')
set(gca, 'YDir','reverse')
end
function[c,f,s]=pdefun(x,t,u,DuDx)
c = 0.18;
f = 1e-8*DuDx;
s = -1e-9*u;
end
function u=icfun(x)
u = 0;
end
function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
u0 = 280;
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
end

Scott Lines

Scott Lines (view profile)

on 22 Mar 2019
Thank you very much, that works perfect. One follow up question if you have time, it would be much appreciated.
If I wanted to change the top boundary condition from 280 to 0 after a certain amount of time, is there an easy way to implement this?
Torsten

Torsten (view profile)

on 22 Mar 2019
Yes, call pdepe twice: first from t=0 to t=t_switch, then from t=t_switch to t=tend. And use the solution obtained at the end of the first time interval as initial profile for the second time interval.

Scott Lines (view profile)

on 22 Mar 2019

Okay I've tried a fair bit to get it to work, but I can't seem to pass the end results to the initial condition for the first row in t_switch, it sounds so simple but something is clearly not working. I've tried extracting the end resutls in a different variable and declaring it as global and bringing it as directly, I believe I'm close, any further help would be greatly appreciated.
m = 0;
tspan = linspace(0,365*5*86400,201);
tswitch = linspace(tspan(end),365*20*86400,201);
xmesh = linspace(0,20,101);
sol = pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tspan);
sol2 = pdepe(m,@pdefun2,@icfun2,@bcfun2,xmesh,tswitch);
u = sol(:,:,1);
u2 = sol2(:,:,1);
figure(1), plot(u(end,:),xmesh)
ylabel('Depth (m)')
xlabel('Oxygen Concentration (g/m3)');
set(gca, 'XAxisLocation', 'top')
set(gca, 'YDir','reverse')
w = u(end,:); % tried to just extract it to a separate variable and that didn't work.
function[c,f,s]=pdefun(xmesh,tspan,u,DuDx)
c = 2;
f = 1e-8*DuDx;
s = -1e-7*u;
end
function u=icfun(x)
u = 0.01;
end
function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,tspan)
u0 = 280;
pl = ul-u0;
ql = 0;
pr = ur;
qr = 0;
end
function[c,f,s]=pdefun2(xmesh,tswitch,u,DuDx)
c = 2;
f = 1e-8*DuDx;
s = -1e-7*u;
end
function u2=icfun2(xmesh)
global u
u2 = u1(end, find(u==xmesh,1,'first')) % also tried to call it this way and no luck.
end
function[pl,ql,pr,qr]=bcfun2(xl,ul,xr,ur,tswitch)
u0 = 0;
pl = ul-u0;
ql = 0;
pr = ur;
qr = 0;
end

Torsten

Torsten (view profile)

on 22 Mar 2019
function main
m = 0;
xmesh = linspace(0,20,101);
tspan = linspace(0,365*5*86400,201);
u0 = 280.0;
icarg = @(x) 0.01*ones(size(x));
sol = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan);
w = sol(end,:);
plot(xmesh,w)
tspan2 = linspace(tspan(end),365*20*86400,201);
u0 = 0.0;
icarg = @(x)interp1(xmesh,w,x);
sol2 = pdepe(m,@pdefun,@(x)icfun(x,icarg),@(xl,ul,xr,ur,t)bcfun(xl,ul,xr,ur,t,u0),xmesh,tspan2);
w2 = sol2(1,:);
hold on
plot(xmesh,w2)
end
function [c,f,s] = pdefun(xmesh,tspan,u,DuDx)
c = 2;
f = 1e-8*DuDx;
s = -1e-7*u;
end
function u = icfun(x,icarg)
u = icarg(x);
end
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t,u0)
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
end