Solving partial differential system

Hello,

I am trying to sovle a modelling problem, however pdepe will not solve a mixed PDE system (as far as I know) with the following PDE's:

With the following initial conditions:

And boundary conditions:

A similar question was asked here: https://uk.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations this talks about the use of discretisation, something I am unsure of and I am unable to figure out how to use this code (suggested by Torsten) with my system.

Any help or codes would be greatly appreciated.

[Torsten]

After a long discussion with Bill Greene, we finally agreed that pdepe can be used to solve the system from above although the second equation does not include spatial derivatives.

As boundary condition, we suggest using dq/dx = 0 at both boundaries.

[Daniel Head]

How would I go about splitting the first PDE into the c,f,s format? I have tried, but I intially end up with c being non diagonal, by diagnolising this matrix I essentially lose the final term of the equation.

[Torsten]

Insert the expression for dq/dt from the second equation.

[Daniel Head]

Here is my code that I'm using, the main issue I have is I want to define q as u2 but I am unsure how to define it:

function Modelling1

m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);


figure;
surf(x,t,u1);
title('u1(x,t)');
xlabel('Distance x');
ylabel('Time t');

figure;
surf(x,t,u2);
title('u2(x,t)');
xlabel('Distance x');
ylabel('Time t');

% --------------------------------------------------------------------------

function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1];
Q=2.5; %mL min^-1
A_c=1.96e-5; %cm^2
u=Q/A_c;
e_b=0.368;
v=u/e_b;
A=35.13;
d_p=50e-6; %m
D_L=A*0.5*d_p*v;
k_max=0.18; %min^-1
S_1=0.6245;
S_2=2.071;
q_sat=94.72; %mg mL^-1
H_ref=246.8;
pH=7.5;
pH_ref=6;
N=0.5;
H_c=H_ref*(pH/pH_ref)^N
qst=H_c/(1+(H_c/q_sat));
qe=qst;
km=k_max*(S_1+(1-S_1)*(1-(qst/q_sat))^(S_2^2));
f = [D_L; 0] .* DuDx;
s = [-v*DuDx; km*(qe-u2)];

% --------------------------------------------------------------------------

function u0 = pdex4ic(x)
u10=1.2; %mg mL^-1
u20=0; %mg mL^-1
u0 = [u10; u20];

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [u1-pu2; ul(2)];
ql = [-1/v; 0];
pr = [1/DL; 0];
qr = [0; 1];

[m_p]

can u send me your final correct code i am facing some problem ?

Answer 1

[Torsten]

s= [-v*DuDx(1)-psi*km*(qe-u(2)); km*(qe-u(2))];

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
Q=2.5; %mL min^-1
A_c=1.96e-5; %cm^2
u=Q/A_c;
e_b=0.368;
v=u/e_b;
pu2 = ?;

pl = [ul(1)-pu2; 0.0];
ql = [-1/v; 1.0];
pr = [0; 0];
qr = [1.0; 1.0];

[Daniel Head]

For pu2, what am I definining here? If I set this to a number I am given the error not enough input arguments with regards to the use of psi. Sorry about the multitude of questions, I greatly apreciate your help.

[Torsten]

In pdex4pde, set psi to a number and replace s by the vector I wrote.

In pdex4bc, set pu2 to a number and use the function as written.

[Daniel Head]

I am assuming psi is Φ and defined it as such (I had forgotten it initally sorry).

Within pde:

e_p=50e-6;
psi=1-e_p/e_p;
s= [-v*DuDx(1)-psi*km*(qe-u(2)); km*(qe-u(2))];

Within bc:

e_p=50e-6;
psi=1-e_p/e_p;
pl = [ul(1)-psi*u2(2); 0.0];

Now I get a new error: Index exceeds the number of array elements (1). This is with regards to defining s.

[Torsten]

Do you see any array u2 in pdex4bc ? I don't. pu2 is your c_in in the set of equations from above - you'll have to set it to a fixed value.

And note that psi=0 since you set psi=1-e_p/e_p.

[Daniel Head]

I've changed pu2 to my c_in (I misunderstood what it was representing) and I've changed psi to be psi=(1-e_b)/e_b, i.e. not equalling 0 and I had used the wrong value, but still get the same error with regards to defining s within pdex4pde.

[Torsten]

In pdex4pde, you overwrite the two-element solution variable vector u by the setting

u=Q/A_c;

[Daniel Head]

I have changed this variable name and now the system solves! Thank you!