Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

I get this error in writing the code.... The code I wrote is as follows

function pdeop global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a; m = 0; x = linspace(0,1.5,31); t = linspace(0,3600,3601); sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3=sol(:,:,3); u4=sol(:,:,4);

% -------------------------------------------------------------- function [c,f,s] = pdeoppde(x,t,u,DuDx) c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol]; f = [-G;L;-G;L] .*u; Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2); ytstar=0.0023*exp(6.0058*u(4)); M=Cs*(ysstar-u(1)); N=Ct*(ytstar-u(3)); s = [M;-M;N;-N]; end % -------------------------------------------------------------- function u0 = pdeopic(x) u0 = [0;0.5;0;0.5]; end % -------------------------------------------------------------- function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t) pl = [ul(1);0;ul(3);0]; ql = [0;0;0;0]; pr = [ur(1)-0.5;0;ur(3)-0.5;0]; qr = [0;0;0;0]; end end

Can anyone please suggest me a solution?

[Alan Weiss]

Please edit this and use the {} Code button so we might be able to read it.

Alan Weiss

MATLAB mathematical toolbox documentation

[Torsten]

... and best include the PDE system with boundary and initial conditions in mathematical notation as pdf document such that we are able to compare it to your code.

Best wishes

Torsten.

[Widitha Samrakoon]

My code contains two functions named paracal.m which calculated parameters to be used in the pde system and then the main function named pdeop.m. First paracal.m is run and then pdeop.m is run. First I will paste the code of paracal.m and then pdeop.m. The error occurs when running pdeop.m.

function paracalc
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
A=pi*0.0254^2*0.25;
a=1710;
L=7.2888*10^-5;
G=2.4296*10^-3;
Rhog=842.77;
Rhol=875;
void=0.86;
Usg=G/(Rhog*A);
Usl=L/(Rhol*A);
hl=(1.5688*Usl^(1/3))*(1+4.1435*Usg^0.05)*(1-exp(-88.9560*Usl^(0.4)))^(2/3);
Kos=(1278102.274+6643.8785*Usg^-0.5)^-1;
Kot=(1620456.645+10029.7854*Usg^-0.5)^-1;
aeff=a*(1-exp(-88.9560*Usl^0.4));
Ugeff=1.6444*Usg/(1-hl);
Uleff=1.6444*Usl/hl;
end

Now I will paste the code of pdeop.m which gives the error

function pdeop
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
m = 0;
x = linspace(0,1.5,31);
t = linspace(0,3600,61);
sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3=sol(:,:,3);
u4=sol(:,:,4);
function [c,f,s] = pdeoppde(x,t,u,DuDx)
c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol];
f =[-G;L;-G;L].*u+[-0.001;0.001;-0.001;0.001].*DuDx;
Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2);
ytstar=0.0023*exp(6.0058*u(4));
M=Cs*(ysstar-u(1));
N=Ct*(ytstar-u(3));
s = [M;-M;N;-N];
function u0 = pdeopic(x)
u0 = [0;0.5;0;0.5]; 
end
function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t)
pl = [ul(1);0;ul(3);0]; 
ql = [0;0;0;0]; 
pr = [0;ur(2)-0.5;0;ur(4)-0.5]; 
qr = [0;0;0;0]; 
end
end

Please note that my system of PDEs has four elliptical equations but the problem is all four of them do not have second order spatial derivatives. Therefore in my system, the flux vector doesn't have a DuDx term but U term only. By a search in the internet I got to know that pdepe function cannot work with such equations.However in order to get rid of this, as suggested in an answer to a similar question, I included a second order spatial derivative to all four equations which you can see as a separate column vector in the flux vector as [-0.001;0.001,-0.001,0.001]. I hoped this would solve the problem. But the error message keeps popping out. I am really helpless. Please help.

I am attaching a pdf file of my PDE system with the initial conditions and boundary conditions.

Answer 1

[Bill Greene]

I see two problems.

First, your boundary conditions are incorrectly defined. They should be:

pl = [ul(1);ul(2);ul(3);ul(4)]; 
pr = [ur(1);ur(2)-0.5;ur(3);ur(4)-0.5];

(ql and qr are correct). All four of your equations are hyperbolic which pdepe is not designed to handle. But, you have added some artificial diffusion (the small second derivative terms). That is a good approach but because of the negative signs on some of the coefficients you have de-stabilized the system instead stabilizing it. Try this for the vector of diffusion coefficients, instead:

eps=.001;
[eps eps eps eps]'

[Widitha Samrakoon]

Thank you Bill for the hint. However, this system only has four know boundary conditions ( which should be enough to solve since the order of spatial derivatives is one in each). ul(2), ul(4) and ur(1) and ur(3) are not known. that's why I put 0 for both p and q in place of those variable. The objective of solving this system is finding the variation of those values with time....However, I will adjust my diffusion derivative and check and reach you back. Thank you

[Bill Greene]

OK, I took a guess at what your intention was with respect to the boundary conditions. pdepe requires boundary conditions at both ends and p=q=0 is not acceptable. I suggest p=0 and q=1 at the ends where you really don't want a BC. This will set your "artificial" flux to zero there but have minimal impact on your "real" PDE.

Answer 2

[Torsten]

As Bill already pointed out, pdepe is not designed to solve PDEs without 2nd derivative terms.

You will have to discretize the spatial first derivatives (e.g. by a first-order upwind scheme) and solve the resulting system of ordinary differential equations using ODE15S, e.g. . Look up "method of lines" for more details.

Best wishes

Torsten.

[Widitha Samrakoon]

I tried as Bill said and I get unrealistic answers. However I will try to discretize the spatial first derivatives using first order upwind method and comment about the results

[Torsten]

If you have the NAG Toolbox, use

http://www.nag.co.uk/numeric/MB/manual64_24_1/html/D03/d03plf.html

Best wishes

Torsten.