Steady-state solution for system of parabolic PDEs?

Question

Eric Roden 2022-11-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1842313-steady-state-solution-for-system-of-parabolic-pdes

评论： Torsten 2022-11-5

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-11-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1842313-steady-state-solution-for-system-of-parabolic-pdes#answer_1090173

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

11 个评论
显示 9更早的评论隐藏 9更早的评论

Eric Roden 2022-11-3

编辑：Eric Roden 2022-11-3

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 2022-11-4

编辑：Torsten 2022-11-4

在 MATLAB Online 中打开

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

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 2022-11-4

编辑：Torsten 2022-11-4

在 MATLAB Online 中打开

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 2022-11-5

test_pdepe.m

pdepe implementation

Torsten 2022-11-5

编辑：Torsten 2022-11-5

在 MATLAB Online 中打开

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 2022-11-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1842313-steady-state-solution-for-system-of-parabolic-pdes#answer_1091963

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2022-11-5

My pleasure.

请先登录，再进行评论。

Steady-state solution for system of parabolic PDEs?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Steady-state solution for system of parabolic PDEs?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论