Solving partial differential system

2019 3 24
1 个回答
回答已采纳
更新时间：2019 10 11
8 次查看（30 天）

0 个投票

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.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Torsten
Torsten 2019-3-25
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
Daniel Head 2019-3-25
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 2019-3-25
Insert the expression for dq/dt from the second equation.
Daniel Head
Daniel Head 2019-3-25
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
m_p 2019-10-11
can u send me your final correct code i am facing some problem ?

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2019-3-25
编辑：Torsten 2019-3-25

1 个投票

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

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Daniel Head
Daniel Head 2019-3-25
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
Torsten 2019-3-25
编辑：Torsten 2019-3-25
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
Daniel Head 2019-3-25
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
Torsten 2019-3-25
编辑：Torsten 2019-3-25
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
Daniel Head 2019-3-25
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 2019-3-25
In pdex4pde, you overwrite the two-element solution variable vector u by the setting
u=Q/A_c;
Daniel Head
Daniel Head 2019-3-25
I have changed this variable name and now the system solves! Thank you!

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 PDE Solvers 的更多信息

产品

版本

R2018b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by