Steady-state solution for system of parabolic PDEs?

Is there an established method for finding the steady-state solution for a system of parabolic PDEs? I know how to use the Matlab ODE solvers to obtain transient-state solultions via Method of Lines (MoL), but am wondering if someone has come up with a reliable, robust strategy for going directly to a steady-state solution. As an example, think of a system of two PDEs describing the diffusion/reaction behavior of two different solutes along a 1-dimensional axis. I have tried using fsolve to obtain such a solution for where all the dydt's in the MoL are set equal to zero. This approach worked in some cases, but failed in others. Any suggestions will be much appreciated.

Answer 1

[Torsten]

Since this gives a system of second-order ODEs in general, the usual code to try would be "bvp4c".

[Eric Roden]

Dear Torsten, thanks for your reply. However, this is not a second-order ODE, but rather a PDE discretized in space and solved over time as a system of ODEs.

[Eric Roden]

Dear Torsen:

Wait, I see what you're saying! And yet, my skills in this type of coding are not really that well developed. I wonder if you could be so kind as to help me set the problem up for use with bvp4c. The problem I'm working on is 1-D spherical (simplified to radial) diffusion-reaction, which looks like the following:

where D is the diffusion coefficient, r is the radius, and the term is an enzyme kinetic function.

setting for steady state gives

The boundary conditions are fixed concentration at the surface of the spherical domain, i.e. at where R is the radial length (note that solute C diffuses into and is consumed within the sphere), and at the inner boundary.

Any assistance you may be able to provide with the coding will be much appreciated. I should mention here that I have been able to solve the time-dependent version of this problem using both pdepe as well a standard finite-difference approach implemented with ode15s. But I have never attemped to solve this kind of problem as a steady-state boundary value problem.

Kind regards, Eric

[Torsten]

CR = 1.0;
nu = 1.0e-2;
D = 1e-4;
K = 3;
%
rstart = 1e-8;
rend = 1.0;
rmesh = linspace(rstart,rend,100);
solinit = bvpinit(rmesh, @(r)guess(r,rend,CR));
%
sol = bvp4c(@(r,C)bvpfcn(r,C,nu,D,K), @(Ca,Cb)bcfcn(Ca,Cb,CR), solinit);
%
plot(sol.x, sol.y(1,:), '-o')

%
function dCdr = bvpfcn(r,C,nu,D,K)
  dCdr = [C(2);-2/r*C(2)+nu/D*C(1)/(K+C(1))];
end
%
function res = bcfcn(Ca,Cb,CR)
  res = [Ca(2);Cb(1)-CR];
end
%
function g = guess(r,R,CR)
  g = [0;0];
  if r==R
    g(1) = CR;
  end
end

[Eric Roden]

Dear Torsten:

Thanks ever so much for your help! I was able to implement the code no problem. Without meaning to abuse your kindness, I wonder if you could lead me one step further by showing me how to solve a system of two diffusion-reaction PDEs (see below). The first PDE is the same as before, and the second PDE represents another compound that is also consumed within the hypothetical sphere, where its consumption is inhibited by the presence of . The term represents an "inhibition function" that suppreses consumption of by .

Boundary conditions:

r = R,

r = 0,

Kind regards, Eric

[Torsten]

Exactly the same as the code above - only the vectors returned from bvpfcn, bcfcn and guess become 4x1 instead of 2x1.

Make an attempt. If you have problems, ask a specific question.

[Eric Roden]

I implemented the new model, and compared the results to a solution obtained with ode15s, which employs a well-tested central difference method of lines coding structure that I have used for many years. For some reason(s) I don't understand, the two approaches are producing substantially different results for C2. I think I have the bvp4c code correct, but perhaps there is an error? Note that I changed the boundary conditions and various parameters to match those for an actual problem I'm working on (diffusion-reaction of oxygen (C1) and nitrate (C2) in a hypothetical soil aggregate). Thanks very much for looking at this when/if time permits.

[Torsten]

Your discretization in "rgradfaqdunl" is wrong.

