Solving a system of PDE using pdepe
30 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
nir livne
el 30 de Jun. de 2022
Comentada: Bill Greene
el 1 de Jul. de 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
1 comentario
Bill Greene
el 1 de Jul. de 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?
Respuesta aceptada
Torsten
el 30 de Jun. de 2022
Editada: Torsten
el 30 de Jun. de 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 comentarios
Torsten
el 1 de Jul. de 2022
- 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.
Más respuestas (0)
Ver también
Categorías
Más información sobre PDE Solvers en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!