PDEPE boundary condition outputting odd results.

3 次查看（过去 30 天）
Scott Lines
Scott Lines 2019-3-22
回答： Scott Lines 2019-3-22
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 个评论
显示 2更早的评论
Torsten
Torsten 2019-3-22
编辑：Torsten 2019-3-22
ok, then
pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;
Scott Lines
Scott Lines 2019-3-22
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2019-3-22
编辑：Torsten 2019-3-22
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 个评论
Scott Lines
Scott Lines 2019-3-22
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 2019-3-22
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.

请先登录，再进行评论。

更多回答（2 个）

Scott Lines
Scott Lines 2019-3-22
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 个评论
Torsten
Torsten 2019-3-22
编辑：Torsten 2019-3-22
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

请先登录，再进行评论。

Scott Lines
Scott Lines 2019-3-22
Thank you very much, I really appreciate the effort and time you've taken to help me :)

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by