Steady-state solution for system of parabolic PDEs?

7 vues (au cours des 30 derniers jours)
Eric Roden
Eric Roden le 2 Nov 2022
Commenté : Torsten le 5 Nov 2022
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.

Réponse acceptée

Torsten
Torsten le 3 Nov 2022
Since this gives a system of second-order ODEs in general, the usual code to try would be "bvp4c".
  11 commentaires
Eric Roden
Eric Roden le 5 Nov 2022
pdepe implementation
Torsten
Torsten le 5 Nov 2022
Modifié(e) : Torsten le 5 Nov 2022
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

Connectez-vous pour commenter.

Plus de réponses (1)

Eric Roden
Eric Roden le 5 Nov 2022
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

Produits


Version

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by