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

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

Answer 1

[Bill Greene]

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.

[Yannan]

Hello! Thank you so much for your response!! I've clicked your link but it shows 404 not found there, I think there is something wrong for the web, could you plz fix it or give me another link? Again, thank you very much!

[Bill Greene]

Sorry, I think the link is fixed now.

[Yannan]

Hello! The link works, thank you so much for the help. And I got a quick question, the only things I need to do are: simply replacing the pdepe in my code to pde1dm and reset the time interval longer, the pde1dm will give me the graph of u and x right?

[Torsten]

@Yannan

You cannot impose a condition as u = 0.1 at r = 0. You must always have a corresponding area for that a Dirichlet condition can have an effect, and the area in a spherical geometry at r = 0 is 0.

[Bill Greene]

Yes, just change pdepe to pde1dm and make sure pde1dm is on your matlab path.

As far changing the time interval, I don't know anything about the physics of your problem so can't say whether that makes sense or not. I would double-check the value of the diffusion coefficient.

[Yannan]

For 10min here is just I'm trying to see if the pdepe works. Because the pdepe has the limitation of left boundary condition, the graph it shows is strange. The time can be 10 mins, 1 day, 5 days, 10 days..etc.!

[Yannan]

@Torsten Hi, Torsten, thank you for your response! This question is not such heat transfer question with Dirichlet condition. It is more likely that I have a constant concentration surface at r=0, and the solution diffuses to the space and I want to figure out the concentration profile of solution with time and distance.

[Torsten]

The area of the surface over which your substance could diffuse is 0 at r = 0. So there is nothing that could diffuse into the sphere over it.

[Yannan]

@Torsten

Ok, I think I understand what is your meaning! I think your point is also a question and really important. In fact, my code is a simplified way to think about the whole question. The entire question is this:

Here I have a spherical implanted with radius = R (constant), and concentration = C0 = 0.1, the solution diffuses from the surface of the spherical implanted to the infinity area out of the implanted. The diffusion equation is the PDE I gave above. Could you please tell me how to solve the PDE according to this boundary conditions:

u(R,t)=0.1, u(infinity,t)=0

and inital conditions:

u(x,0)=0

Thank you for your patience!

[Torsten]

Make two coordinate transformations:

Transforming the spherical diffusion equation by r' = r/R, t'=t gives

du/dt = 1/R^2 * 1/r^2 d/dr (r^2*D*du/dr) , u(r=1,t) = c0, u(Inf,t) = 0

Transforming this equation again by r' = 1/r, t' = t gives (please check)

du/dt = r^4/R^2 * D* d^2u/dr^2, u(0,t) = 0, u(1,t) = c0

or

R^2/r^4 * du/dt = D * d^2u/dr^2, u(0,t) = 0, u(1,t) = c0

This equation can now be solved in a cartesian frame (m=0). Instead of 0, you will have to use some small number (e.g. 1e-8) to start with at the left end.

A solution at r of the last equation corresponds to the solution at r' = R/r of the first (original) equation.

[Yannan]

Ok, Let me try for that! Thank you so much, I will let you know if I have futher question! :)

[Yannan]

@Torsten

Hello, Torsten! Could you please explain more how you do the coordinate transformation? I didn't really get the processes like why we set up r'=r/R and t'=t and we can transform the diffusion equation to the second one. Where does the 1/R^2 come from? Thank you!

[Torsten]

