Solving a system of PDE using pdepe

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]

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]

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

[SARSKOLIN FOSSO]

Thank you Mr. Torsten for your help in this community. Please help me to solve numerically this coupling of PDE and ODE with the conditions given in the pdf named coupling1. I have tried pdepe but no results found.

[Torsten]

@SARSKOLIN FOSSO

Please open a new question since this seems to have nothing to do with the question asked by @nir livne

I will delete your comment here in the course of the day.

[nir livne]

Thanks a lot for the answer! I'll fix the boundary conditions (still trying to wrap my head around how to set those in pdepe).

I'm reading your lines of code and I can't find what went wrong with mine - what changes did you make?

By the way, I left it out in my original question, but beta is supposed to be a constant addition to c(x,t)

Again, thanks a ton for your help!

[Torsten]

1. You forgot to include the globals in the script part of your code.
2. s in pdefun has changed. I assumed beta = c0 and thus set s(2) = -k*u(1)*(u(2)-beta).
3. 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]

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!