MATLAB Answers

PDEPE boundary condition outputting odd results.

2 views (last 30 days)
Scott Lines
Scott Lines on 22 Mar 2019
Answered: 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
Torsten
Torsten on 22 Mar 2019
ok, then
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
Scott Lines
Scott Lines 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.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 22 Mar 2019
Edited: Torsten 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

  2 Comments

Scott Lines
Scott Lines 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 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.

Sign in to comment.

More Answers (2)

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

Torsten
Torsten 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

Sign in to comment.


Scott Lines
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.

Sign in to answer this question.

Products


Release

R2018b

Translated by