Solving a system of PDE using pdepe

Question

nir livne 2022-6-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1751260-solving-a-system-of-pde-using-pdepe

评论： Bill Greene 2022-7-1

采纳的回答： Torsten

在 MATLAB Online 中打开

Hi,

I'm trying to solve system of 2 PDE's. It is a one-dimensional problem (cylinderical coordinates with symmetry):

with the following boundary conditions:

,

(R stands for r=R which is the boundary of the domain).

I'm keep getting index errors but I can't figure out why, I've been stuck for a while now... Specifically, the 'pdefun' is keeping me stuck right now with the following error:

Index exceeds the number of array elements. Index must not exceed 1.

Error in AggSim>pdefun (line 35)

f = [Dn, alpha*(u(1)/u(2)); 0, Dc]* dudx;

Any help would be appreciated! Thanks in advnace :)

%% constants and space/time variables
L = 0.5;
dL = 0.001;
x = 0:dL:L;
t_steps = 100;
t_f = 1;
t = linspace(0, t_f, t_steps);
m = 1;
alpha = 10^-3;
Dn = 4 * 10^-6;
Dc = 9 * 10^-6;
k = 10^-10;
pH0 = 5.5;
beta = 0.1 * 10^-(pH0);
%% solve and plot
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')
%% function defs
function u0 = icfun(x)
global pH0
u0 = [1, 10^-pH0];  
end
function [c,f,s] = pdefun(x, t, u,dudx)
global alpha Dn Dc k beta
c = [1; 1];
f = [Dn,  alpha*(u(1)/u(2)); 0, Dc]* dudx;
s = [0; beta -k*u(1)*u(2)];
end
function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
global pH0
pL = [1, 1];
qL = [0; 0];
pR = [1; 0];
qR = [0; uL(2)-10^-pH0];
end

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

Bill Greene 2022-7-1

In your equations for the boundary conditions you show some derivatives with respect to t. Are those really supposed to be derivatives with respect to r?

请先登录，再进行评论。

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

Answer 1

Torsten 2022-6-30

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1751260-solving-a-system-of-pde-using-pdepe#answer_997465

编辑：Torsten 2022-6-30

在 MATLAB Online 中打开

The error is solved, but I think the boundary conditions at the right end can't be set within pdepe.

The condition set at the moment (by me) is not what you want.

I assumed beta = c0 in your code.

%% constants and space/time variables

global alpha Dn Dc k beta

global pH0

L = 0.5;

dL = 0.001;

x = 0:dL:L;

t_steps = 100;

t_f = 10000;

t = linspace(0, t_f, t_steps);

m = 1;

alpha = 10^-3;

Dn = 4 * 10^-6;

Dc = 9 * 10^-6;

k = 10^-10;

pH0 = 5.5;

beta = 0.1 * 10^-(pH0);

%% solve and plot

sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);

u1 = sol(:,:,1);

u2 = sol(:,:,2);

surf(x,t,u1)

title('u_1(x,t)')

xlabel('Distance x')

ylabel('Time t')

%% function defs

function u0 = icfun(x)

global pH0

u0 = [1; 10^-pH0];

end

function [c,f,s] = pdefun(x, t, u,dudx)

global alpha Dn Dc k beta

c = [1; 1];

f = [Dn, alpha*(u(1)/u(2)); 0, Dc]*dudx;

s = [0; -k*u(1)*(u(2)-beta)];

end

function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)

global pH0

pL = [0; 0];

qL = [1; 1];

pR = [0; uR(2)-10^(-pH0)];

qR = [1; 0] ;

end

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2022-7-1

You forgot to include the globals in the script part of your code.
s in pdefun has changed. I assumed beta = c0 and thus set s(2) = -k*u(1)*(u(2)-beta).
The boundary condition setting (pL qL, pR, qR) is substantially different from your settings. At the moment, the setting at r=R is c = 10^(-pH0) and D_n*dn/dr - alpha*n/c * dc/dr = 0. I guess that with your definition of f in pdefun, it is impossible to set dn/dr = 0 at r=R in pdepe.

nir livne 2022-7-1

Actually, your setting of D_n*dn/dr - alpha * n/c * dc/dr is what I meant to set. Sorry for the confusion I've created! I meant the flux of n to vanish at r=R.

You've been exteremly helpful! Thank you!

请先登录，再进行评论。