Error in solving system of two reaction-diffusion equations?

Hello there,

I have a system of two reaction-diffusion equations that I want to solve numerically (attached is the file). However, it doesn't resemble with the standard system used in pdepe.m example. The problems I have are:

(1) I don't know how to incorporate it and write c, f, s for my system. As per my knowledge the problem is with the extra term $\gamma \partial u/\partial x$.

(2) I tried coding it up but I got an error:

Warning: Failure at t=2.806210e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.551115e-17) at time t. > In ode15s (line 730) In pdepe (line 289) In pde_bee (line 67)

Here is the code I wrote:

   m = 0; % not sure what m should be for my system of equations
  x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
  t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure;
surf(x,t,u1);
title('u1(x,t)');
xlabel('Distance x');
ylabel('Time t');

figure;
surf(x,t,u2);
title('u2(x,t)');
xlabel('Distance x');
ylabel('Time t');

% --------------------------------------------------------------------------

function [c,f,s] = pdex4pde(x,t,u,DuDx)
delta = 0.5;
gamma = 0.1;
alpha = 0.1;
beta = 0.1;
muB = 0.2;
kB = 0.2;
deltaB = 0.1;
Du1Dx = DuDx(1);
%Du2Dx = DuDx(2);
u1 = u(1); % u in my eqns
u2 = u(2); % v in my eqns
c = [1; 1];
f = [delta*Du1Dx-gamma*u1; 0] .* DuDx;
F = -alpha*u1 + beta*u2;
s = [F; muB*kB*u2^2/(kB^2+u2^2)-deltaB*u2-F];

% --------------------------------------------------------------------------

function u0 = pdex4ic(x)
u0 = [1; 0];

% --------------------------------------------------------------------------

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

Thank you so much for your attention.

Answer 1

[Torsten]

"pdepe" is not designed to solve mixtures of partial differential equations and ordinary differential equations (the equation for u is a partial differential equation, the equation for v is an ordinary differential equation).

You will have to discretize the first equation in space on your own and solve the complete resulting system of ordinary differential equations using ODE15S.

Look up "method of lines" for more details.

Best wishes

Torsten.

[Rose]

@Torston Thank you very much! Is there a sample code that you would suggest? I found material on the analytical part (i.e. discretization) but didn't find any sample code that would be similar to my problem.

[Rose]

I also saw another similar question where you suggested NAG toolbox. Is that something that I could also use? Here is the question.

[Torsten]

If you have access to the NAG toolbox, you can use the cited function (d03plf) for your problem because it is able to solve coupled ordinary differential equations.

But discretizing your first equation using the method-of-lines is not that hard - I'd give it a try.

Best wishes

Torsten.

[Rose]

I don't have access to the toolbox. So I discretized the system. I haven't yet implemented it using ode15s. But I believe it will be primarily a for loop for the grid points.

[Torsten]

Try this:

x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

n = numel(x);
u0 = ones(n,1);
v0 = zeros(n,1);
y0 = [u0 ; v0];

delta = 0.5;
gamma = 0.1;
alpha = 0.1;
beta = 0.1;
muB = 0.2;
kB = 0.2;
deltaB = 0.1;

[T Y] = ode15s(@(t,y)fun(t,y,x,n,delta,gamma,alpha,beta,muB,kB,deltaB),t,y0);

function DyDt=fun(t,y,x,n,delta,gamma,alpha,beta,muB,kB,deltaB)

u = zeros(n,1);
v = zeros(n,1);
DuDt = zeros(n,1);
DvDt = zeros(n,1);
DyDt = zeros(2*n,1);

xhalf = zeros(n-1,1);
DuDx = zeros(n-1,1);
D2uDx2 = zeros(n-1,1);

u = y(1:n);
v = y(n+1:2*n);

xhalf(1:n-1)=(x(1:n-1)+x(2:n))/2;