u(r,t) = v(r',t'), r' = r/R, t' = t ->
du/dt = dv/dr' * dr'/dt' + dv/dt' * dt'/dt = dv/dr' * 0 + dv/dt' * 1 = dv'/dt'
du/dr = dv/dr' * dr'/dr + dv/dt' * dt'/dr = dv/dr' * (1/R) + dv/dt' * 0 = 1/R * dv/dr'
d^2u/dr^2 = d/dr (du/dr) = d/dr (1/R * dv/dr') = 1/R * d^2v/dr'^2 * dr'/dr = 1/R * d^2v/dr'^2 * 1/R = 1/R^2 * d^2v/dr'^2
->
1/r^2 * d/dr(r^2*D*du/dr) = D*(2/r * du/dr + d^2u/dr^2) = D*(2/(r'*R) * 1/R * dv'/dr' + 1/R^2 * dv^2/dr'^2) =
D/R^2 * (2/r' dv/dr' + dv^2/dr'^2 ) = 1/R^2 * 1/r'^2 * d/dr' (r'^2*D*dv/dr')
->
dv/dt' = 1/R^2 * 1/r'^2 * d/dr' (r'^2*D*dv/dr')

A solution of the new PDE in(r',t')

dv/dt' = 1/R^2 * 1/r'^2 * d/dr' (r'^2*D*dv/dr') v(1,t') = c0, v(Inf,t) = 0

corresponds to a solution of the old PDE

du/dt = 1/r^2 * d/dr (r^2*D*du/dr) u(R,t) = c0, u(Inf,t) = 0

in (r=r'*R,t)

Now try the same with the new PDE and the transformation r' = 1/r, t' = t.

[Yannan]

@Torsten

Hi, Torsten! Thank you so much for the help and the method of coordinate transformation is so clever!

I've tried to apply the same thing on the transformation r'=1/r, t'=t. I found the r^4 should be 1/r^4, here is my process to make sure I didn't do something stupid:

[Torsten]

r' is the old spatial coordinate - you must express the new equation in the new coordinate "a" since c = c(a,t''):

dc/dt'' = D/R^2 * a^4 * d^2c/da^2

But besides this: well done.

[Yannan]

Ahh, I got it, thank you so much for the correction!! Instead of solving the complex m=2 PDE by using the pdepe, this time we only have m=0 in pdepe general format. So we only need to use the pdepe to calculate this with the new boundary conditions and plot the graph. That's genius!

[Torsten]

Yes, but as noted, the PDE has a singularity at a = 0. You must choose a small positive value for "a" at the left, e.g. aleft = 1e-6. This means you solve your original spherical PDE over the interval r = R up to r = R/a and set the boundary condition c = 0 at R/a.

[Yannan]

That is the last question I have! For pdepe we know, the boundary condition format is p(r,u,t)+q(r,t)f(r, t, u, dudr) =0 and we should define the pl,ql and pr, qr for the left and right boundary conditions. How we apply the initial condition and boundary conditions?

function u0 = diffic(x) %initial condition

u0 = 0;

end

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

pl = ul; % <-- how we define c(1e-6, t) = 0?

ql = 0;

pr = ur - 0.1;

qr = 0;

end

[Torsten]

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]

I understand! Thank you so much!!

[Yannan]

@Torsten

Hi Torsten! Sorry to keep distrubing you. I considered some further situations for this device. As I mentioned above, this device has a radius of R and the case we considered before is assuming the spherical device has infinity stable concentration C0 and the diffusion started at the surface of the device. But what if the total concentration is C0, as the diffusion occurs, the concentration C0 will decrease in the device with time as

u=C0*e^(-bt)

There is also a diffusion occurred in the device with diffusivity Di, and the tissue out side has diffusivity Dt. They share the same diffusion equation. Is this possible to get the concentration profile of u and x?

Thank you so much for the help.

[Torsten]

On the one hand you say that the concentration u in the device decreases by an exponential law, on the other hand you say that the diffusion equation in the device is valid with diffusivity Di.

Does that mean you want to compute two different cases: one where the concentration CR at r = R is given by CR=C0*e^(-bt), the other where you assume the diffusion equation in two regions: one from 0 to R and the other from R to Inf with different diffusivities and initial conditions C(0<=r<=R)=C0 and C(r>R) = 0 ?

[Yannan]

Yes! That's the exact thing I want to figure out!

It's just a more realistic situation for this device. As we know the device can carry C0 solution in total, and when it starts to release the solution into the tissue, the concentration will decrease as time passes. Therefore, the other diffusion appears in the device as the device concentration goes down.

I thought the general diffusion equation will not change but the initial condition becomes:

function u0 = diffic(x)

if x <= 1e-2 % the radius of device

u0 = 0.1;

else

u0 = 0;

end

end

And the diffusion equation becomes:

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

c = 1/x^4;

s = 0;

if x <= 1e-2

f = Di/L^2 * dudx; % Di=1e-13

else

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

end

end

and the boundary conditions will keep the same?

[Torsten]

It's not possible to work as you (we) did for the last problem where 1 was r = L and 0 was r = Inf because you only have one real boundary condition now, namely C = 0 at r = Inf.

So I suggest you work in the original physical frame [0,Inf] where you replace Inf by an adequate number XL >L and solve from 0 to L and L to XL with different diffusivities. Your function "diffpde" should work in this case, but with c = 1, f = Di(Dt) * dudx and s = 0. You will also have to choose m = 2 and the boundary conditions as dC/dr = 0 at r = 0 and C = 0 at r = XL. r = L should be a mesh coordinate in these simulations. The initial condition is as you set in your code above.

[Yannan]

I was consider whether we can set the left boundary condition very very small like r=1e-9, u=0 to transform the question to what we do before?

[Torsten]

But you don't have a boundary condition at r = L (resp. r = 1) as we had before.

We only have boundary conditions at r = 0 and r = Inf.

At r = L, we have two transmission conditions inside the simulation domain:

C(r->L-) = C(r->L+)
Di * dC/dr(r->L-) = Dt * dC/dr(r->L+)

[Yannan]

Yes, that's true. Ok, I got it! But return to the question I asked on the top of web. When I used m=2, the matlab will ignore the left boundary condition, if we want to set up the dc / dr = 0 at r = 0. How we can set up the boundary condition of pr, qr, pl and ql for the pdepe?

[Yannan]

@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
end

function u0 = diffic(x)                                % initial condition
      if x <= 1e-2
        u0 = 0.1;
    else
         u0 = 0;
     end
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]

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
end

function u0 = diffic(x,L)                                % initial condition
      if x <= L
        u0 = 0.1;
    else
         u0 = 0;
     end
end

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

[Yannan]

But for this case, the graph seems to turn 90 degrees straightly and doesn't look like the curve in tissue we got before any more. Is that because the diffusivity changed rapidly from the device to the tissue?

I'm also curious that if we want to draw the 3D graph with u, x, and t. How I set up the mesh to draw it?

Thank you for your patience again!

[Torsten]

Is that because the diffusivity changed rapidly from the device to the tissue?

Try it - set Di = Dt.

I'm also curious that if we want to draw the 3D graph with u, x, and t. How I set up the mesh to draw it?

Examples are provided under

https://uk.mathworks.com/help/matlab/ref/pdepe.html

[Yannan]

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
end

function u0 = diffic(x,L)                                % initial condition
      if x <= L
        u0 = 0.1;
    else
         u0 = 0;
     end
end

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

[Torsten]

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]

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
end

function u0 = diffic(x,L)                                % initial condition
      if x <= L
        u0 = 0.1;
    else
         u0 = 0;
     end
end

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