Solving a system of PDE using pdepe
28 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commento
Bill Greene
il 1 Lug 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?
Risposta accettata
Torsten
il 30 Giu 2022
Modificato: Torsten
il 30 Giu 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 Commenti
Torsten
il 1 Lug 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.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su PDE Solvers in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!