Left Boundary Condition of the pdepe for solving Spherical Coordinates Diffusion Equation

13 次查看（过去 30 天）

Yannan 2023-8-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004377-left-boundary-condition-of-the-pdepe-for-solving-spherical-coordinates-diffusion-equation

评论： Yannan 2023-8-6

Hello! I'm trying to figure out the 1D spherical cooordinates diffusion equation by using the pdepe in MATLAB. My PDE function is:

du/dt = D/x^2 * d/dx(x^2 * du/dx)

where D is the diffusivity, a constant. u is the concentration of the solution. x is the distance and t is time.

The initial condition is u(x,0)=0 and u(0,0)=C0, where C0 is a constant concentration.

The Boundary conditions is u(0,t) = C0=0.1 and u(100,t)=0.

I wrote the left boundary condition by pl=ul-0.1 because it should be started at 0.1 and decreased as the distance increases. However, it seems that I can't modify the left boundary condition because of the MATLAB setting, but I want to know if there are some methods to change that setting and make the pdepe work?

Thank you!

Here is what I have got right now:

clear;clc;close all
L = 1e-2;
x = linspace(0,L,20);                           % 1e-2m, divided into 20 unit distance
t = linspace(0,600,10);                          % 10min
m = 2;                                          % dc/dt=x^-m * d/dx(x^m * f(x, dc/dx)), for this case, m=2
sol = pdepe(m,@diffpde,@diffic,@diffbc,x,t);    % pde solve program, enter the required m, pde function, 
                                                % inital condition, boundary condition and variables
u = sol(:,:,1);                                 % extract solution from the sol matrix
plot(x,u,x,u,'*')                               % plot concentration u and distance x graph
ylim([0, 0.1]);
xlabel('Distance x')
ylabel('Solution u')
title('Solution profiles at several times')
function [c,f,s] = diffpde(x,t,u,dudx)          % Set up function parameter of the general form:
    c = 1;                                      
    s = 0;
    f = 6e-11 * dudx;                           % diffusivity = 6e-11
end
function u0 = diffic(x)
    if x == 0
        u0 = 0.1;
    else
        u0 = 0;
    end
end
function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t,u)  % boundary condition, pl is the p function value for the left boundary
pl = ul-0.1;                                      % ql is the q function for the left; pr is the p function for the right;
                                                  % qr is the q function for the right
ql = 0;                                           % pl comes from u(0,t)=C0=0.1 --> u(0,t)-0.1=0
pr = ur;                                        
qr = 0;                                           % should be u(0,t)=0, p(ur)-(0 * f(dudx)) = 0, pr=ur=0, qr=0
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

采纳的回答

Bill Greene 2023-8-4

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004377-left-boundary-condition-of-the-pdepe-for-solving-spherical-coordinates-diffusion-equation#answer_1283252

编辑：Bill Greene 2023-8-4

When the model geometry is cylindrical or spherical, the

mathematically well-posed boundary condition at r=0 is the

symmetry condition. And, as you have discovered, pdepe automatically

imposes this condition regardless of what you specify. The only way to

remove this restriction is to modify the pdepe code itself.

As an alternative, I have written a PDE solver, pde1dm, which accepts

essentially the same input as pdepe but does not require a symmetry

boundary condition at r=0; it simply applies whatever BC the user has

defined. I made this programming choice because I decided, specifically, that

the type of BC you are trying to impose may be reasonable in some cases.

I ran pde1dm on your code by simply replacing "pdepe" with "pde1dm".

The code runs without problem, imposing your BC at r=0.

However, I noted that your diffusivity coefficient is very low and that after only

ten minutes, the material hasn't diffused very far.

If you want to try pde1dm, it can be downloaded here.

34 个评论
显示 32更早的评论隐藏 32更早的评论

Torsten 2023-8-4

编辑：Torsten 2023-8-4

在 MATLAB Online 中打开

If you defined the xmesh as linspace(1e-6,1,100), e.g., both of your above functions are correct.

Note that you have to put the 1/a^4 into the c-coefficient for the definition of your PDE.

L = 1e-2;

farpoint = 10;

D = 6e-11;

c0 = 0.1;

xleft = L/farpoint;

xright = 1.0;

x = linspace(xleft,xright,7000);

t = linspace(0,100*24*3600,10); % 100 days

m = 0; % dc/dt=x^-m * d/dx(x^m * f(x, dc/dx)), for this case, m=2

sol = pdepe(m,@(x,t,u,dudx)diffpde(x,t,u,dudx,D,L),@diffic,@(xl,ul,xr,ur,t,u)diffbc(xl,ul,xr,ur,t,c0),x,t); % pde solve program, enter the required m, pde function,

% inital condition, boundary condition and variables