You don't use the Dirichlet-condition at r=R.

As you can see from your plots, you reverse the boundary conditions.

In the expressions

L(1) = D*cxx(1) + (2*D/r)*cx(1);

L(i) = D*cxx(i)+(2*D/r)*cx(i);

you use the fixed value r=1 in the division, but you must divide by x(i).

You seem to approximate the steady-state solution by using ODE15S artificially until steady-state is reached. "fsolve" made problems here ?

test1

The solution was obtained on a mesh of 37 points.
The maximum residual is  8.920e-04. 
There were 20571 calls to the ODE function. 
There were 325 calls to the BC function. 

test2

155 successful steps
17 failed attempts
2345 function evaluations
10 partial derivatives
50 LU decompositions
335 solutions of linear systems

function [] = test1
C1R = 4e-4;
C2R = 5e-4;
nu1 = 1.8e-4;
nu2 = 3.6e-5;
D = 0.054;
K1 = 1e-6;
K2 = 2e-4;
IC1 = 1e-6;
%
rstart = 1e-8;
rend = 1.0;
r = linspace(rstart,rend,100);
solinit = bvpinit(r, @(r)guess(r,rend,C1R,C2R));
%
options = bvpset('RelTol',1e-3,'AbsTol',1e-9,'Stats','on');
sol = bvp4c(@(r,C)bvpfcn(r,C,nu1,nu2,D,K1,K2,IC1), @(Ca,Cb)bcfcn(Ca,Cb,C1R,C2R), solinit, options);
%
x = sol.x;
C1 = sol.y(1,:);
C2 = sol.y(3,:);
figure
plot(x, C1, x, C2)
%set(gca,'xdir','reverse')
legend('C1','C2')
end
%
function dCdr = bvpfcn(r,C,nu1,nu2,D,K1,K2,IC1)
dCdr = [C(2); -2/r*C(2) + nu1/D*C(1)/(K1+C(1)); C(4); -2/r*C(4) + IC1/(IC1+C(1))*nu2/D*C(3)/(K2+C(3))];
end
%
function res = bcfcn(Ca,Cb,C1R,C2R)
res = [Ca(2);Cb(1)-C1R;Ca(4);Cb(3)-C2R];
end
%
function g = guess(r,R,C1R,C2R)
g = [0;0;0;0];
if r==R
    g(1) = C1R;
    g(3) = C2R;
end
end

function [] = test2
global nx npde neqn r dx x C1R C2R nu1 nu2 D K1 K2 IC1
r = 1;
nx = 100;
dx = r/nx;
x=linspace(dx,r,nx);
npde = 2;
neqn = npde*nx;
C1R = 4e-4;
C2R = 5e-4;
nu1 = 1.8e-4;
nu2 = 3.6e-5;
D = 0.054;
K1 = 1e-6;
K2 = 2e-4;
IC1 = 1e-6;

y0(1:npde:(neqn-(npde-1))) = C1R;
y0(2:npde:(neqn-(npde-2))) = C2R;

tspan=[0 1000];
options=odeset('RelTol',1e-3,'AbsTol',1e-9,'InitialStep',1e-6,'Stats','on','Events',@eventfun);
[t,y]=ode15s(@f,tspan,y0,options);
C1=(y(end,1:npde:(neqn-(npde-1))));
C2=(y(end,2:npde:(neqn-(npde-2))));
figure
plot(x, C1, x, C2)
legend('C1','C2')
end

function dydt = f(t,y)
global npde neqn C1R C2R nu1 nu2 D K1 K2 IC1
C1 = y(1:npde:(neqn-(npde-1)));
C2 = y(2:npde:(neqn-(npde-2)));
LC1 = rgradfaqdunl(C1R,C1,D);
LC2 = rgradfaqdunl(C2R,C2,D);
dydt(1:npde:(neqn-(npde-1))) = LC1 - (nu1.*C1./(K1+C1))';
dydt(2:npde:(neqn-(npde-2))) = LC2 - (IC1./(IC1+C1).*nu2.*C2./(K2+C2))';
dydt = dydt';
end

