Spatial discretization has failed. (pdex1)

I am trying to use Matlabs pdex1 function to solve the below:

du/dt=d/dx(-u/2(2-u)

u0=0.75+0.25sech(x/0.5)^2

and u tends to 0.75 as x tends to infinity and minus infinity.

function pdex1

m = 0; x = linspace(0,1,200); t = linspace(0,2,50);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

u = sol(:,:,1);

surf(x,t,u)

title('Numerical solution computed with 20 mesh points.')

xlabel('Distance x')

ylabel('Time t')

figure

plot(x,u(end,:))

title('Solution at t = 2')

xlabel('Distance x')

ylabel('u(x,2)')

% -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx)

c = 1;

f = -(u/2)*(2-u);

s = 0;

% -------------------------------------------------------------- function u0 = pdex1ic(x)

u0 = 0.75+.002*sech((x/0.5)^2);

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

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

pl = .002*sech((xl/0.5)^2)-0.75;

ql = 0;

pr = .002*sech((xr/0.5)^2)-0.75;

qr = 0;

But I keep receiving the following error message and to be honest I don't really understand what it means!

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

Answer 1

[Bill Greene]

That error message is somewhat misleading.

There are two problems in your code. The first is that your boundary condition function is incorrect. From from your IC function, I believe you want

pl = ul-.002*sech((xl/0.5)^2)-0.75;
ql = 0;
pr = ur-.002*sech((xr/0.5)^2)-0.75;
qr = 0;

This change will allow you to obtain a solution; however the solution will be incorrect as you will clearly see. The error message you are getting relates to this problem. The pdepe function is designed to solve partial differential equations that are second-order in x. That is, it expects that f in your pdeex1pde function will include a term with DuDx. One work-around to this is to simply include a small DuDx term in your f-coefficient like this:

alpha = 1e-3;
f = -(u/2)*(2-u) + alpha*DuDx;

You have to experiment with the value of alpha. Too small a value and the solution contains the spurious oscillations. Too large and the solution differs significantly from the correct one.