u = sol(:,:,1); % extract solution from the sol matrix

figure(1)

plot(x,u) % plot concentration u and distance x graph

ylim([0, c0]);

xlabel('Distance L/x')

ylabel('Solution v')

title('Solution profiles at several times')

x = L./x;

idx = x<10*L;

x = x(idx);

u = u(:,idx);

figure(2)

plot(x,u)

ylim([0, c0]);

xlabel('Distance x')

ylabel('Solution u')

title('Solution profiles at several times')

function [c,f,s] = diffpde(x,t,u,dudx,D,L) % Set up function parameter of the general form:

c = 1/x^4;

s = 0;

f = D/L^2 * dudx; % diffusivity = 6e-11

end

function u0 = diffic(x)

u0 = 0;

end

function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t,c0) % boundary condition, pl is the p function value for the left boundary

pl = ul; % ql is the q function for the left; pr is the p function for the right;

% qr is the q function for the right

ql = 0; % pl comes from u(0,t)=C0=0.1 --> u(0,t)-0.1=0

pr = ur - c0;

qr = 0; % should be u(0,t)=0, p(ur)-(0 * f(dudx)) = 0, pr=ur=0, qr=0

end

Yannan 2023-8-6

在 MATLAB Online 中打开

@Torsten

Hi, Torsten! I did something like this but the graph looks wired actually.

L = 1e-2;

Di = 1e-13; % Diffusivity of device = 1e-13

Dt = 6e-11; % Diffusivity of tissue = 6e-11

c0 = 0.1;

xleft = 0;

xright = 10 * L;

x = linspace(xleft,xright,5000);

t = linspace(0,8640000,20); % 100 days

ci = c0 * 0.1 * exp(-t); % original concentration in implanted device c=0.1

m = 2;

sol = pdepe(m,@(x,t,u,dudx)diffpde(x,t,u,dudx,Di,Dt,L),@diffic,@(xl,ul,xr,ur,t,u)diffbc(xl,ul,xr,ur,t,c0),x,t); % pde solve program, enter the required m, pde function,

% inital condition, boundary condition and variables

u = sol(:,:,1); % extract solution from the sol matrix

figure(1)

plot(x,u) % plot concentration u and distance a=R/r graph

ylim([0, 0.1]);

xlabel('Distance x')

ylabel('Solution u')

title('Solution profiles at several times')

function [c,f,s] = diffpde(x,t,u,dudx,Di,Dt,L) % Diffusion equation

c = 1;

s = 0;

if x <= 1e-2

f = Di * dudx;

else

f = Dt * dudx;

end

function u0 = diffic(x) % initial condition

if x <= 1e-2

u0 = 0.1;

else

u0 = 0;

end

function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t,c0) % Boundary conditions

pl = 0;

ql = 1;

pr = ur;

qr = 0;

end

Torsten 2023-8-6

编辑：Torsten 2023-8-6

在 MATLAB Online 中打开

Your code looks fine to me. You can check whether xright is big enough by integrating the concentration over the computational domain. It should remain constant, i.e. the right boundary should have no influence on the results near r = L.

By the way: You can check the results obtained in our first step by integrating from L to Inf with m = 2 and boundary conditions C(L) = C0 and C(inf) = 0.

L = 1e-2;

Di = 1e-13; % Diffusivity of device = 1e-13

Dt = 6e-11; % Diffusivity of tissue = 6e-11

c0 = 0.1;

xleft = 0;

xright = 100 * L;

x = linspace(xleft,xright,5000);

x1 = linspace(xleft,L,500);

x2 = linspace(L,xright,4500);

x = [x1,x2(2:end)];

t = linspace(0,1000*24*3600,20); % 1000 days

ci = c0 * 0.1 * exp(-t); % original concentration in implanted device c=0.1

m = 2;

sol = pdepe(m,@(x,t,u,dudx)diffpde(x,t,u,dudx,Di,Dt,L),@(x)diffic(x,L),@diffbc,x,t); % pde solve program, enter the required m, pde function,

% inital condition, boundary condition and variables

u = sol(:,:,1); % extract solution from the sol matrix

figure(1)

plot(x,u) % plot concentration u and distance a=R/r graph

xlim([0 0.1])

ylim([0, 0.1])

xlabel('Distance x')

ylabel('Solution u')

title('Solution profiles at several times')

function [c,f,s] = diffpde(x,t,u,dudx,Di,Dt,L) % Diffusion equation

c = 1;

s = 0;

if x <= L

f = Di * dudx;

else

f = Dt * dudx;

end

function u0 = diffic(x,L) % initial condition

if x <= L

u0 = 0.1;

else

u0 = 0;

end

function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t) % Boundary conditions

pl = 0;

ql = 1;

pr = ur;

