Solving a system of PDE using pdepe

2 views (last 30 days)
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 Comment
Bill Greene
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?

Sign in to comment.

Accepted Answer

Torsten
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
  5 Comments
nir livne
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!

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by