MATLAB Answers

0

PDEPE boundary condition outputting odd results.

Asked by Scott Lines on 22 Mar 2019
Latest activity Answered by Scott Lines on 22 Mar 2019
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

  4 Comments

Show 1 older comment
ok, then
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
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.

Sign in to comment.

Products


Release

R2018b

3 Answers

Answer by Torsten
on 22 Mar 2019
Edited by Torsten
on 22 Mar 2019
 Accepted Answer

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

  2 Comments

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?
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.

Sign in to comment.


Answer by Scott Lines 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

  1 Comment

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

Sign in to comment.


Answer by Scott Lines on 22 Mar 2019

Thank you very much, I really appreciate the effort and time you've taken to help me :)

  0 Comments

Sign in to comment.