Asked by Scott Lines
on 22 Mar 2019 at 2:24

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

Answer by Torsten
on 22 Mar 2019 at 8:30

Edited by Torsten
on 22 Mar 2019 at 8:30

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

Scott Lines
on 22 Mar 2019 at 8:34

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
on 22 Mar 2019 at 8:49

Sign in to comment.

Answer by Scott Lines
on 22 Mar 2019 at 11:15

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
on 22 Mar 2019 at 11:46

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 at 11:56

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## Torsten (view profile)

Direct link to this comment:https://fr.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results#comment_684214

## Scott Lines (view profile)

Direct link to this comment:https://fr.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results#comment_684223

## Torsten (view profile)

Direct link to this comment:https://fr.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results#comment_684226

## Scott Lines (view profile)

Direct link to this comment:https://fr.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results#comment_684232

Sign in to comment.