function [L] = rgradfaqdunl(c0,c,D)
% Compute radial transport flux for an aqueous component via central differences with Dirichlet
% upper and Neumann lower boundary conditions
global r nx dx
cx(1) = (c(2)-c0)/(2*dx);
cxx(1) = (c(2)-2*c(1)+c0)/dx^2;
L(1) = D*cxx(1) + (2*D/r)*cx(1);
for i=2:nx-1
    cx(i) = (c(i+1)-c(i-1))/(2*dx);
    cxx(i) = (c(i+1)-2*c(i)+c(i-1))/dx^2;
    L(i) = D*cxx(i)+(2*D/r)*cx(i);
end
cxx(nx) = (c(nx-1)-c(nx))/dx^2;
L(nx) = 2*D*cxx(nx);
end

function [x,isterm,dir] = eventfun(t,y)
dy = f(t,y);
x = norm(dy) - 1e-9;
isterm = 1;
dir = 0;
end

[Torsten]

Here is a correct discretization taken from

https://dl.acm.org/doi/pdf/10.1145/355644.355649

R = 1;
nr = 101;
r = linspace(0,R,nr).';
C1R = 4e-4;
C2R = 5e-4;
nu1 = 1.8e-4;
nu2 = 3.6e-5;
D = 0.054;
K1 = 1e-6;
K2 = 2e-4;
IC1 = 1e-6;

y0 = zeros(2*nr,1);
y0(nr) = C1R;
y0(2*nr) = C2R;

tspan=[0 1000];
[t,y]=ode15s(@(t,y)f(t,y,r,nr,D,K1,K2,IC1,nu1,nu2),tspan,y0);
C1 = y(end,1:nr);
C2 = y(end,nr+1:2*nr);
figure
plot(r, [C1; C2])
legend('C1','C2')

function dydt = f(t,y,r,nr,D,K1,K2,IC1,nu1,nu2)
  C1 = y(1:nr);
  C2 = y(nr+1:2*nr);
  dC1dt = zeros(nr,1);
  dC2dt = zeros(nr,1);
  dC1dt(1) = 3/((r(1)+r(2))/2)^3*((r(1)+r(2))/2)^2*D*(C1(2)-C1(1))/(r(2)-r(1))-...
             nu1*C1(1)/(K1+C1(1));
  dC2dt(1) = 3/((r(1)+r(2))/2)^3*((r(1)+r(2))/2)^2*D*(C2(2)-C2(1))/(r(2)-r(1))-...
             IC1/(IC1+C1(1))*nu2*C2(1)/(K2+C2(1));
  dC1dt(2:nr-1) = 3./(((r(2:nr-1)+r(3:nr))/2).^3 - ((r(1:nr-2)+r(2:nr-1))/2).^3).*...
                 (((r(2:nr-1)+r(3:nr))/2).^2 .* D .* (C1(3:nr)-C1(2:nr-1))./(r(3:nr)-r(2:nr-1))-...
                  ((r(1:nr-2)+r(2:nr-1))/2).^2 .* D .* (C1(2:nr-1)-C1(1:nr-2))./(r(2:nr-1)-r(1:nr-2)))-...
                  nu1*C1(2:nr-1)./(K1+C1(2:nr-1));
  dC2dt(2:nr-1) = 3./(((r(2:nr-1)+r(3:nr))/2).^3 - ((r(1:nr-2)+r(2:nr-1))/2).^3).*...
                  (((r(2:nr-1)+r(3:nr))/2).^2 .* D .* (C2(3:nr)-C2(2:nr-1))./(r(3:nr)-r(2:nr-1))-...
                  ((r(1:nr-2)+r(2:nr-1))/2).^2 .* D .* (C2(2:nr-1)-C2(1:nr-2))./(r(2:nr-1)-r(1:nr-2)))-...
                  IC1./(IC1+C1(2:nr-1)).*nu2.*C2(2:nr-1)./(K2+C2(2:nr-1));
  dydt = [dC1dt;dC2dt];
