PDEPE boundary condition outputting odd results.

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          

[Torsten]

Please post the boundary conditions you are trying to set in a mathematical notation.

If you want to set u to be 280 at the left boundary point,

pl = ul - 280;
ql = 0.0;

e.g.

[Scott Lines]

[Torsten]

ok, then

pl = ul - u0;
ql = 0;
pr = ur;
qr = 0;

[Scott Lines]

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.

Answer 1

[Torsten]

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]

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]

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.

Answer 2

[Scott Lines]

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]

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  

Answer 3

[Scott Lines]

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