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

Hi everyone.

I am trying to solve two coupled PDEs with pdepe MATLAB. However, I have some dificcualties regarding this problem. I appericiate if anyone can help me with this matter.

So the problem I'm trying to solve is:

du1/dt = (1/rho) * du2/dy
0 = -u2 + eta * (abs(du1/dy)+eps)^(n-1)*du1/dy

As you can see there is no transient term in the second equation.

I also know that I can simply put the second equation into the first equation, however, I need to have the outputs for both u1 and u2, and that's why I am splitting it in this form.

The boundry conditions are:

u1_l = 0
u2_r = 1
u1(t=0) = 0
u2(t=0) = 0 %if anyone cares

I have coded this problem as follows:

function [c,f,s] = pdefun(x,t,u,dudx)

global rho eta n

c = [1; 0];
f = [u(2)/rho; 0];
s = [0; eta*((abs(dudx(1))+eps)^(n-1))*dudx(1) - u(2)];

end

with the following BCs and ICs:

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)

pl = [ul(1)-0; 0];
ql = [0; 1];
pr = [0; ur(2)-1];
qr = [1; 0];

end

function u0 = pdeic(x)

u0 = [0; 0];

end

And the main code for solving equation is:

global rho eta n H tf

rho = 100;
eta = 1;
n = 1.5;
H = 1;
tf = 20;

x = linspace(0,H,101);
t = linspace(0,tf,1001);

m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

By using these codes, I don't get a correct solution. Does anyone have any ideas?

Thank you,

[Mohammadamin Mahmoudabadbozchelou]

Hi Andrei,

Thanks for the response. Actually k and eta are the same thing. I forget to change it in the actual equation. Now it's correct.

I also have changed time and distance intervals, but sadly no luck!

[Bill Greene]

You say "I don't get a correct solution". Specifically, what problem are you referring to? I am assuming that the title of your post has nothing to do with this problem?

[Mohammadamin Mahmoudabadbozchelou]

At first I got the error I mentioned in the title. But after fixing some minor errors, I managed to get a solution out of it. Sorry for the confusion.

Let me ellaborate more on the specific problem. The solutions that are coming out this code are u1 and u2. As I mentioned above, one could easily put the second equation into the first equation and ask for the solution of u1 by it self. These two solutions for u1 are not the same. The second method is correct, beacuase I can check it out with other method, however, the first one gives me wrong solution. What do you think migh be the reason?

[Bill Greene]

Can you post your pdepe code for the single-equation solution that you believe is correct? I assume your posted code for the two-PDE version is the current version you are running?

[Mohammadamin Mahmoudabadbozchelou]

Yes, the posted code is two-PDE version. Here is the single equation solution:

function [c,f,s] = pdefun(x,t,u,dudx)

global rho eta n

c = 1;
f = (eta/rho)*((abs(dudx)+eps)^(n-1))*dudx;
s = 0;

end

function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)

global V n eta rho

pl = ul-0;
ql = 0;
pr = (-(V/eta)^(1/n));
qr = 1;

end

function u0 = pdeic(x)

u0 = 0;

end

clear all
close all
clc

global rho eta n H V tf

rho = 100;
eta = 1;
n = 1.5;
H = 1;
V = 1;
tf = 20;

x = linspace(0,H,101);
t = linspace(0,tf,1001);

m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

[Bill Greene]

Where did this bizarre boundary condition come from? This looks nothing like the BCs in the two-PDE version so it isn't surprising the solutions are different.

[Mohammadamin Mahmoudabadbozchelou]

From the second equation. I need to switch the second BC in the two-version PDE to a BC for u(1), right? If I want u(2) be 1 in the right boundry, I use the second equation to transfer this BC to a condistion for u(1), which is in the form I wrote in the code.