Time-dependent point source for diffusion with pdepe

Hi Everyone,

I was hoping someone could give me a hand with this diffusion problem. In brief, I'm trying to simulate two-dimensional diffusion which includes a time-dependent source function to reproduce some experimental results, however I've encountered challenges...

I have no problem simulating diffusion alone, however the problem arises in introducing the source for t>0 at x = 0. Currently I've been trying to do the following:

 function [c,f,s] = pdex4pde(x,t,u,DuDx, xSpan, tSpan)
 c = [1; 1];
 f = [D1; D2] .* DuDx;                           % diffusion term
 k1 = 0.02;                                      % degradation rate for species 1
 k2 = k1/4;                                      % degradation rate for species 2
 sourcePosition = 0;       %micrometers
 releaseTime = 1;          %seconds
 initialConc = 0.001
  if x == sourcePosition  & t <= releaseTime ;
        n = (0.221 * releaseTime ) + 1.662;     % source function parameter
        k = 0.138 * releaseTime  ^ 1.53);      % source function parameter

        releaseFunction = (initialConc .* (tSpan.^n) ./ ((k.^n) + (tSpan.^n)));                   % sourcefunction
        releaseRates = diff(releaseFunction);                                                     % source function differential

        if t == 0                              % at t = 0, there is nothing in the system
            Source=0;
        else                                   % for t>0, diffusing species is introduced into the system
            Source = interp1(tSpan(2:end), releaseRates, t, 'linear', 'extrap');
        end
   else
        Source= 0;
   end
F1 = -k1*u(1) + Source;
F2 = k1*u(1)-k2*u(2);
s = [F1; F2];
end

Any suggestions would be immensely appreciated, I've been struggling with this for a while now. I'm do not have a strong math background, so the simpler the better if possible. If there's any other information what is needed, please let me know.

Thanks in advance!

[Torsten]

Is the source position a boundary point of your domain of integration ?

[Lukas]

Switches are generally hard to simulate. I would either use events or separate calls of the ode or pde solver. Your initial value of Source is only zero at exactly t=0. This will not have any effect on the simulation result, as the timespan in inf. small. Better omit this statement.

[Nicholas Mikolajewicz]

Torsten, yes my boundaries are 0,x and my source is at x = 0.

currently I have my boundaries set as follows:

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

I've also tried something like this with no luck:

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [ul -pdexSource(t); 0];
ql = [0;1];
pr = [0; 0];
qr = [1;1];
 end

Lukas, thanks I'll omit the source = 0 at t= 0 statement. Also, would the ode be called directly into my pdex4pde(x,t,u,DuDx, xSpan, tSpan) function? I currently call my source function into pdex4bc, but the way that pdepe is discretizing the time steps isn't allowing me to accurately introduce the diffusing species at a defined rate (perhaps I'm missing something here?). I've also looked into events, but to be frank I don't really understand how to implement them in this scenario.

Nick

[Torsten]

If the source is at the boundary, it has to be incorporated in the boundary condition function pdex4bc and removed from the pdex4pde function.

If you know the value u1* of the source at time t, maybe a boundary condition of the form

-D1*du1/dx = ka*(u1-u1*)

(ka = ka-value) is adequate.

Knowing more about the process you are trying to describe might help.

Best wishes

Torsten.

[Nicholas Mikolajewicz]

Ok great, would you mind guiding me through this a bit more. The source is intended to model the transient release of a signalling mediator from a cell. The intracellular concentration c(t=0) = 0.001 is released into the extracellular space at a rate of dc/dt = P(t)*c(t-1) where P(t) is the time-dependent release rate(per second, s^-1). In the model, I have pdexSource representing dc/dt at x=0, and currently the rate of release from the intracellular space is dependent on the intracellular concentration, not the extracellular concentration(i.e., u1)

Setting up the boundary with these conditions is what confuses me. From my understanding, the boundary condition you proposed -D1*du1/dx = ka*(u1-u1*) is dependent on the amount already released. But since the release rate is dependent on the amount within the cell that hasn't been released yet, how would I set up the boundaries to reflect that?

This is my guess but I feel it's wrong in some why...

D1*du1/dx = pdexSource

pl = -pdexSource
ql = 1
pr = 0
ql = 1

thanks again,

Nick

[Torsten]

I don't understand the physical background - so I can't give you further advice. But your solution with

pl = -pdexSource;
ql = 1.0;

can't be correct because D1*du1/dx has unit mol/(m^2*s) wheras pdexSource has unit mol/(m^3*s), I guess.

Best wishes

Torsten.

[Nicholas Mikolajewicz]

Hi Torsten,

what other information would help understand the problem?

Here are all the details I can provide: the source is positioned at x=0, at the center of a cylinder. The radius of the cylinder is 500, such that x=0:500. There is no reflection at the boundaries of the cylinder. As time t tends to infinity, the concentration c tends to 0. The amount of diffusing species introduced at x=0 is time dependent, and decreases exponentially with time.

my biggest problems are that (i) I don't know how to introduce the source properly and (ii) I need to compare the solution concentrations to those I've measured experimentally, however the time and distance steps influence the outcome concentrations. (iii) the concentrations I introduce are in terms of moles per volume...I'm not sure how this affects interpretation of the solution if it assumes only radial diffusion?

With respect to my first problem, I think I found a work around where I introduce the source as an initial condition, simulate the problem for one time step and use the corresponding solution after one time step as the initial condition in a subsequent simulation. In each subsequent simulation, in addition to the solution from the previous simulation I also add the time-dependent source to obtain a new initial condition, and so on. This approach has given me the the diffusion profiles that I've expected, however the absolute concentrations are way off and I have no idea how to correct this.

If I'm asking for too much, please let me know so I stop bothering you, however help would be greatly appreciated as I'm completely lost on this right now.

Thank you,

Nick