end

[Eric Roden]

Dear Torsen:

Thank you for continuing to work on this problem. I am not familiar with the differencing scheme adapted from the TOMS code (TOMS494). Is it some kind of higher-order differencing? If so I would like to learn more about this, although I found the Sincovec & Madson (1975) paper pretty difficult to follow. Meanwhile, as another check on things, I implemented the test problem in pdepe (code attached), and it produced results identical to my ode15s implementation. Not sure what this means, but if the actual constitutive equation for 1-D spherical diffusion-reaction system I'm using is incorrect, then of course pdepe will also produce incorrect results with the same input. For the record, the equation I have implemented is from a book called "Diagenetic Models and Their Implementation" by B.P. Boudreau (Springer), p. 63-65. A PDF of this book is available on ResearchGate. Any further thoughts you may have this will be appreciated.

Kind regards, Eric

[Eric Roden]

pdepe implementation

[Torsten]

If so I would like to learn more about this, although I found the Sincovec & Madson (1975) paper pretty difficult to follow.

It's simply the implementation of formula (3.1(c)) for r>0 and (3.2(c)) for r=0.

Meanwhile, as another check on things, I implemented the test problem in pdepe (code attached), and it produced results identical to my ode15s implementation.

I get the same results from pdepe as from the two other codes (bvp4c and ode15s).

You always set the Dirichlet boundary condition at r=0 instead of r=R. This is not possible in a spherical coordinate system.

global nx npde neqn r dx x C1R C2R nu1 nu2 D K1 K2 IC1
r = 1;
nx = 100;
dx = r/nx;
x=linspace(dx,r,nx);
npde = 2;
neqn = npde*nx;
C1R = 4e-4;
C2R = 5e-4;
nu1 = 1.8e-4;
nu2 = 3.6e-5;
D = 0.054;
K1 = 1e-6;
K2 = 2e-4;
IC1 = 1e-6;

tspan=(0:1000);
options = odeset('RelTol',1e-3,'AbsTol',1e-9,'NormControl','off','InitialStep',1e-7);
u = pdepe(2,@pdefun,@pdeic,@pdebc,x,tspan,options);
C1 = u(:,:,1);
C2 = u(:,:,2);

figure
plot(x,C1(end,:),x,C2(end,:));
ylim([0 5.5e-4])
legend('C1','C2')

function [c,f,s] = pdefun(x,t,u,dudx)
global r nu1 nu2 D K1 K2 IC1
c(1)=1;
c(2)=1;
f(1) = D*dudx(1);
f(2) = D*dudx(2);
s(1) = - nu1*u(1)/(K1+u(1));
s(2) = - IC1/(IC1+u(1))*nu2*u(2)/(K2+u(2));
c=c';
f=f';
s=s';
end

function u0 = pdeic(x)
global C1R C2R
u0(1) = C1R;
u0(2) = C2R;
u0=u0';
end

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
global C1R C2R
pr(1) = ur(1) - C1R;
pr(2) = ur(2) - C2R;
qr(1) = 0;
qr(2) = 0;
pl(1) = 0;
pl(2) = 0;
ql(1) = 1;
ql(2) = 1;
pl=pl';
ql=ql';
pr=pr';
qr=qr';
end

function [x,isterm,dir] = pdevents(m,t,xmesh,umesh)
du = fnss(t,y);
x = norm(dy) - 1e-9;
isterm = 1;
dir = 0;
end

Answer 2

[Eric Roden]

Dear Torsten:

I finally got it! I obviously did not understand properly the spherical diffusion equation, and was thus not implementing the correct central difference approximations in the solvers. I wound-up going with the specific finite difference expressions given on p. 320 of the Boudreau (1996) book as opposed to the general ones in Sincovec and Madsen (1975), but these produced results identical to those obtained with the latter. And the results from the three different solvers (ode15s, pdepe, and bvp4c) are consistent. I am very grateful to you for helping me on my way.

Kind regards, Eric

[Torsten]

My pleasure.