Solving a system of PDE using pdepe

2 views (last 30 days)
nir livne on 30 Jun 2022
Commented: Bill Greene on 1 Jul 2022
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
Bill Greene on 1 Jul 2022
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?

Torsten on 30 Jun 2022
Edited: Torsten on 30 Jun 2022
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
nir livne on 1 Jul 2022
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!

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by