Why does pdepe throw an error depending on boundary condition parameters?

Hello,

I want to simulate 1D wave propagation by solving the (acoustic) wave equation using pdepe.

Here the mathematical problem:

And here the corresponding code:

%% Define PDE system

x = linspace(0,10,50); %descretized space

t = linspace(0,0.1,50); %descretized time

m = 0;

sol = pdepe(m,@homWaveEq,@homWaveEq_ic,@homWaveEq_bc,x,t);

p = sol(:,:,1); % acoustic pressure

%% plot time evolution

figure

surf(x,t,p)

xlabel('Distance x')

ylabel('Time t')

title('p(t)')

%% functions

function [c,f,s] = homWaveEq(x,t,u,dudx) % Equation to solve

K = 142000; % bulk modulus air

rho = 1.2; % density air

c = [1; 1];

f = [0;

K/rho*dudx(1)];

s = [u(2);

0];

end

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

function u0 = homWaveEq_ic(x) % Initial Conditions

u0 = [0; 0];

end

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

function [pl,ql,pr,qr] = homWaveEq_bc(xl,ul,xr,ur,t) % Boundary Conditions

freq = 40;

pl = [ul(1)-sin(2*pi()*freq*t); ul(2)-2*pi()*freq*cos(2*pi()*freq*t)];

ql = [0; 0];

pr = [0; 0];

qr = [1; 1];

end

The code works well for freq <= 40.

However, when freq > 40, I get the following error message:

"Error using pdepe (line 293)

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux

term involving spatial derivative."

Has anyone an idea what the reason is for that and how to avoid it? I assume the error lies in the ode15s solver.

[Torsten]

I'm surprised that pdepe can solve this equation at all. Usually, the boundary condition to the right for u1 (which gives 0=0, thus something undetermined) leads to numerical problems.

Your equation is hyperbolic in nature and pdepe is designed to solve parabolic-elliptic problems.

You should use a solver for hyperbolic problems, e.g. CLAWPACK available at

https://www.clawpack.org/

Answer 1

[Bill Greene]

Interesting.

The ODE solver(ode15s) is numerically approximating the

iteration matrix for the PDE and, due to roundoff

error, deciding this is an ODE system that it cannot solve.

One way to get around this problem is to insure that your

initiation and boundary conditions are consistent. This is

often not strictly necessary, as you have noticed, but is

usually a good idea.

If you define freq in your main program and replace your

initial condition function with

function u0 = homWaveEq_ic(x) % Initial Conditions
    if x==0
      u0=[0; -2*pi*freq];
    else
      u0 = [0; 0];
    end
  end

you can obtain solutions for larger driving frequencies.

[Paul Maurerlehner]

That makes sense. Thank you very much!