qr = 0;

end

Yannan 2023-8-6

在 MATLAB Online 中打开

Hi Torsten @Torsten, this is what I got if I set Di = Dt. The graph becomes nearly disappeared

For the mesh to draw the 3D graph, I'm curious that where we can use that ci = C0 * exp(-t) in the whole program. Becuase as you see the pdepe doesn't need the ci for the boundary condition or initial condition, I think it can be only used in plotting I guess?

L = 1e-2;

Di = 6e-11; % Diffusivity of device = 1e-13

Dt = 6e-11; % Diffusivity of tissue = 6e-11

c0 = 0.1;

xleft = 0;

xright = 100 * L;

x = linspace(xleft,xright,5000);

x1 = linspace(xleft,L,500);

x2 = linspace(L,xright,4500);

x = [x1,x2(2:end)];

t = linspace(0,86400000,20); % 100 days

ci = c0 * 0.1 * exp(-t); % original concentration in implanted device c=0.1

m = 2;

sol = pdepe(m,@(x,t,u,dudx)diffpde(x,t,u,dudx,Di,Dt,L),@(x)diffic(x,L),@diffbc,x,t); % pde solve program, enter the required m, pde function,

% inital condition, boundary condition and variables

u = sol(:,:,1); % extract solution from the sol matrix

figure(1)

plot(x,u) % plot concentration u and distance a=R/r graph

xlim([0 0.1])

ylim([0, 0.1])

xlabel('Distance x')

ylabel('Solution u')

title('Solution profiles at several times')

function [c,f,s] = diffpde(x,t,u,dudx,Di,Dt,L) % Diffusion equation

c = 1;

s = 0;

if x <= L

f = Di * dudx;

else

f = Dt * dudx;

end

function u0 = diffic(x,L) % initial condition

if x <= L

u0 = 0.1;

else

u0 = 0;

end

function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t) % Boundary conditions

pl = 0;

ql = 1;

pr = ur;

qr = 0;

end

Torsten 2023-8-6

The graph becomes nearly disappeared

Use t = linspace(0,864000,20);

I'm curious that where we can use that ci = C0 * exp(-t)

Is ci the mean concentration in the implanted device ? And ci does not depend on Di ? That's hard to believe.

Yannan 2023-8-6

在 MATLAB Online 中打开

Oh, that is the point I think wrong. You're correct, the ci should depend on Di, I thought it just like something decay, but the diffusion equaiton has given the concentration equation already. Sorry, I thought it on the wrong way.

This is what I got by setting t = linspace(0,864000,20). It looks pretty! So it is surely related to the value of the diffusivity and time interval?

L = 1e-2;

Di = 6e-11; % Diffusivity of device = 1e-13

Dt = 6e-11; % Diffusivity of tissue = 6e-11

c0 = 0.1;

xleft = 0;

xright = 100 * L;

x = linspace(xleft,xright,5000);

x1 = linspace(xleft,L,500);

x2 = linspace(L,xright,4500);

x = [x1,x2(2:end)];

t = linspace(0,864000,20); % 100 days

ci = c0 * 0.1 * exp(-t); % original concentration in implanted device c=0.1

m = 2;

sol = pdepe(m,@(x,t,u,dudx)diffpde(x,t,u,dudx,Di,Dt,L),@(x)diffic(x,L),@diffbc,x,t); % pde solve program, enter the required m, pde function,

% inital condition, boundary condition and variables

u = sol(:,:,1); % extract solution from the sol matrix

figure(1)

plot(x,u) % plot concentration u and distance a=R/r graph

xlim([0 0.1])

ylim([0, 0.1])

xlabel('Distance x')

ylabel('Solution u')

title('Solution profiles at several times')

function [c,f,s] = diffpde(x,t,u,dudx,Di,Dt,L) % Diffusion equation

c = 1;

s = 0;

if x <= L

f = Di * dudx;

else

f = Dt * dudx;

end

function u0 = diffic(x,L) % initial condition

if x <= L

u0 = 0.1;

else

u0 = 0;

end

function [pl,ql,pr,qr] = diffbc(xl,ul,xr,ur,t) % Boundary conditions

pl = 0;

ql = 1;

pr = ur;

qr = 0;

end

请先登录，再进行评论。

类别

Mathematics and Optimization Partial Differential Equation Toolbox Geometry and Mesh

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

产品

MATLAB

版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by

Left Boundary Condition of the pdepe for solving Spherical Coordinates Diffusion Equation

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

34 个评论
显示 32更早的评论隐藏 32更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Left Boundary Condition of the pdepe for solving Spherical Coordinates Diffusion Equation

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

34 个评论 显示 32更早的评论隐藏 32更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

34 个评论
显示 32更早的评论隐藏 32更早的评论