## Solving partial differential system

on 24 Mar 2019 at 20:19
Latest activity Commented on by Daniel Head

on 25 Mar 2019 at 12:30

### Torsten (view profile)

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.

Show 1 older comment

on 25 Mar 2019 at 9:52
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

### Torsten (view profile)

on 25 Mar 2019 at 9:59
Insert the expression for dq/dt from the second equation.

on 25 Mar 2019 at 10:23
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];

R2018b

### Torsten (view profile)

on 25 Mar 2019 at 10:36
Edited by Torsten

### Torsten (view profile)

on 25 Mar 2019 at 10:52

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];

on 25 Mar 2019 at 12:09
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

### Torsten (view profile)

on 25 Mar 2019 at 12:23
In pdex4pde, you overwrite the two-element solution variable vector u by the setting
u=Q/A_c;