Error in solving system of two reaction-diffusion equations?

64 views (last 30 days)

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.

Accepted Answer

Torsten
Torsten on 6 Dec 2017
"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.
  9 Comments
ZIYI LIU
ZIYI LIU on 8 Oct 2020
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

Sign in to comment.

More Answers (1)

Bill Greene
Bill Greene on 14 Dec 2017
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.
  3 Comments
Bill Greene
Bill Greene on 15 Dec 2017
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
Torsten on 15 Dec 2017
Edited: Torsten on 15 Dec 2017
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.

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!