DuDx(1) = 0.0;
D2uDx2(1) = 1.0/(xhalf(1)-x(1))*(u(2)-u(1))/(x(2)-x(1));
for i=2:n-1
  DuDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(u(i+1)-u(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(u(i)-u(i-1)))/(x(i+1)-x(i-1));
  D2uDx2(i) = ((u(i+1)-u(i))/(x(i+1)-x(i))-(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
end

for i=1:n-1
  F = -alpha*u(i)+beta*v(i);
  DuDt(i) = -gamma*DuDx(i)+delta*D2uDx2(i)+F;
  DvDt(i) = v(i)*(muB*kB*v(i)/(kB^2+v(i)^2)-deltaB)-F;
end
F = -alpha*u(n)+beta*v(n);
DuDt(n) = 0.0;
DvDt(n) = v(n)*(muB*kB*v(n)/(kB^2+v(n)^2)-deltaB)-F;

DyDt = [DuDt;DvDt];

end

Best wishes

Torsten.

[Jason Coumans]

Hello Torsten.

Out of curiosity which Finite Difference scheme is this? It looks like a central difference but I can't quite figure it out.

Thank you! -J

[Torsten]

The discretization schemes are described in

https://opus4.kobv.de/opus4-zib/files/505/TR-93-14.pdf

Page 9-11 explains the discretization of the first spatial derivative.

Best wishes

Torsten.

[SDO]

[ZIYI LIU]

Hello Torsten,

What if the boundary conditions of u contain dvdt?

I have a really similar system to this one, but my boundary conditions of u is dvdt, and I wrote my code but cannot get the correct result.

Thanks!

Best,

Ziyi

Answer 2

[Bill Greene]

The flux term in your pdex4pde function doesn't match the equations you show in your attached pdf. Assuming the pdf is correct, change f to this:

f = [delta*Du1Dx-gamma*u1; 0];

That change will also resolve the convergence problems you are observing.

[Torsten]

Although I don't recommend to use "pdepe" for the problem because of the above mentionned reason, your setting of the flux term is also incorrect because at x=0, Rose wants to set du/dx = 0, not delta*du/dx-gamma*u=0.

To account for this, one has to put the convective part of the flux into the source term:

c=[1;1];
f=[delta;0].*DuDx;
s=[-gamma*DuDx(1)+F;muB*kB*u2^2/(kB^2+u2^2)-deltaB*u2-F];

u(2) does not need boundary conditions; so setting pl(2), ql(2), pr(2) and qr(2) is somewhat arbitrary and (in my opinion) will most probably lead to erroneous results.

Best wishes

Torsten.

[Bill Greene]

I could not find any specification of boundary conditions in the problem description so it is unclear to me exactly what the desired boundary conditions are. It is not uncommon to set the total flux (both diffusive and convective components) equal zero. So that seems just as reasonable to me as setting du/dx=0.

Regarding the setting of boundary conditions for u(2), since the flux is identically zero in equation 2, setting qr(2)=ql(2)=1 and pr(2)=pl(2)=0 has no effect, adverse or otherwise, on the solution. Its just a simple way to satisfy pdepe's requirement for boundary conditions on all dependent variables.

In general, pdepe is far more versatile than many people (including Mathworks documentation writers) understand.

[Torsten]

From a physical standpoint, setting the total flux to zero at x=0 makes no sense because the convective flux of u must leave the system. Otherwise, u will accumulate at x=0.

It's correct that setting qr(2)=ql(2)=1 and pr(2)=pl(2)=0 has no effect on the solution for v if v was the only solution variable. Nonetheless, the boundary values that v attains over time will influence u. And this not only at the boundary, but over the complete domain because of the fluxes of u in the spatial coordinate. And since u and v are coupled, this in turn will also influence v over the complete domain.

I admit that I don't know about the strength of this effect, but the uncertainty alone would keep me from using this tool for the above problem.

Best wishes

Torsten.